The first experiment is a free evolution of reservoirs in the absence of input, i.e. u_i = 0. We define the network activity as A(t)=∑ixi(t)/N, with N the number of neurons in the reservoir, and A∈[0, 1]. A is also the proportion of excited neurons: A = 0 if the network is extinguished (x_i = 0 ∀i), A = 1 if the network is saturated (x_i = 1 ∀i). At the initial state, we randomly force 20% of neurons to an up state (x_i = 1), i.e. A(t = 0) = 0.2. After a transient regime of 1, 000 time steps, the reservoir reaches a steady state where we perform statistics. In the following, A will refer to the activity measured in that steady state (see Supplementary material S9 for a more formal definition). For each value of σ^⋆, we perform statistics on 100 randomly generated reservoirs (see Supplementary material S7.1 for more details on the experiment).

In the following, a bar over a variable (.)¯ represents an average over time for a given reservoir, while the brackets 〈.〉 represents an average over different randomly generated reservoirs. In the first analysis, we calculate the time-average steady activity Ā for a given reservoir and its time-variance δA2¯, where we define δA = A−Ā. We average these quantities over the reservoirs to give 〈Ā〉 and 〈δA2¯〉 for each value of σ^⋆. Next, we evaluate the average and variance over reservoirs of the binary entropy H_b, or BiEntropy (Croll, 2014) of the time-dependent activity. Compared to the Shanon entropy, the advantage of this metric is that it can discriminate ordered from disordered strings of binary digits. It has been used in machine learning (Mamun et al., 2016; Zhou and Zeng, 2022), but to our knowledge, this is the first time in reservoir computing. The binary entropy varies between 0 for fully ordered bit-streams and 1 for fully disordered ones. We compute the BiEntropy of the binarized time dependence of the steady activity for each reservoir (for the exact definition of all the metrics, see Supplementary material S10).

The time-average activity 〈Ā〉 as a function of σ^⋆ is shown in Figure 3A, for both signs of σ^⋆ (Figure 2). The green dashed line represents the value obtained for μ = 0, i.e. σ → ∞. The perfect balance in excitation (b = 0) results in half of the neurons being activated 〈Ā〉 = 0.5. The variance 〈δA2¯〉 vs. σ^⋆ is shown in Figure 3C. For the lowest values of |σ^⋆|, the reservoirs are frozen (zero variance) either extinguished (for σ^⋆ < 0) or saturated (σ^⋆> 0). This corresponds to reservoirs being respectively purely inhibitory (b = −1) or excitatory (b = 1). Already at the level of the statistics of the activity, there is a clear difference between both signs of σ^⋆: for σ^⋆ < 0 there is a threshold in σ^⋆ (vertical dashed line at σ^⋆~−0.7) above which the average activity and its variance rise abruptly and simultaneously. In contrast, for σ^⋆> 0, there is a wide region where no dynamic is detected (zero variance), yet the network is not saturated but its activity decays continuously. The variance starts rising at σ^⋆> 4 (vertical dash-dotted line).

The average BiEntropy 〈H_b〉 vs. |σ^⋆| is plotted in Figure 3B for σ^⋆ < 0 and Figure 3D for σ^⋆> 0 (left scale, black stars on both plots). These two plots zoom in the vicinity of the phase transition, as statistics are stationary elsewhere. There is a continuous transition between a fully ordered phase (〈H_b〉 = 0) and a fully disordered one (〈H_b〉 = 1). Since the BiEntropy is a measure of order, these results suggest that the transition we observe is related to the apparition of chaos in the reservoir above a critical value of σ^⋆ (Lewin and Bak, 1993; Bertschinger and Natschläger, 2004; Seifter and Reggia, 2015; Kuśmierz et al., 2020). The variance of the BiEntropy 〈δHb2〉 is shown in Figure 3B for σ^⋆ < 0 and Figure 3D for σ^⋆> 0 (right scale, orange circles on both plots). It is zero when either in the ordered or disordered phase and spikes at the transition. Its maximum coincides with 〈H_b〉≃0.5: the variance of BiEntropy captures the edge of chaos as a balance between order and disorder. More striking, it also coincides with the position at which the variance of the activity rises (vertical dashed lines). The peak of 〈δHb2〉 thus provides a clear definition of the position of two critical points: σc⋆≃-0.66 and σc⋆≃4.0, which correspond respectively to critical balances b_c≃−0.87 (94% of inhibitory synapses) and b_c≃0.19 (60% of excitatory synapses). Moreover, the transition between order and disorder is much wider for σ^⋆> 0. This asymmetry between both signs of σ^⋆ is a property of our model since θ(−x)≠±θ(x) in Eq.(1). From now on, we will refer to the critical points as the point where the maximum of BiEntropy variance is obtained, and we will define the critical regions as the regions with 〈δHb2〉≠0.

Similar to (Krauss et al., 2019b) and (Metzner and Krauss, 2022), the existence of phases separated by critical points as a parameter is varied is reminiscent of the phase diagrams drawn in thermodynamics. If we associate the state of a neuron, 0 or 1, to the state of an Ising spin, either down or up, then 〈Ā〉 corresponds to the average magnetization per spin of the network and 〈δA2¯〉 to the variance of its fluctuations, i.e., magnetization noise. At equilibrium, it is proportional to the magnetic susceptibility according to the fluctuation-dissipation theorem (Callen and Welton, 1951). The total magnetization plays the role of an order parameter, and the transition order is obtained by considering discontinuities, as a function of temperature, of the order parameter and its derivatives with respect to the external field (Landau and Lifshitz, 1980). Here we observe that the average activity is always continuous as a function of σ^⋆. At the same time, 〈δA2¯〉 is continuous for σ^⋆> 0 but shows a discontinuity at the critical point for σ^⋆ < 0. This strongly suggests that the two “phase transitions” are of a different type.

We submitted our reservoirs again to a free evolution without input [u_i(t) = 0]. For each value of σ^⋆, we created 100 reservoirs with randomly tossed weight matrices. Each reservoir is run 100 times, with a different random initial state, of activity A(t = 0) = 0.2 (for more details, see Statistics of reservoirs in Supplementary material S7).

We classify the attractor obtained in the steady-state activities, as proposed in (Krauss et al., 2019a). We categorize the activity signals into one of the four types of attractors (see Supplementary material S11 for a grounded justification of each category):

We determine the attractor obtained at the steady state for each reservoir and initial condition. We then compute the distribution of attractors for each value of σ^⋆ obtained overall the initial conditions of all reservoirs. The statistics are thus computed on 10,000 steady activities for each σ^⋆.

Figures 4A–F provide examples of attractors, encoded in the colors, for each reservoir (x-axis) and each initial condition (y-axis) for different values of σ^⋆. The left column (blue-bordered boxes) corresponds to values below the critical point (vertical blue lines on Figures 4G, H), the center column (gray-bordered boxes) to values at the critical point (vertical gray lines), and the right column (red-bordered boxes) to values above the critical points (vertical red lines). The upper row displays negative σ^⋆ values, while the lower row features positive σ^⋆ values.

Away from the critical point, a dominant color is observed, meaning that reservoirs exhibit a dominant attractor. Steady activities are predominantly extinguished and fixed on the left side of the critical point (Figures 4A, D) and irregular at the right (Figures 4C, F). Close to the critical points (Figures 4B, E), there is an increase in the diversity of attractors, as previously observed (Karimipanah et al., 2017).

Figures 4G, H show the statistical distribution of all obtained attractors vs. |σ^⋆|. As expected from the previous analysis, there is no attractor diversity on the far left and right of the plots, as we obtain one primary attractor. Dead (blue line) or fixed (green line) attractors are found for low values of |σ^⋆|, and their proportion decays slowly across the transition. Within the critical region coexist all attractors in various proportions. Chaotic attractors start to appear precisely at the transition (vertical gray lines), while the domain where cyclic attractors exist coincides with the critical region of nonzero BiEntropy variance (Figures 3B, D). The point at which cyclic attractors are most present is also precisely σc⋆. These results corroborate what we inferred in the previous section: on the disordered phase |σ⋆|> |σc⋆|, attractors are irregular, while the ordered phase is characterized by fixed or dead attractors. We note an asymmetry between both sides of the transition: irregular attractors appear only in the disordered phase. From the point of view of the attractors, both signs of σ^⋆ lead to similar behaviors, except again, that the transition region is much wider for σ^⋆> 0.

Our results suggest that the critical points enhance sensitivity to the initial states and configuration of the weights, explaining the reservoir-to-reservoir variance and increase in dynamic diversity. Reservoirs around the negative σc⋆ (Figure 4B) possessed a distribution of attractors with far more variety than the one with positive σc⋆ (Figure 4E), further reinforcing the idea that the sign of σ^⋆ produces two distinct types of critical regimes. We quantified this in the Supplementary material S12 by computing the entropy of reservoir attractor distributions plotted in Supplementary Figure 9. We interpret that result by suggesting that inhibition might be a key factor for enhancing dynamic diversity.

For the purpose of reservoir design, our findings suggest that with both critical points, most reservoirs possess an attractor obtained predominantly in most trials, independent of the initial state. The statistics of dominant reservoir attractors are presented in Supplementary Figure 9, and found to be similar to the one in Figure 3. The presence of vertical color lines in Figures 4A–F means that, in most cases, the behavior of the reservoirs does not depend on the initial state, even in the critical region (this is more thoroughly shown in the Supplementary material S12). As a consequence, a dominant attractor can be associated with each reservoir, irrespective of the initial condition.

We evaluate the performance of the networks to execute two fundamental tasks: memory and prediction. Each reservoir receives an input u(t), and the readout target is T(t) = u(t+δ), equal to the input shifted in time by δ time steps. δ < 0 corresponds to a memory task, and δ> 0 to a prediction task. For each value of σ^⋆, we use 100 reservoirs, and each reservoir is run five times, with a different random tossing of the input weight matrix (more detail on the training procedure in Supplementary material S8.4).

The first task consists of memorizing a purely random signal (i.e., uncorrelated white noise), and since there is absolutely no correlation in the input, only memorization is involved. Figure 5A illustrates this task for one value of δ = −6, with white noise as input u, and the target T.

For the second task, we explore the ability of the reservoir to predict a time series, δ = 10 time steps in the future. The input is the well-known Mackey-Glass time series, as it is a common benchmark of this type of task (Hajnal and Lörincz, 2006; Goudarzi et al., 2016; Canaday et al., 2018, among others), notably testing the ability to infer non-linear dynamics. The signal regularity is controlled by the parameter τ, see Figure 5B, ranging from periodic with τ = 5, to chaotic for τ = 28 (more information on the experiments in Supplementary material 8.3).

The performance of a reservoir is measured by computing the correlation product Corr(y, T) between the output y and the target T. A perfect match corresponds to a correlation of one, while a random output gives a zero correlation. An individual reservoir performance score is then obtained by averaging over the initial conditions. Each individual reservoir is associated with its dominant attractor, and the statistics of the performance of reservoirs are performed separately for each attractor.

The average performance is plotted as a function of |σ^⋆| in Figures 5C, E for the memory task and in Figures 5D, F for the prediction task. The left column (Figures 5C, D) corresponds to σ^⋆> 0, and the right column (Figures 5E, F) to σ < 0. The color of the lines corresponds to the attractor.

For σ^⋆ < 0 (Figures 5C, D), performance increases over a very wide range of σ^⋆, both for memory and prediction. This range includes the critical region (gray hatched area) but is vastly broader. Thus, being within the critical region is absolutely not mandatory to perform well. A shaded area in Figure 5 indicates the spreading of the results. There is none in plots Figures 5C, D, meaning that all reservoirs perform exactly the same for a given σ^⋆. Moreover, as σ^⋆ is increased through the critical region, the dominant attractors change [see zooms in plot Figures 5C, D], but surprisingly, there is no discontinuity in the performance. Indeed, inside the gray area, even though the four attractor categories are present, their respective performance all align. This strongly suggests that the dynamics of the reservoir, as measured in the absence of input, is irrelevant for the performance. Only the value of σ^⋆ matters. In both tasks, the average performance decreases monotonically with increasing difficulty via τ and δ. In the memory task, we register a dip in performance close to the critical point. This goes against the common assumption that the edge of chaos is optimal for memory (Natschläger et al., 2005). In the prediction task, the peak of performance roughly coincides with the critical region, except for the greater difficulty, where the peak is slightly on the left.

The picture is very different for σ^⋆> 0 (Figures 5E, F). First, the region in which some level of performance is observed is comparable to the critical region observed in free-running reservoirs. Second, there is substantial variability in performance across different reservoirs, as indicated by the large shaded areas. Despite this variability, there is an overall dependence of performance on σ^⋆. The average performance of distinct dominant attractor categories is much noisier. However, despite being more noisy, the average performance of the distinct attractor categories aligns again, so there is still no evidence that the attractor category has any significant impact on the reservoir's performance. We observe that performance decreases as the difficulty of the memorization task increases, but interestingly, this trend appears to be inverted in the prediction task.

Once again, the two signs of σ^⋆ give rise to different behavior. In particular, networks with σ^⋆> 0 memorize better and are less reliable than those with σ^⋆ < 0 but have poorer prediction capability. Yet, in all cases, attractors do not seem to be correlated to performance, as the top performance can be found in any of the four attractor categories.

Our results somewhat challenge the common assumption that the edge of chaos is optimal for performance and suggest that this is true for reservoirs with a majority of excitation but not necessarily with a majority of inhibition. Reservoirs with negative σ^⋆ exhibit very reliable performances with very low reservoir-to-reservoir variability over a range in σ^⋆ much broader than the critical region. Since reservoirs behave the same, in practice, it is sufficient to generate one, with σ^⋆ at the left of the critical region. However, if the goal is optimal memorization, it is wiser to choose σ~0.4 in the critical region and try different reservoirs until finding a good one, which requires training and testing.

We analyze the absolute value of the performance of all reservoirs independently of the control parameter. For each reservoir, we study its performance in the memory task as a function of its performance in the prediction tasks (see Section 4.1 methodology). For this, we fixed values of δ (memory) and τ (prediction), and we chose three levels of difficulty: simple (τ = 5, δ = −2), average (τ = 20, δ = −6) and difficult (τ = 28, δ = −10).

Figure 6 displays the performance of reservoirs as colored dots. As in the previous section, each color corresponds to the dominant attractor of the reservoir. The vertical axis represents the performance in the prediction task (Mackey-Glass), and the horizontal axis that of the memory task (white noise). The columns correspond to the sign of σ^⋆, negative on the left, positive on the right. The rows correspond to the degree of difficulty, from simple (top) to difficult (bottom).

For σ^⋆ < 0 (left column), there is an apparent relationship in performance between the two tasks. The degree of difficulty roughly acts as a scaling factor on the curves, but the correlation is always clear. The reservoirs which are good at predicting do memorize well as well. There is, however, a reentrant region of irregular reservoirs (red) which are the best at memorizing but are not optimal for predicting, in particular for the intermediate degree of difficulty (see the red loop on Figure 6B). This point corresponds to the maximum on the right of the dip observed in Figure 5C. Interestingly, in all difficulties, reservoirs at the critical point (encircled dots) create a narrow area with a good overall performance. This picture refines our previous analysis of the impact of attractors on performance, as it seems that extinguished reservoirs can perform slightly better at memory than the others, and the same for the irregular reservoirs at prediction. Nonetheless, the difference between attractors is rather small, and it could be argued that is insignificant.

The picture is very different for σ^⋆> 0 (right column). Performances are distributed as clouds of points. In Figure 5, we observe a large reservoir variability in both tasks. Some reservoirs are suitable for one task and bad for the other one. For intermediate and difficult tasks, some reservoirs outperform the best ones with σ^⋆ < 0 in both memorization and prediction. Reservoirs at the critical point are found on the right of the plots, i.e. they promise good memorization but are nonetheless widely spread, especially in prediction (vertical axis). Overall, the distinct attractors occupy the space in overlapping and indistinguishable clouds; this confirms the previous analysis that attractors do not play a role in performance.

Correlations in performances are very different for both signs of σ^⋆. Choosing a reservoir with σ^⋆~−0.66 ensures a good, general-purpose reservoir but with suboptimal performance. In contrast, going into the positive side of σ^⋆ may lead to the best reservoirs in a given task or even better general purpose reservoirs, but this comes at a price: those gems cannot be found by the statistical analysis we have performed on their free running activity.