Frequency-Based Reservoir Computing Arthur S Powanwe∗ Department of Mathematics, Western University, London, ON, Canada (Dated: July 28, 2026) Reservoir computing has emerged as an efficient machine learning framework for predicting time series generated by dynamical systems. In contrast to other machine and deep learning approaches, a reservoir computer trains only the output layer via linear regression, leaving the reservoir (recurrent layer) untrained. This simplification makes reservoir computers easier to train and more amenable to experimentation. However, because current reservoirs consist of networks of randomly connected nodes and require the optimization of numerous hyperparameters, a framework that precisely explains how reservoir computing operates and how it can be optimized remains missing. Here, we propose a frequency-based reservoir inspired by the brain’s oscillatory dynamics and its hierarchy of timescales. The frequency-based reservoir can be interpreted as an ensemble of inde-
arXiv:2607.24420v1 [stat.ML] 27 Jul 2026
pendent oscillatory units, each processing a portion of the input’s frequency content. This allows us to understand the reservoir’s internal behavior by modeling it as a single unit driven by an external input. Borrowing from the theory of a nonlinear oscillator forced by complex periodic inputs, we found that units of the frequency-based reservoir selectively amplify and store specific input frequencies, which are then used for prediction. The frequency-based reservoir performs as well as or better than equivalent random reservoirs. Furthermore, the frequency-based approach can be optimized to improve short-term prediction, a property that random reservoirs lack. Finally, we show that the frequency-based reservoir can also predict complex spatiotemporal dynamics. Our results show that reservoir computing can be designed using brain properties and theoretical insights borrowed from the physics of forced nonlinear oscillators.
2 I.
INTRODUCTION
Reservoir computers are a class of neural networks that consist of three main modules: an input layer that maps the input to the reservoir; a recurrent layer (the reservoir) that stores information about the input; and an output layer that maps the reservoir to a target [1–3]. The input and recurrent layers are initialized randomly and kept fixed, whereas only the output layer is trained, typically using a simple linear regression. This provides a significant advantage for reservoir computers over other machine learning models that require training the recurrent layer. Despite their simple training, reservoir computers have been shown to perform surprisingly well at forecasting time series generated by chaotic dynamical systems [4, 5]. Understanding how the reservoir encodes information about the input to generate a prediction is crucial for designing optimal reservoir computers. Unfortunately, apart from initializing the recurrent layer at the edge of chaos, there is currently no methodology for designing optimal reservoir computers. Building reservoir computers that clarify how the recurrent layer processes the input to generate a prediction will not only enhance interpretability but also enable effective implementation.
Recent works have focused on making reservoir computers more interpretable by incorporating brain-like properties. This includes building biologically motivated reservoirs that contain excitatory and inhibitory neurons [6], with connections between nodes subject to synaptic plasticity. It was demonstrated that controlling the excitatoryinhibitory (E-I) balance in a brain-inspired reservoir could significantly boost its performance [7]. Similarly, increasing the number of connections between neurons has been shown to improve accuracy. Recently, oscillations have been introduced as a critical substrate for computations and learning [8–10]. Oscillations have been observed in the brains of several species and across diverse frequency bands [11]. They usually exhibit a substantial variability, appearing as short, random epochs of synchrony [12, 13]. They are believed to be involved in many cognitive processes [14], including communication [15–17], attention[18], and working memory [19, 20]. However, it remains poorly understood how neural oscillations contribute to computations and learning. Recent reservoir computing approaches that use oscillations have been shown to outperform traditional random reservoirs [8] and to be more interpretable [10]. Nevertheless, those reservoirs still employ a random recurrent layer, which makes them difficult to optimize. Here, we introduce a reservoir computing approach based on the diversity of oscillatory behavior across frequency bands and brain regions. The reservoir comprises excitatory and inhibitory populations, each oscillating at a distinct frequency. It uses the frequency property of units for computation; hence the name ”frequency-based reservoir computer”. In addition, the frequency-based reservoir computer is endowed with the hierarchy of timescales [21–27], an important and widespread property of brain computation.
The frequency-based reservoir computing approach is motivated by the oscillatory dynamics of the Wilson-Cowan model [28, 29]. Each unit consists of two nodes, each representing a population of excitatory or inhibitory neurons. The model was designed to ensure that each unit oscillates at a specific frequency, behaves as a brain region with its own timescale, and exhibits the noisy behavior inherent in brain recordings. In contrast to previous works [8–10], the coupling between the units is not chosen randomly. Furthermore, we found that coupling between the units was optional, suggesting that computations happen at the unit level. With this choice, the frequency-based reservoir can be viewed as an ensemble of small, independent reservoirs (each comprising two nodes) that process a portion of the
3 input’s frequency content. Borrowing from the theory of a nonlinear oscillator forced by a complex input, we found that the reservoir units selectively amplify frequency components of the input that are near their intrinsic frequencies. A property that is lacking in random reservoirs. We further found that frequency-based reservoirs perform as well as, or better than, equivalent random reservoirs and can be optimized to improve short-term prediction performance. Finally, we extend the predictive capabilities of frequency-based reservoirs to complex spatiotemporal chaotic inputs. Taken together, our results show that brain-inspired reservoir computing can be designed to be fully interpretable and optimized for a specific task.
II.
RESULTS
Frequency-based reservoir computing model
We consider the dynamics of the reservoir computer obeying the following equation
r(t+1) = (1 − α)r(t) + α tanh Ar(t) + Win u(t) .
(1)
The matrix A describes the connectivity among the nodes and is usually chosen to be a random matrix (Fig.1-(a)). The matrix Win maps the lower-dimensional input u(t) to the higher-dimensional recurrent layer and is also a random matrix. The coefficient α is a leaky rate that controls how quickly the variable r(t) incorporates the reservoir’s memory of past values. The goal in reservoir computing is to train an output layer Wout so that the dynamics of the reservoir follow a target signal y(t). In this work, we are interested in the prediction of the next step y(t)=u(t+1) of the input u(t). During training, we minimize the following quantity.
∥Wout r(t) − y(t)∥2
(2)
This is done using a linear regression. During testing, we use the target as input, that is we add a feedback loop (Fig.1-(a)) so that the input is now Wout r(t) and the reservoir becomes a self-evolving dynamical system that projects to the target via Wout . Random reservoirs (Fig.1-(a)) have been shown to perform very well in chaotic timeseries forecasting [2]. However, the reservoirs are often regarded as a black boxes, lacking interpretability and being difficult to optimize. The magic relies on hyperparameter optimization, a process that is also poorly understood and requires substantial computational resources rather than a transparent methodology. Here, we propose a reservoir computing framework that does not rely on random connections among nodes, but rather on the intrinsic oscillatory dynamics of each unit (Fig.1). Our reservoir framework is motivated by the Wilson-Cowan model, a well-known model for representing the dynamics of coupled populations of excitatory and inhibitory neurons in the brain [28, 29]. A unit describes a system of two connected nodes, one node representing a population of excitatory neurons (red disc in Fig.1-(b)) and the other a population of inhibitory neurons (blue disc in Fig.1-(b)). We are interested in the oscillatory behavior of each unit in isolation and when coupled to other units. In the original Wilson-Cowan model (Supplementary Information), the E and I populations have a recurrent coefficient that captures interactions
4
FIG. 1. Two reservoir computing approaches with the same architecture but distinct recurrent layers. (a) In a random reservoir, the recurrent layer is initialized by randomly choosing couplings among the nodes. (b) In a frequency-based reservoir, the recurrent layer consists of coupled excitatory and inhibitory nodes that can be coupled or not. Each unit has intrinsic oscillatory dynamics with a specific frequency used for computation. For both reservoir types, the matrix Win that maps the input to the recurrent layer is chosen randomly, and the output matrix Wout that maps the recurrent layer dynamics to the target is obtained through linear regression.
within each population. They receive constant external input, and coupling between units occurs only from excitatory populations. Two distinct types of coupling were chosen, a ring-graph coupling where each unit is only coupled to its nearest neighbor (see dashed red line Fig.1-(b) here dashed means optional) and a distance-dependent graph coupling where the strength of the coupling between units decreases as the distance between them increases. However, we found that the coupling mainly shaped each unit’s intrinsic dynamics. Inspired by the original Wilson-Cowan dynamics, we proposed a version that follows the common reservoir-computing dynamics (Equation 1). Because computation happens at the unit level (see Supplementary Information), the coupling between the units is optional (dashed red line Fig.1-(b)). We do not claim that coupling does not matter; it mainly shapes the intrinsic behavior of each local unit. In addition, for each coupled network, we can find a corresponding uncoupled one with similar properties.
In the frequency-based reservoir, the recurrent matrix A is chosen so that its eigenvalues are complex conjugate pairs (Fig.2. Similarly, the eigenvalues of the Jacobian matrix J = (1 − α)I + αA of Equation.1 are also complex conjugate pairs (Fig.2). The recurrent matrix is block-diagonal, with each block representing the intrinsic dynamics of a unit.
5
FIG. 2. Spectra of the recurrent A and Jacobian J matrices for different values of the coefficient α The y axis shows the imaginary part and the x axis the real part of the eigenvalues.
ν1
ω1
0
0 ···
0
0 ···
0
0
−ω1 0 0 . .. A= 0 0 . .. 0 0
ν1
0
0 ···
0
0 ···
0
0
ν2
ω2 · · ·
0
0 ···
0
0 −ω2 ν2 · · · .. .. .. . . . . . .
0 .. .
0 ··· .. . ···
0 .. .
νk
0 0 0 .. . 0 0 .. . ωN νN
0
0
0 ···
ωk · · ·
0
0 .. .
0 .. .
0 · · · −ωk νk · · · .. . .. . . . . · · · .. .
0 .. .
0
0
0 ···
0
0 ···
νN
0
0
0 ···
0
0 · · · −ωN
(3)
The coupling between the units is optional; therefore, we did not include it here. However, in the original WilsonCowan formulation, we considered distinct coupling matrices. Our recurrent matrix is not randomly generated; it can be selected appropriately based on the input we aim to forecast. We arranged the units so that the first unit has the lowest angular frequency, ωmin , and the last unit, N, has the largest angular frequency, ωmax . With this choice, each unit is a small reservoir of two connected nodes. The frequency-based reservoir is an ensemble of N small reservoirs (each having two nodes). Furthermore, the reservoir is endowed with a hierarchy of timescales [21–27] (units with large indices k have larger angular frequencies ωk and appear early in the hierarchy. In contrast, those with small k have smaller angular frequencies ωk and appear later), a key property of brain computation. Finally, we added white Gaussian noise to the dynamics to model the stochasticity of neural recordings (see Methods). To investigate the ability of the frequency-based reservoir computing to predict chaotic time series, we chose two dynamical systems known to generate chaotic dynamics, namely the Mackey-Glass [30] and the Lorenz-63 chaotic model [31]. We evaluated
6
FIG. 3. Closed-loop forecasting of the Mackey-Glass and the Lorenz-63 timeseries and the corresponding attractors. (a) Short-term prediction of the Mackey-Glass timeserie using the random and (b) the frequency-based reservoirs. (c) Long-term prediction of the Mackey-Glass attractor for both the random and frequency-based reservoirs.(d) Short-term prediction of the Lorenz-63 timeserie from the random and (e) Frequency-based reservoirs. (f ) Long-term prediction of the Lorenz-63 attractor for both the random and frequency-based reservoirs. For all panels, the red data refer to the real data, while the blue and green data refer to the predictions from the random and frequency-based reservoirs, respectively.
random and frequency-based reservoirs under these two inputs and trained them to predict their next steps. This enables us to obtain the output matrix used during closed-loop testing to forecast future values. To maintain a fair comparison between the random and frequency-based reservoirs, we used the same Hyperparameters for both. The results are shown in Fig.3, which displays only the data from the closed-loop forecasting phase. We found that the frequency-based reservoir performed as well as, or better than, the random reservoirs under the same conditions. For the example in Fig.3, we did not aim to optimize both reservoirs, and the hyperparameters were selected manually. The goal was to compare random and frequency-based reservoirs under identical conditions. We consider a network size of N=301 units (601 nodes) for the frequency-based reservoir. Note that for the first unit with an imaginary part equal to zero, we only chose one node instead of two to form a unit. The values of the component of the recurrent matrix A are shown in Fig.2, and the coefficient α was set to be 0.5. The frequency-based reservoir is a promising candidate for time-series forecasting. By design, it has a very low number of components. For a matrix of size 2Nx2N, the number of non-zero components is 4N, and the number of distinct non-zero components is 2N. In addition, each unit is independent of the others, and computation occurs at the unit level. This stands in clear contrast to a random reservoir, in which computations arise from the collective behavior of the nodes in the recurrent layer. The frequency-based reservoir, therefore, offers a unique opportunity to understand how the reservoir operates and which computations occur within the recurrent layer. We can clearly understand the computations inside the recurrent layer by studying the dynamics of a single unit driven by the input.
7 Selective-frequency amplification of multi-frequency periodic inputs as a mechanism for information storage
To understand how the Frequency-based reservoir encodes and stores information about the input, let’s study the dynamics of a single unit driven by a one-dimensional input u(t). We write the dynamics of the Frequency-based reservoir for a single unit as follows r1 (t+1) = (1 − α)r1 (t) + α tanh νr1 (t) + ωr2 (t) + ϵ1 u(t)
(4)
r2 (t+1) = (1 − α)r2 (t) + α tanh − ωr1 (t) + νr2 (t) + ϵ2 u(t)
(5)
To have a better understanding of how a single unit processes an arbitrary input u(t), we first consider a periodic input with multiple frequencies defined as
u(t) =
L X
Al sin(ωl t)
(6)
l=1
The single unit of the frequency-based reservoir is a nonlinear oscillator that oscillates at a frequency given by the imaginary part of the corresponding Jacobian. The Equation. 4-5 and Equation.6 represent the dynamics of a nonlinear oscillator forced by a complex periodic input. Understanding the computation at a single unit is a physics problem borrowed from the theory of forced nonlinear oscillators. We investigate this question by analyzing the power spectrum of the system of Equation.4-5 as the imaginary part of A varies according to Fig.2. We are interested in the log(P ower) as a function of the frequency and the imaginary part ω. We first consider a periodic input described in Equation. 6. The input is made of L=5 fundamentals frequencies fl = 1, 2, 3, 4, 5Hz so that ωl = 2πfl . The amplitude of each fundamental frequency is given by Al and is all scaled by a factor ϵ. The corresponding Power spectrum is given in Fig.4-(a) and shows the contribution of each frequency to the input u(t). The 4 Hz frequency component contributes the most, followed by the 5 and 2 Hz frequencies, then the 1 Hz frequency, and finally, the 3 Hz frequency is very weak due to the low value of A3 . For low amplitude values ϵ = 0.05, a single unit responds to all the strong frequency components in the input independently of the value of its own intrinsic frequency (related to the imaginary part ω). However, a unit powerfully amplifies the input’s frequency components near its inherent frequency. The mechanism is similar to nonlinear resonance but involves multiple frequency components. Therefore, a unit oscillating at a given frequency will store more information about frequency components of the input that are nearby. We refer to this mechanism as selective frequency amplification. For intermediate amplitude values ϵ = 0.1, the selective frequency amplification persists. Still, the unit receives more information overall about all frequency components, including the element with the lowest contribution, A3 . For a higher amplitude value ϵ = 0.5, the selective amplification mechanism is destroyed, and the intrinsic oscillator frequency being near a frequency component of the input does not lead to additional amplification. The unit responds preferentially to the element with the highest amplitude. Therefore, components with higher amplitude are more amplified. In addition, we observe the emergence of higher-order harmonics due to nonlinear effects, resulting in increased amplitudes. We consider a second periodic signal generated by the Mackey-Glass dynamical system with τ = 16. The corresponding power spectrum Fig.4-(e) showed that the periodic Mackey-Glass can be seen as a periodic input with multiple fundamental frequencies of distinct amplitude, as in the first case Equation.6. By also scaling the Mackey-Glass by a factor ϵ, the analogy with
8
FIG. 4. Selective frequency amplification of multi-frequency periodic inputs. (a) Power spectrum of the multifrequency periodic input described in Equation.6. (b) log(Power) of a single unit, precisely Equation.4 for low amplitude values ϵ = 0.05. (b) log(Power) for intermediate amplitude values ϵ = 0.1. (c) log(Power) for large amplitude values ϵ = 0.1. (e) Power spectrum of the periodic Mackey-Glass (τ = 16) dynamical system. (f ) log(Power) of a single unit with the MackeyGlass as input for low amplitude values ϵ = 0.05.(g) log(Power) for intermediate amplitude values ϵ = 0.5. (h) log(Power) for large amplitude values ϵ = 0.5. For all the panels concerned (b),(c),(d),(f ),(g) and (h) the coefficients ν and ω were varied as in Fig.2.
Equation.6 is complete. For low amplitude values ϵ = 0.05, the selective amplification is observed as previously described for the first periodic input. For intermediate amplitude values Fig.4-(g) selective frequency amplification persists and for high amplitude values Fig.4-(h) selective frequency amplification is destroyed. Selective amplification occurs for periodic inputs with low to intermediate amplitudes and is destroyed when the amplitude becomes large. Can selective frequency amplification be the mechanism at work in the frequency-based reservoir? To answer this question, we need to clarify a few points: 1) The mechanism we have described is for a periodic system that shows a Power spectrum with a discrete and finite number of frequencies. Chaotic systems, on the other hand, exhibit a continuous spectrum in which an infinite number of frequencies contribute to the dynamics. However, even in that case, some chaotic systems still exhibit frequencies that dominate their dynamics. It is then possible that the same mechanism generalizes for chaotic systems.
Selective frequency amplification generalizes to a one-dimensional chaotic system
Our goal is to understand how the frequency-based and random reservoirs process external inputs. We consider the chaotic Mackey-Glass input used in Fig.3-(a-c). The delay in the Mackey-Glass equation introduces additional dimensions, rendering it an effective infinite-dimensional system. However, we refer to the dimensions as the number of
9
FIG. 5. Selective frequency amplification of the chaotic Mackey-Glass input. (a) Power spectrum of the chaotic Mackey-Glass (τ = 17) dynamical system. (b) log(Power) of the frequency-based reservoir. (c) log(Power) of the random reservoir.(d) Power spectra of the first node in the frequency-based (green) and random (blue) reservoirs. (e) Power spectra of node 153 in the frequency-based (green) and random (blue) reservoirs. (f ) Power spectra of node 197 in the frequency-based (green) and random (blue) reservoirs.
variables simulated to generate the input. By this definition, the Mackey-Glass system is one-dimensional. In contrast to the periodic case (τ = 16), the power spectrum of the chaotic case (τ = 17) is continuous, meaning that an infinite number of frequencies contribute to the dynamics of the Mackey-Glass. Nevertheless, a finite number of frequency components contribute the most as shown in Fig.5-(a). We first consider the frequency-based reservoir. Each unit is uncoupled. Precisely, the dynamics of each unit in the recurrent layer are given by Equation.4-5 with ν,ω replaced by νk and ωk . In addition, the inputs amplitude ϵ1,2 in Equation.4-5 are replaced by the corresponding coefficients Win (k). A total of 301 (N=601 nodes) of these units are considered with their parameters νk and ωk covering all the spectrum shown in Fig.2. The nodes in the reservoir were arranged so that the first node corresponds to the lowest intrinsic frequency (lowest value of ωk ) and the last node to the largest intrinsic frequency (most significant value of ωk ). This arrangement facilitated interpretation and aligned with the single-unit case. We compute the power spectrum of each node, and we display the log(power) as a heat map in Fig.5-(b). We found that nodes in the reservoir selectively amplify frequencies near their intrinsic frequencies, a mechanism we have identified as selective frequency amplification in the context of multifrequency periodic inputs. For the chaotic input here, the selective frequency amplification generalizes to the continuous power spectrum of the Mackey-Glass input (Fig.5-(b)). We further investigate the power spectra of individual nodes. The first node has the lowest intrinsic frequency and amplifies the dominant frequency component, the lowest component of the input, while the fastest components are dampened (Fig.5-(d)). Note that the Mackey-Glass equation has a strong frequency component that is always amplified independently of the node. Node
10
FIG. 6. Selective frequency amplification of a multi-dimensional chaotic system (a) Power spectra of the x,y and z components of the chaotic Lorenz-63 dynamical system.(b) log(Power) heat map of the frequency-based reservoir.(c) log(Power) heat map of the random reservoir. (d) Power spectra of node 4 in the frequency-based (green) and the random (blue) reservoirs. (e) Power spectra of node 152 in the frequency-based (green) and random (blue) reservoirs. (f ) Power spectra of node 155 in the frequency-based (green) and the random (blue) reservoirs.
153 has an intermediate intrinsic frequency and amplifies intermediate- to high-frequency components of the input, but not low-frequency components. Finally, node 197 exhibits a high intrinsic frequency and selectively amplifies high-frequency components of the input while attenuating lower-frequency components. For the random reservoir (Fig.5-(c)), there is no identifiable mechanism. All nodes process the input similarly. In contrast to the frequencybased reservoir, all the nodes preferentially process the largest frequency components present in the input. This is clearly evident in the power spectra of nodes 1, 153, and 197, which all exhibit a strong value at the same frequency.
Selective frequency amplification generalizes to a multi-dimensional chaotic system
We now focus on the case of the three-dimensional Lorenz-63 chaotic system used as the input to the frequency-based and random recurrent layers of the reservoir computers ( Fig.3-(d,e,f )). The power spectra of the x,y variables are similar and concentrated at lower frequencies. They are broadband, with no clear dominant component. In contrast, the z component shows a strong peak at intermediate frequency values, characteristic of a dominant periodicity (Fig.6-(a)). Overall, the Lorenz-63 dynamical system exhibits low, intermediate (dominant), and large frequencies. We analyze how the frequency-based and random reservoirs process the Lorenz-63 dynamical system. As in the one-dimensional input case, each unit in the frequency-based reservoir is decoupled from other units. Their dynamics are described by the Equation.4-5 with ν,ω replaced by νk , ωk covering the spectrum described in Fig.2. We found
11 that the inputs to each odd and even component of a unit were W1in (1, :, k)x + W1in (:, 2, :, k)y + W1in (:, 3, k)z and W2in (1, :, k)x + W2in (:, 2, :, k)y + W2in (:, 3, k)z respectively. Moreover, the x and y components are similar, and the effective input contribution were approximated to (W1in (1, :, k) + W1in (:, 2, :, k))(x+y)W1in (:, 3, k)z and (W2in (1, :, k) + W2in (:, 2, :, k))(x+y) + W2in (:, 3, k)z respectively. This allows us to determine the amount of input from the x-, y, and z-components that each unit receives. The situation, nevertheless, is similar to the multi-frequency input in Equation 6. The mostly lower frequency components of the x and y variables will be scaled by the amplitude (W1in (1, :, k) + W1in (:, 2, :, k)) and (W2in (1, :, k) + W2in (:, 2, :, k)) while the frequency components of the z variable will be scaled by W1in (:, 3, k) and W2in (:, 3, k). We found that the units in the reservoir selectively amplified frequency components close to their intrinsic frequencies while dampening other frequencies (depending on their amplitudes). The frequency-based reservoir employs the selective frequency-amplification mechanism (Fig.6-(c)). Precisely, the node 4 in the Frequency-based reservoir selectively amplifies lower frequency components ( mostly x and y variables) of the Lorenz-63 input (see Fig.6-(d)) and dampens higher ones. This can be explained by the fact that the frequency of the fourth node is low, the amplitudes to the x and y variables are strong enough, and the amplitude to the z variable is weak. Nodes (152 and 155) with intermediate or higher intrinsic frequencies selectively amplify the strong frequency component of the z variable and the higher frequency components of the three variables. Therefore, the selective frequency amplification mechanism generalizes to multidimensional chaotic systems. In the random recurrent layer, we found no discernible structure or transparent mechanism for processing input frequency. As in the one-dimensional case, all nodes responded similarly to the three-dimensional input.
Selective frequency amplification as a mechanism for chaotic time series prediction
The frequency-based reservoir processes temporal information differently from the random reservoir via a selective frequency-amplification mechanism. The question becomes how the reservoir takes advantage of this mechanism to make an accurate prediction of chaotic dynamics, as we showed in Fig.3. In the reservoir computing framework, only the readout layer, Wout , is trained using linear regression. The information stored in the reservoir is crucial for accurate predictions. Previous attempts to explain how the reservoir stores information about the input have failed because of its random nature. The frequency-based reservoir does not have this curse. Each unit in the reservoir stores a portion of the input information in the frequency domain. Specifically, units with lower intrinsic frequencies will store lower-frequency components of the input, units with intermediate intrinsic frequencies will store intermediatefrequency components, and units with higher intrinsic frequencies will store higher-frequency components. In addition, strong frequency components of the input are stored in all the units. The frequency-based recurrent layer acts as a set of many small reservoirs. Each small reservoir is a unit (a layer of two nodes) that stores part of the input information. This is similar to an ensemble of small reservoirs, each processing a part of the input in the frequency domain. When the reservoir is sufficiently large and the stimulation duration is sufficiently long, all information about the input is stored in the recurrent-layer units. A simple linear combination Wout r(t) of all the units should be enough to recover all the information about the input and therefore make an accurate prediction. In addition, the frequency-based reservoir is endowed with a hierarchy of timescales, an essential property for brain computation. This makes the selective frequency amplification an excellent mechanism for brain computation.
12
FIG. 7. Optimizing the frequency-based reservoir for short time prediction performance of the Lorenz-63. (a) Short time prediction of the Lorenz evaluated in Lyapunov units (1 Lyapunov unit =0.9). (b) Spectrum of the default connectivity matrix A (black) and the corresponding Jacobian J (green)for α = 0.5 and N=150.(c) Spectrum of the optimized connectivity matrix A (black) and the corresponding Jacobian J (green) for α = 0.5. (d) Log(power) heat map of the frequencybased reservoir for the default connectivity in (a). (e) Log(power) heat map of the frequency-based reservoir for the optimized connectivity in (c). (f ) Mean and standard deviation over 500 input initial conditions. The green bar represents both the default and optimized frequency-based reservoir connectivity, whereas the blue represents the random reservoir.
Optimizing the Frequency-based reservoir
The main issue with random reservoirs is the difficulty of optimizing their performance for a specific task or target. The frequency-based reservoir selectively stores information about the input in each of its units. Choosing the intrinsic coefficients α, νk , and ωk for each unit is key to achieving frequency-based reservoir performance. So far, we have chosen a default configuration in which the units have intrinsic frequencies ranging from ω1 = 0 to the largest value, ωN . This configuration ensures that the network is provided with all the necessary intrinsic frequencies to store all the frequency components of any given input. However, this configuration may not be optimal for specific inputs. For example, some inputs may have frequency components concentrated in a lower band, and the reservoir may not require (or require only a few) units with higher intrinsic frequencies. On the other hand, some inputs may have a broader frequency band, requiring a reservoir with more units that have high intrinsic frequencies. Designing the reservoir so that its predictions can track the chaotic trajectory for a sufficiently long period is particularly challenging, especially with limited training data. This is known as a short-term prediction task (Fig.7-(a)). We investigate the performance of the frequency-based reservoir on a short-time prediction task and compare it with a random reservoir. We use the chaotic Lorenz-63 system as the input. We set the reservoir size to N=301 nodes and the training set size to 1500 steps only. We found that under the default connectivity (Fig.7-(b)), the frequency-based reservoir
13
FIG. 8. Forecasting of spatiotemporal chaotic dynamical systems by a frequency-based reservoir computing. (a) real data obtained after simulating the KS spatiotemporal chaotic systems and using it as input to the frequency-based reservoir. (b) Prediction from the frequency-based reservoir during closed-loop testing.(c) Error between the real data (a) and the prediction (b). (d) Long-term prediction of the KS chaotic attractor (the third x3 (y-axis) and first x1 (x-axis) components of the x variable are used to build the attractor) by the frequency-based reservoir (green).
performs worse than the random reservoir (Fig.6-(f )). The default connectivity leads to an overrepresentation of units with low intrinsic frequencies and an overprocessing of low-frequency components of the input (Fig.6). This may explain the poor performance of the frequency-based reservoir under the default configuration. To address the overrepresentation of low intrinsic frequencies in the default configuration, we reshape the connectivity so that half of the units with low intrinsic frequencies now have the same higher intrinsic frequencies as the remaining units. This new configuration allows the reservoir to efficiently process high-frequency components of the input (Fig.7-(e)) and significantly overperforms the random reservoir (Fig.7-(f )). This simple optimization technique demonstrates that, unlike a random reservoir, a frequency-based reservoir can be optimized to predict a specific input.
The Frequency-based reservoir can forecast spatiotemporal chaotic systems
We found that the frequency-based reservoir performs very well across a broad class of low-dimensional chaotic systems (see the supplementary information). We are now interested in how the frequency-based reservoir performs on more complex systems, such as spatiotemporal chaos. We used the Kuramoto-Shivanski as input. The results of the closed-loop forecasting are shown in Fig.8. The frequency-based reservoir can predict the spatiotemporal chaotic trajectory of the Kuramoto-Shivanski (Fig.3-(b)) for a sufficiently longer epoch of time (Fig.3-(c)) as well as its corresponding attractor (Fig.3-(d)). Our results demonstrate the efficacy of frequency-based reservoir computing for predicting complex spatiotemporal data.
14 Properties of the frequency-based reservoir
The performance of random reservoirs depends on the optimization of several hyperparameters, including the input amplitude, the reservoir spectral radius, and the leaky rate α. This is usually done using several techniques, such as random search or Bayesian optimization [32]. Hyperparameter optimization requires substantial computational resources and relies on no practical methodology or physical principles beyond minimizing the error between the ground truth and the prediction. The process also contributes to the black box nature of random reservoirs. In a frequency-based reservoir, the hyperparameters have physical meanings and can be selected based on the theory of forced nonlinear oscillators. The leaky rate α shapes the eigen spectrum of the reservoir (see Fig.2) and therefore participates in the selective frequency amplification mechanism. When α = 1, the reservoir is an echo state network, and the eigen spectrum is directly given by the recurrent matrix A. The specific choice of the reservoir components ωk and νk is set by the value of the spectral radius ρ. Specifically, the spectral radius determines the maximum frequency of the reservoir. Higher spectral radius values correspond to higher maximum frequencies, whereas lower values correspond to lower maximum frequencies (see Methods). When α is decreased, the eigen spectrum is shrunk, leading to an overall reduced eigen spectrum (see Fig.2). Therefore, the leaky rate α and the spectral radius both account for the reservoir’s intrinsic frequencies. The network input amplitude is given by the values of the input matrix Win , which is scaled by a coefficient ϵ (see Methods). We found that lower and intermediate values of ϵ yielded accurate predictions, whereas higher values of ϵ yielded inaccurate predictions. The inaccuracy of predictions for large values of ϵ is consistent with the destruction of the selective frequency amplification mechanism (see Fig.4). Thus, the three hyperparameters can be interpreted in the context of the theory of forced nonlinear oscillator, and their choice can be motivated by that same theory. This is in clear contrast to a random reservoir, where it’s difficult to interpret the hyperparameters. A recent work [32] has shown that noise is an additional hyperparameter that can significantly enhance the performance of a random reservoir. We also found that the frequency-based reservoir’s performance was significantly noise-dependent. As in the random case [32], we identified an optimal noise level for each task that yielded the best performance. We explain this resonant behavior by the role of noise in sustaining neural rhythms. In fact, noise is crucial for neural rhythms observed in vivo and modeled by linearly stable dynamical systems [12, 33]. Adding noise to the frequency-based reservoir allows fast rhythms that would otherwise vanish [34, 35] to contribute to input processing. In contrast, lower rhythms that would have remained coherent will be perturbed by such noise [35]. Thus, the role of noise in the frequency-based reservoir can also be interpreted through the lens of neural rhythm theory. The frequency-based reservoir provides an alternative approach in which the hyperparameter optimization step can be replaced by a transparent methodology based on oscillatory theory.
III.
DISCUSSION
Many systems in physics [36], ecology[37–39], biology[40–42], and neuroscience[11] exhibit complex oscillatory dynamics. Harnessing the oscillatory behavior of systems to perform practical computations and enable predictions has recently attracted attention. This goes from using the oscillatory behavior of Kuramoto-like dynamics for diverse computations [43, 44] to optimizing networks of coupled harmonic oscillators for classification tasks [9, 10]. In this paper, we make an additional contribution to these existing results. We have proposed a reservoir computing based
15 on the oscillatory behavior of populations of coupled excitatory and inhibitory neurons. In contrast to random reservoirs, in which nodes are connected at random, the frequency-based reservoir comprises excitatory-inhibitory units that are not necessarily coupled. In frequency-based reservoirs, computations occur at the unit level, and couplings among units are optional. The key insight of frequency-based reservoirs is that each unit oscillates at a specific intrinsic frequency and timescale, selectively amplifying the frequency components of the input near that frequency. By increasing the unit’s intrinsic frequency from low to high, each unit will store a portion of the input information during stimulation. If the number of units and the duration of the stimulation are sufficient, a simple linear combination should be enough to recover all the information contained in the input and therefore reconstruct it accurately.
The frequency-based reservoir introduced here offers several advantages in terms of biological realism. First, neural rhythms recorded across multiple brain regions and frequency bands are believed to be involved in many cognitive processes, including communication via coherence [15, 16], attention [18], memory, and learning [14]. Second, it has been shown that the brain possesses a hierarchy of time scales [21–26], in which early sensory regions process input rapidly. In contrast, other regions interact with inputs on much slower timescales. Like the brain, our frequency-based reservoir exhibits intrinsic oscillatory dynamics across different frequency bands and processes inputs according to a hierarchy of timescales. Additionally, the Wilson-Cowan model is a well-established model for brain dynamics, used to simulate oscillatory neural recordings across various species and frequency bands. The selective frequency amplification mechanism identified here may reflect a biologically plausible way the brain processes and stores information to support accurate future predictions. It is biologically plausible that the coordinated activity of neurons (here reflected in the oscillatory dynamics of each Wilson-Cowan unit) in each brain region plays an essential role in processing and storing information about sensory inputs. Moreover, the frequency-based reservoir provides a transparent interpretation of what is happening inside the reservoir, a property that is invaluable for efficient optimization.
Random reservoirs have demonstrated strong performance in predicting chaotic time series [5] with a minimal number of learnable parameters. However, they suffer from the curse of interpretability since the reservoir is seen as a black box. Previous works have sought to shed light on the black-box nature of random reservoirs by designing more interpretable and biologically realistic reservoirs [7] or by connecting random reservoirs to nonlinear autoregressive processes [45]. The frequency-based reservoir leverages the intrinsic oscillatory behavior and the hierarchy of timescales observed across brain areas. Our approach is grounded in the theory of nonlinear oscillators driven by complex, multifrequency periodic inputs. We showed that the dynamics of the frequency-based reservoir could be interpreted as those of a single forced nonlinear oscillator [46]. By varying the oscillator intrinsic frequency and timescale, we can infer which input frequency components are processed and stored at a given node/unit of the reservoir. This provides interpretability for the frequency-based reservoir, which, to our knowledge, has not been observed in prior works. The frequency-based reservoir has a small number of connections. For a reservoir with N units, the frequency-based reservoir has 2N-1 nodes and 4N-1 connections. It is equivalent to an ensemble of small, independent reservoirs (each with two nodes) that process a portion of the same input. The absence of couplings among the units in the reservoir makes it readily amenable to experimental realization and to applications in neuromorphic computing. This enables the design and fabrication of reservoirs with oscillatory units based on the selective frequency-amplification mechanism described here.
16
A recent line of works has used the harmonic oscillator as a basic unit for computations. Results from networks of coupled harmonic oscillators are promising [47, 48] and indicate that they can outperform neural networks that lack intrinsic oscillatory dynamics. Harmonic oscillatory recurrent networks (HORN) have been shown to outperform non-oscillatory networks in pattern recognition tasks [9]. The HORN is based on a large number of randomly selected parameters and uses a random connectivity matrix. It was trained using backpropagation through time. In contrast to our frequency-based reservoir, the complexity of HORN and its training procedure may obscure the exact computational mechanisms at play. A reservoir computing approach based on networks of coupled random oscillators has been proposed in [8]. The architecture is similar to the HORN architecture with a random connectivity matrix. As in our case, only the readout matrix is trained. However, their reservoir is significantly different from our frequency-based reservoir. Our reservoir is fully interpretable, does not rely on the random choice of intrinsic parameters or a random connectivity matrix. In contrast with the reservoir in [8], the dynamics of each unit can be fully understood independently of those of other units. We have proposed a mechanism that our frequency-based reservoir uses for computation and prediction. Another recent work [10] has used reservoir computing based on the damped harmonic oscillator for brain-rhythm classification and handwritten and spoken digit classification. They proposed a resonant mechanism that their reservoir uses for classification. Although their resonant behavior shares some similarities with selective frequency amplification, our mechanism is more general. The selective frequency amplification operates by amplifying a spectrum of continuous frequencies near the oscillator’s intrinsic frequency. It extends beyond periodic inputs and generalizes to high-dimensional chaotic dynamics. Moreover, our frequency-based reservoir can be designed and optimized to predict a specific input. A property that has not been shown in previous oscillatory-based neural networks.
The dynamics of the damped harmonic oscillator and the Wilson-Cowan model can overlap in the parameter regime in which the latter has complex-conjugate eigenvalues with negative real parts. However, the two models are different. A single Wilson-Cowan neuron is a nonlinear dynamical system, and its response to inputs can be highly complex. A single damped harmonic oscillator is a linear dynamical system. The nonlinearity is introduced through the coupling with other units [8, 9]. While a single Wilson-Cowan unit is already nonlinear, in most damped-harmonic-oscillator neural networks, the nonlinearity is emergent. Thus, computations in our frequency-based recurrent network happen at the unit level, and coupling between the units is optional. This opens interesting questions about computations with the frequency-based reservoir. For example, what is the minimum number of units required to predict a given input? Which input frequency components should be accurately processed by the frequency-based reservoir to increase the performance in short/long-term prediction? The answer to the first question can be anticipated for periodic inputs as 2N ≤ 2p, where N is the number of units and p the number of harmonics of the input. For chaotic inputs, the question remains open. For the second question, the optimization procedure (Fig.7) suggests that the frequency-based reservoir should accurately process the fast components of the chaotic inputs to improve short-term prediction.
Another distinction between a frequency-based reservoir and a neural network based on a damped harmonic oscillator concerns the role of noise. While noise was optional [10] or seen as a robustness measure[9], in coupled damped harmonic oscillator networks, we found that noise was critical for accurate short and long-term prediction in the
17 frequency-based neural network. The effect of noise was observed for chaotic inputs with fast (e.g., the Mackey-Glass) and low-frequency components. Previous work has identified a noise-induced resonance mechanism in random reservoirs [32]. Here, we argue that noise is critical because it induces oscillations in units with a fast intrinsic frequency and perturbs rhythms in units with a low inherent frequency [35]. In fact, units with large intrinsic frequencies decay quickly, and noise acts as a random force that sustains them [12, 13, 34]. The noise-induced rhythms are filtered noise signals suitable for predicting the fast components of the inputs. On the other hand, units with small intrinsic frequencies remain coherent for a very long epoch of time, similar to self-sustained or limit-cycle oscillations. Such coherent, self-sustaining rhythms could be less flexible in efficiently processing and storing external inputs. Noise may act as a force that breaks the self-sustained or coherent nature of such low-frequency oscillations [35], enabling them to process and store information about external inputs efficiently.
IV.
METHODS
Frequency-based recurrent layer
The construction of the recurrent matrix A in Equation 1 was motivated by the oscillatory dynamics of a network of connected Wilson-Cowan units (see Supplementary Information). The dynamics of a single Wilson-Cowan unit have been extensively studied [28, 29]. In the regime in which its complex-conjugate eigenvalues have negative real parts, the dynamics decay to zero. However, adding noise sustains the rhythms [34, 35, 49]. This is especially relevant for fast rhythms. We found that noise was essential for predicting chaotic time series. For that reason, the reservoir’s input was corrupted by noise. This leads to the following effective dynamics for the frequency-based reservoir.
r(t+1) = (1 − α)r(t) + α tanh Ar(t) + Win u(t) + σξ(t) ,
(7)
where ξ(t) is a time-dependent Gaussian white noise with zero mean and unit variance, the coefficient σ controls its intensity. The matrix A was initialized by first choosing the spectral radius ρ, and then selecting N angular frequencies (or imaginary parts) ωk equally spaced from a minimum value 0 ≤ ωmin to a maximum value ωmax ≤ ρ. p The eigenvalue real parts were then deduced as νk = ρ2 − ωk2 . Furthermore, we arranged the units so that the first unit has the lowest angular frequency, ωmin , and the last unit, N, has the largest angular frequency, ωmax . With this choice, each unit is a small reservoir of two connected nodes and a spectral radius of ρ. The frequency-based reservoir is an ensemble of N small reservoirs (each having two nodes) with the same spectral radius, ρ. Furthermore, the frequency-based reservoir is endowed with a hierarchy of timescales (units with large indices k have larger angular frequencies ωk and appear earlier in the hierarchy). In contrast, those with small k have smaller angular frequencies ωk and appear later ), a key property of brain computation [21–26]. The Jacobian of the Equation. 7 is obtained by linearization, leading to J = (1 − α)I + αA,
(8)
where I is an 2Nx2N identity matrix. The corresponding eigenvalues are complex conjugate pairs given as (1 + α) +
18 ανk ± jαωk , with k=1,...,N and j =
√
−1 the pure imaginary number.
The oscillatory frequency of each unit is obtained as
fk =
αωk . 2π
(9)
The coefficient α scales the value of the intrinsic frequency. When α = 1, we have an echo state and A and J coincide. In contrast to a random reservoir, where the role of α is speculative and often associated with reservoir memory, here we can clearly identify α as a factor that scales the network’s intrinsic frequencies, thereby contributing to the selective frequency-amplification mechanism. The matrix Win represents a feedforward layer that maps the noise-corrupted input u(t) + σξ(t) to the recurrent layer. Its size is (d, 2N-1) where d is the dimension of the input and 2N-1 that of the reservoir. Precisely, it was chosen as Win = ϵ(2*U(d, 2N-1)-1 ), where U is a (d, 2N-1) matrix with components generated from a uniform distribution between 0 and 1. This makes the values of Win random numbers uniformly generated between −ϵ and ϵ. Win is the only matrix that was randomly generated. With this formulation, the effective noise intensity is ϵσ scaled by the corresponding random value generated by the matrix 2*(U(d,2N-1)-1 ). The output matrix Wout was obtained by minimizing the quantity ∥Wout r(t) − y(t)∥2 . Previous work in reservoir computing has employed Tikhonov regularization to mitigate overfitting. Here, we impose a different constraint. In addition to minimizing the quantity ∥Wout r(t) − y(t)∥2 we consider the solution Wout that also minimizes the norm ∥Wout ∥2 . In Matlab 2024a, we used the command lsqminnorm to obtain Wout . During testing, we introduce a feedback loop by replacing the input u(t) with Wout r(t). The dynamics of the reservoir during closed-loop testing are a self-evolving discrete dynamical system given by
r(t+1) = (1 − α)r(t) + α tanh Ar(t) + Win Wout r(t) .
(10)
During training and closed-loop testing, the even nodes of the reservoir were further squared. This is a common practice when the reservoir uses the tanh nonlinearity. We applied this approach to all inputs in the main text.
Forecasting the Mackey-Glass chaotic system
The input in Fig.3-(a-c) was generated by the Mackey-Glass dynamical system [30] whose dynamic is given as
du βu(t − τ ) = − γu(t), dt 1 + uk (t − τ )
(11)
with the following coefficients β = 0.2,γ = 0.1, k = 10 and τ = 17. We used a recurrent layer with N = 301 units, which corresponds to 2N-1=601 units. The input was simulated using a fourth-order Runge-Kutta scheme with a step size of dt = 0.01. A total of T = 5.103 + 5.103 + 105 samples were generated. The first 5.103 were discarded as transient, the next 5.103 were used as input during the reservoir washup period, and the remaining 105 were used for training. The self-evolving reservoir during the closed-loop testing period was simulated for T = 2.105 steps to
19 construct the attractor in Fig.8-(c). Further, the noise intensity was chosen as σ = 10−3.5 , α was set to 0.5 and the input strength ϵ = 0.05. The same input configuration was used for the random reservoir. The random reservoir was generated as a random matrix with symmetrically distributed coefficients and degree q. The value of q was chosen as q = 0.00525 to ensure that the random and frequency-based reservoirs have approximately the same number of connections.
Forecasting the Lorenz-63 chaotic system
The Lorenz-63 chaotic input was generated using the following dynamical system [31]
dx = 10(y(t) − x(t)) dt dy = x(t)(28 − z(t)) − y(t) dt dz 8 = x(t)y(t) − z(t) dt 3
(12) (13) (14)
and the corresponding input was constructed as u = [x(t), y(t), x(t)]. In Fig.8-(d-f ) we used the Lorenz-63 as input. The recurrent layer size was the same as in the Mackey-Glass case. The Lorenz-63 was simulated with a Fourth-order Runge-Kutta algorithm. All the other training parameters were almost the same as in the Mackey-Glass case, apart from the noise strength σ = 10−7 and the training length T = 2.105 . For Fig.7, we significantly reduced the training size to T = 1500 steps, the network size was also reduced to N = 301, and the noise intensity was modified to σ = 10−6 . The same parameters were chosen for the random recurrent layer, with the degree q = 0.0067.
Forecasting the spatiotemporal Kuramoto-Sivashinki chaotic system
In Fig.8, the input we chose was the Kuramoto-Sivashinsky. A one-dimensional model of instabilities in flames is described as ∂u ∂2u ∂4u ∂u =− 2 − −u ∂t ∂x ∂x4 ∂x
(15)
The variable u = u(x, t) describes a spatiotemporal behavior. We simulate on the spatial dimension Equation.15 x x from 0 to L = 22. We chose sinusoidal initial condition u(x, 0) = sin 2π and periodic boundary condition. L 1 The temporal and spatial resolutions were set to dt = 0.25 and dx = , respectively. The reservoir size was set to 3 N = 2201 and initialized with the default configuration given in Fig.2. The input was corrupted by Gaussian white noise with strength σ = 10−1 . The other parameters were the same as in the Mackey-Glass and Lorenz-63 cases. The training phase consists of a total of T = 2.105 steps. In the closed loop testing, we self-evolve the reservoir for T = 104 steps to generate the attractor shown in Fig.8-(d).
20 Power spectra
All the power spectra in the main text and the supplementary information were generated using the pwelch command in Matlab 2024a. The heat maps were generated by taking the logarithm of the power spectra for all nodes in the recurrent layer.
The Wilson-Cowan frequency-based reservoir
We first developed our analysis and obtained our results by using the Wilson-Cowan model. The discrete dynamics of the Wilson-Cowan recurrent layer are given by
Ek (t + 1) = (1 − ∆tαe )Ek (t) + ∆tβe 1 − Ek (t) f (sEk ) Ik (t + 1) = (1 − ∆tαik )Ik (t) + ∆tβi 1 − Ik (t) f (sIk )
(16) (17)
The total synaptic inputs to the E and I populations are given by
k k sEk (t) = Wee Ek − Wei Ik + hke +
N X
Wkl El (t) +
N X l=1
kp Win up (t) + σek ξek (t)
(18)
kp Win up (t) + σik ξik (t) .
(19)
p=1
l=1 k sIk (t) = Wie Ek − Wiik Ik + hki +
P X
Wkl El (t) +
P X p=1
k k k Wiik Wei ,Wie αe ,αik , βe ,βi , Each unit (Ek , Ik ) has its own oscillatory dynamics described by the coefficients Wee
,hke ,hki ,σek and σik . The biological meanings of these parameters are provided in the supplementary material. The parameter ∆t is a step-size that plays the same role as α. The nonlinear response function f is the sigmoid function 1 . The 13 parameters make the single-unit Wilson-Cowan model difficult to modify and optimize f (x) = 1 + e−x k for those unfamiliar with it. However, we found that the parameters αik and Wee were useful for modifying the intrinsic frequency of each unit (see Supplementary Information). The units are coupled through the matrix W and receive the input u(t) through the feedforward input layer Win . The coupling targets both nodes of each unit, but only the E node projects to other units. In the supplementary information, we performed a detailed analysis of the Wilson-Cowan model described here. We chose two coupling matrices, a ring and a distance-dependent one. For both matrices, the main effect was to shape the unit intrinsic frequencies. This led us to conclude that computation occurs at the unit level and that coupling is an additional parameter that contributes to prediction by shaping each unit’s intrinsic dynamics. In fact, we found that even a network of uncoupled Wilson-Cowan units (W=0) could make accurate predictions. However, each Wilson-Cowan unit has 13 coefficients, making it difficult to optimize for predictive performance. However, the recurrent layer described in the main text is easier to control and optimize than the Wilson-Cowan model. To obtain the output matrix Wout , we collect all the reservoir states as
21
E2k−1 , k = 1, 2, 3, ...N, r(t) = I2k
(20)
and minimize the quantity Equation.2. In contrast to the tanh nonlinearity, the sigmoid nonlinearity does not require squaring the even reservoir state components. During closed-loop testing, the input for the next step is replaced by the prediction from the current step. The Wilson-Cowan reservoir, as described here, implements the selective amplification mechanism to store information about the inputs and uses that information for prediction, as described in the main text (see Supplementary information).
[1] Herbert Jaeger, “The “echo state” approach to analysing and training recurrent neural networks-with an erratum note,” Bonn, Germany: German national research center for information technology gmd technical report 148, 13 (2001). [2] Herbert Jaeger and Harald Haas, “Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication,” science 304, 78–80 (2004). [3] Herbert Jaeger, Mantas Lukoševičius, Dan Popovici, and Udo Siewert, “Optimization and applications of echo state networks with leaky-integrator neurons,” Neural networks 20, 335–352 (2007). [4] Mantas Lukoševičius and Herbert Jaeger, “Reservoir computing approaches to recurrent neural network training,” Computer science review 3, 127–149 (2009). [5] Gouhei Tanaka, Toshiyuki Yamane, Jean Benoit Héroux, Ryosho Nakane, Naoki Kanazawa, Seiji Takeda, Hidetoshi Numata, Daiju Nakano, and Akira Hirose, “Recent advances in physical reservoir computing: A review,” Neural Networks 115, 100–123 (2019). [6] Wolfgang Maass, Thomas Natschläger, and Henry Markram, “Real-time computing without stable states: A new framework for neural computation based on perturbations,” Neural computation 14, 2531–2560 (2002). [7] Keshav Srinivasan, Dietmar Plenz, and Michelle Girvan, “Boosting reservoir computing with brain-inspired adaptive control of ei balance,” Nature Communications 16, 10212 (2025). [8] Andrea Ceni, Andrea Cossu, Maximilian W Stölzle, Jingyue Liu, Cosimo Della Santina, Davide Bacciu, and Claudio Gallicchio, “Random oscillators network for time series processing,” in International Conference on Artificial Intelligence and Statistics (PMLR, 2024) pp. 4807–4815. [9] Felix Effenberger, Pedro Carvalho, Igor Dubinin, and Wolf Singer, “The functional role of oscillatory dynamics in neocortical circuits: a computational perspective,” Proceedings of the National Academy of Sciences 122, e2412830122 (2025). [10] Mark A Kramer, “Brain-inspired interpretable reservoir computing with resonant recurrent neural networks,” arXiv preprint arXiv:2506.17083 (2025). [11] Gyorgy Buzsaki, Rhythms of the Brain (Oxford University Press, 2006). [12] Dajun Xing, Yutai Shen, Samuel Burns, Chun-I Yeh, Robert Shapley, and Wu Li, “Stochastic generation of gamma-band activity in primary visual cortex of awake and anesthetized monkeys,” Journal of Neuroscience 32, 13873–13880a (2012). [13] Samuel P Burns, Dajun Xing, and Robert M Shapley, “Is gamma-band activity in the local field potential of v1 cortex a “clock” or filtered noise?” Journal of Neuroscience 31, 9658–9664 (2011). [14] Xiao-Jing Wang, “Neurophysiological and computational principles of cortical rhythms in cognition,” Physiological reviews 90, 1195–1268 (2010). [15] Pascal Fries, “A mechanism for cognitive dynamics: neuronal communication through neuronal coherence,” Trends in
22 cognitive sciences 9, 474–480 (2005). [16] Pascal Fries, “Rhythms for cognition: communication through coherence,” Neuron 88, 220–235 (2015). [17] Agostina Palmigiano, Theo Geisel, Fred Wolf, and Demian Battaglia, “Flexible information routing by transient synchrony,” Nature neuroscience 20, 1014 (2017). [18] Conrado A Bosman, Jan-Mathijs Schoffelen, Nicolas Brunet, Robert Oostenveld, Andre M Bastos, Thilo Womelsdorf, Birthe Rubehn, Thomas Stieglitz, Peter De Weerd, and Pascal Fries, “Attentional stimulus selection through selective synchronization between monkey visual areas,” Neuron 75, 875–888 (2012). [19] Jun Yamamoto, Junghyup Suh, Daigo Takeuchi, and Susumu Tonegawa, “Successful execution of working memory linked to synchronized high-frequency gamma oscillations,” Cell 157, 845–857 (2014). [20] Mikael Lundqvist, Jonas Rose, Pawel Herman, Scott L Brincat, Timothy J Buschman, and Earl K Miller, “Gamma and beta bursts underlie working memory,” Neuron 90, 152–164 (2016). [21] Uri Hasson, Eunice Yang, Ignacio Vallines, David J Heeger, and Nava Rubin, “A hierarchy of temporal receptive windows in human cortex,” Journal of neuroscience 28, 2539–2550 (2008). [22] John D Murray, Alberto Bernacchia, David J Freedman, Ranulfo Romo, Jonathan D Wallis, Xinying Cai, Camillo PadoaSchioppa, Tatiana Pasternak, Hyojung Seo, Daeyeol Lee, et al., “A hierarchy of intrinsic timescales across primate cortex,” Nature neuroscience 17, 1661–1663 (2014). [23] Caroline A Runyan, Eugenio Piasini, Stefano Panzeri, and Christopher D Harvey, “Distinct timescales of population coding across cortex,” Nature 548, 92–96 (2017). [24] Mehran Spitmaan, Hyojung Seo, Daeyeol Lee, and Alireza Soltani, “Multiple timescales of neural dynamics and integration of task-relevant signals across cortex,” Proceedings of the National Academy of Sciences 117, 22522–22531 (2020). [25] Joshua H Siegle, Xiaoxuan Jia, Séverine Durand, Sam Gale, Corbett Bennett, Nile Graddis, Greggory Heller, Tamina K Ramirez, Hannah Choi, Jennifer A Luviano, et al., “Survey of spiking in the mouse visual system reveals functional hierarchy,” Nature 592, 86–92 (2021). [26] Ana MG Manea, Anna Zilverstand, Kamil Ugurbil, Sarah R Heilbronner, and Jan Zimmermann, “Intrinsic timescales as an organizational principle of neural processing across the whole rhesus macaque brain,” Elife 11, e75540 (2022). [27] Songting Li and Xiao-Jing Wang, “Hierarchical timescales in the neocortex: Mathematical mechanism and biological insights,” Proceedings of the National Academy of Sciences 119, e2110274119 (2022). [28] Hugh R Wilson and Jack D Cowan, “Excitatory and inhibitory interactions in localized populations of model neurons,” Biophysical journal 12, 1–24 (1972). [29] Edward Wallace, Marc Benayoun, Wim Van Drongelen, and Jack D Cowan, “Emergent oscillations in networks of stochastic spiking neurons,” Plos one 6 (2011). [30] Michael C Mackey and Leon Glass, “Oscillation and chaos in physiological control systems,” Science 197, 287–289 (1977). [31] Edward N Lorenz, “Deterministic nonperiodic flow 1,” in Universality in Chaos, 2nd edition (Routledge, 2017) pp. 367–378. [32] Zheng-Meng Zhai, Ling-Wei Kong, and Ying-Cheng Lai, “Emergence of a resonance in machine learning,” Physical Review Research 5, 033127 (2023). [33] Georgios Spyropoulos, Matteo Saponati, Jarrod Robert Dowdall, Marieke Louise Schölvinck, Conrado Arturo Bosman, Bruss Lima, Alina Peter, Irene Onorato, Johanna Klon-Lipok, Rasmus Roese, et al., “Spontaneous variability in gamma dynamics described by a damped harmonic oscillator driven by noise,” Nature communications 13, 1–18 (2022). [34] Arthur S Powanwe and André Longtin, “Determinants of brain rhythm burst statistics,” Scientific Reports 9, 1–23 (2019). [35] Arthur S Powanwe and André Longtin, “Amplitude-phase description of stochastic neural oscillators across the hopf bifurcation,” Physical Review Research 3, 033040 (2021). [36] Stephen Strogatz, “Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering (studies in nonlinearity),” (2001).
23 [37] Alan J McKane, James D Nagy, Timothy J Newman, and Marianne O Stefanini, “Amplified biochemical oscillations in cellular systems,” Journal of Statistical Physics 128, 165–191 (2007). [38] Ganna Rozhnova and Ana Nunes, “Fluctuations and oscillations in a simple epidemic model,” Physical Review E 79, 041922 (2009). [39] Alan J McKane and Timothy J Newman, “Predator-prey cycles from resonant amplification of demographic stochasticity,” Physical review letters 94, 218102 (2005). [40] Albert Goldbeter, “Modelling biochemical oscillations and cellular rhythms,” Current Science , 933–939 (1997). [41] Albert Goldbeter, Biochemical oscillations and cellular rhythms: the molecular bases of periodic and chaotic behaviour (Cambridge university press, 1997). [42] Leon Glass, “Synchronization and rhythmic processes in physiology,” Nature 410, 277–284 (2001). [43] Gabriel B Benigno, Roberto C Budzinski, Zachary W Davis, John H Reynolds, and Lyle Muller, “Waves traveling over a map of visual space can ignite short-term predictions of sensory input,” Nature Communications 14, 3409 (2023). [44] Luisa HB Liboni, Roberto C Budzinski, Alexandra N Busch, Sindy Löwe, Thomas A Keller, Max Welling, and Lyle E Muller, “Image segmentation with traveling waves in an exactly solvable recurrent neural network,” Proceedings of the National Academy of Sciences 122, e2321319121 (2025). [45] Daniel J Gauthier, Erik Bollt, Aaron Griffith, and Wendson AS Barbosa, “Next generation reservoir computing,” Nature communications 12, 5564 (2021). [46] Benjamin Lindner, Kai Dierkes, and Frank Jülicher, “Local exponents of nonlinear compression in periodically driven noisy oscillators,” Physical review letters 103, 250601 (2009). [47] T Konstantin Rusch and Siddhartha Mishra, “Coupled oscillatory recurrent neural network (cornn): An accurate and (gradient) stable architecture for learning long time dependencies,” arXiv preprint arXiv:2010.00951 (2020). [48] T Konstantin Rusch, Ben Chamberlain, James Rowbottom, Siddhartha Mishra, and Michael Bronstein, “Graph-coupled oscillator networks,” in International Conference on Machine Learning (PMLR, 2022) pp. 18888–18909. [49] Arthur S Powanwe and André Longtin, “Brain rhythm bursts are enhanced by multiplicative noise,” Chaos: An Interdisciplinary Journal of Nonlinear Science 31 (2021). [50] Andreas Daffertshofer and Bernadette van Wijk, “On the influence of amplitude on the connectivity between phases,” Frontiers in neuroinformatics 5, 6 (2011). [51] Priscilla E Greenwood, Mark D McDonnell, and Lawrence M Ward, “Dynamics of gamma bursts in local field potentials,” Neural computation 27, 74–103 (2015). [52] Steven H Strogatz, “Nonlinear dynamics and chaos: with applications to physics,” Biology, Chemistry and Engineering , 1 (1994). [53] Crispin W Gardiner et al., Handbook of stochastic methods, Vol. 3 (springer Berlin, 1985). [54] Priscilla E Greenwood, Mark D McDonnell, and Lawrence M Ward, “A kuramoto coupling of quasi-cycle oscillators with application to neural networks,” Journal of Coupled Systems and Multiscale Dynamics 4, 1–13 (2016). [55] Arthur S Powanwe and André Longtin, “Phase dynamics of delay-coupled quasi-cycles with application to brain rhythms,” Physical Review Research 2, 043067 (2020). [56] Roberto C Budzinski, Alexandra N Busch, Samuel Mestern, Erwan Martin, Luisa HB Liboni, Federico W Pasini, Ján Mináč, Todd Coleman, Wataru Inoue, and Lyle E Muller, “An exact mathematical description of computation with transient spatiotemporal dynamics in a complex-valued neural network,” Communications Physics 7, 239 (2024). [57] Nurani Rajagopal Rohan, C Vigneswaran, Sayan Ghosh, Kishore Rajendran, A Gaurav, and V Srinivasa Chakravarthy, “Deep oscillatory neural network,” Scientific Reports 15, 40968 (2025). [58] T Anderson Keller and Max Welling, “Neural wave machines: learning spatiotemporally structured representations with locally coupled oscillatory recurrent neural networks,” in International Conference on Machine Learning (PMLR, 2023)
24 pp. 16168–16189. [59] Naama Barkai and Stanislas Leibler, “Circadian clocks limited by noise,” Nature 403, 267–268 (2000).