In the following, we are analyzing networks composed of N_neu = 100 deterministic neurons with arctangent activation functions. The random matrix of connection weights is set up in a controlled way, so that the density d of non-zero connections as well the balance b between excitatory and inhibitory connections can be pre-defined independently (for details see Section 4). Visualizations of typical weight matrices for different combinations of the statistical control parameters d and b are shown in Figure 1.

We first investigate free-running networks without external input and compute a dynamical phase diagram C_ss(b, d) of the average correlation C_ss = C(s_t, s_t+1) between subsequent system states (Figure 2a, for details see Section 4). The resulting landscape is mirror-symmetric with respect to the line b = 0, due to the symmetric activation functions of our model neurons, combined with the definition of the balance parameter. Apart from the region of very low connection densities with d ≤ 0.1, the phase space consists of three major parts: the oscillatory regime in networks with predominantly inhibitory connections (b ≪ 0, left green area in Figure 2a), the chaotic regime with approximately balanced connections (b ≈ 0, central blue and red area in Figure 2a), and the fixed point regime with predominantly excitatory connections (b ≫ 0, right green area in Figure 2a).

It is important to note that C(s_t, s_t+1) is a root-mean-square (RMS) average over all the N_neu × N_neu pairwise correlations between subsequent neural activations (so that negative and positive correlations are not distinguished), and that these pairwise correlations are properly normalized in the sense of a Pearson coefficient (each ranging between –1 and +1 before the RMS is computed). For this reason, C(s_t, s_t+1) is close to one (green) both in the oscillatory and in the fixed point regimes, where the system is behaving regularly. By contrast, C(s_t, s_t+1) is close to zero (blue) in the high-density part of the chaotic regime, where the time-evolution of the system is extremely irregular. Medium-level correlations (red) are therefore expected in the transition region between these two extreme dynamical regimes, and they are indeed found in the correlation phase diagram for densities larger than ≈0.3 in the form of narrow stripes at the border of the chaotic “valley.” It is however surprising that medium-level correlations also exist across the whole chaotic valley for relatively low densities d ∈ [0.1, 0.3]. Since medium-level correlations are thought to be optimally suited for information processing (Bertschinger and Natschläger, 2004; Natschläger et al., 2005; Legenstein and Maass, 2007; Schrauwen et al., 2009; Büsing et al., 2010; Toyoizumi and Abbott, 2011; Dambre et al., 2012), it is remarkable that this can take place not only at the classical “edge of chaos” (between the oscillatory and the chaotic regime), but also in other (and less investigated) regions of the network's dynamical phase space.

In order to verify the nature of the three major dynamical regimes, we investigate the time evolution of the neural activations for selected combinations of the control parameters b and d. In particular, we fix the connection density to d = 0.5 and gradually increase the balance from b = −0.5 to b = +0.5 in five steps (Figures 2b–f). As expected, we find almost perfect oscillations (here with a period of two time steps) for b = −0.5 (case Figure 2b), at least after the transient period in which the system is still carrying a memory of the random initialization of the neural activations. At b = 0 (case Figure 2d), we find completely irregular, chaotic behavior, and at b = +0.5 (case Figure 2d) almost all neurons reach the same fixed point. However, the cases close to the two edges of the chaotic regime reveal an interesting intermediate dynamic behavior: For b = −0.25 (cases Figures 2c,h), most neurons are synchronized in their oscillations, but some are out of phase. Others show a long-period regular “beating”-like behavior superposed on the oscillations of period two (see the longer time trace in Figure 2h). For b = +0.25 (cases Figures 2e,i), most neurons reach (approximately) a shared fixed point, but some end up in a different, individual fixed point, thus resembling a state of quenched disorder. However, a sub-group of neurons is simultaneously engaged in long-period oscillations (see the longer time trace in Figure 2i).

The apparent irregularity of the neural activations in case Figure 2d does not necessarily imply chaotic behavior. To demonstrate the sensitive dependence of the neural trajectories on the initial condition, we change the activation of only a single neuron at t = 0 by a small amount of 0.1 and re-run the simulation. We find that drastic, system-spanning differences appear between the two time evolutions after about 30 time steps (see Figure 2j).

Moreover, we observe that the memory time τ of the system for the information imprinted by the initialization (that is, the duration of the transient phase) depends systematically on the balance parameter: Deep within the oscillatory regime (b = −0.5, case Figure 2b), τ is short. As we approach the chaotic regime (b = −0.25, case Figure 2c), τ increases, finally becoming “infinitely” long at b = 0 (case Figure 2d). Indeed, from this viewpoint the chaotic dynamics may be interpreted as the continuation of the transient phase. As we move deeper into the fixed point regime (cases Figures 2e,f), the memory time τ is decreasing again.

In the medium and high-density regime of the phase diagram, we find for negative values of the bias parameter mainly oscillations of period two, as the large number of negative weights causes the neurons to switch the sign of their sigmoidal outputs from one time step to the next. However, in the low-density regime, the magnitude of the neuron's total input is reduced and we then find also oscillations with larger periods (case Figure 2g).

Next, we feed into the network a relatively weak external input (with a coupling strength of η = 0.5), consisting of independent normally distributed random signals that are continuously injected to each of the neurons (for details see Section 4).

We find that the external input destroys the medium-level state-to-state correlations C(s_t, s_t+1) in most parts of the chaotic regime, except at the classical edge of chaos (Figure 3a, red). Moreover, the input also brings the state-to-state correlations in the fixed point regime down to a very small value, as now the external random signals are superimposed onto the fixed points of the neurons.

Another important practical factor is the ability of neural networks to store information, i.e., to take up external information at any point in time and to incorporate it into their system state. We quantify this ability of information import by the RMS-averaged correlation C(x_t, s_t+1) between momentary input and subsequent system state. Surprisingly, we find that information import is best, i.e., C(x_t, s_t+1) is large, in the low-density part of the chaotic regime, including the lowest part of the classical edge of chaos (region between chaotic and oscillatory regimes), but also at the opposite border between the chaotic and fixed point regimes (Figure 3b, green and red). We thus come to the conclusion that (at least for weak external inputs with η = 0.5) our network model is simultaneously capable of information import and information processing only in the low-density part of the classical edge of chaos.

To backup this unexpected finding, we also quantify information storage and information import by the average pair-wise state-to-state mutual information I(s_t, s_t+1) (Figure 3c), and the mutual information between the momentary input and the subsequent system state I(x_t, s_t+1) (Figure 3d), respectively. These mutual-information-based measures can also capture possible non-linear dependencies, but are computationally much more demanding (for details see Section 4).

Despite of these drastic differences between the two measures, we obtain practically the same phase diagrams for information import and information storage/processing when we use the RMS-averaged pairwise correlations (Figures 3a,b) and when we use the mutual information (Figures 3c,d). This congruence may simply indicate the absence of higher-order statistical dependencies between subsequent states in our specific RNN system. However, in the context of adaptive stochastic resonance, we have already reported a surprisingly close relation between linear correlation and mutual information for a large range of model systems (Krauss et al., 2017). Taken together, these findings suggest a possible link between information-processing and dynamical approaches to complexity science (Mediano et al., 2021).

Furthermore, we compare the results to a computationally more tractable approximation of the mean pairwise mutual information, where only a sub-population of 10 neurons is included to the evaluation. It also shows the same basic characteristics (Figures 3e,f), implicating the possibility to approximate mutual information in large dynamical systems, where an exhaustive sampling of all joint probabilities necessary to calculate entropy and mutual information is impractical or impossible.

In order to test the robustness of the above results on information import, we re-compute the phase diagram of the correlations between the input and a later system state (Figure 4), now however varying some of the parameters that have been kept at their standard values (w = 0.5, N = 100, Δt = 1, η = 0.5) so far. We obtain results similar to Figure 3b when the fluctuation width w of the Gaussian weight distribution is increased to w = 1 (Figure 4a), when the number of neurons in reduced to N = 50 (Figure 4b), and when the lag-time between input signal and system state is increased to Δt = 2 (Figure 4c). However, when the fluctuation width of the weight distribution is reduced to w = 0.5, which decreases the total neural inputs and therefore brings the system closer to the linear regime, we find that now both edges of chaos become available for information uptake (Figure 4d).

We return to our standard parameters (w = 0.5, N = 100, Δt = 1), but now increase the coupling strength to the random input signals step-wise from η = 0.5 to η = 1 and finally to η = 2 (Figure 5). We observe that by this way also the higher density parts of the chaotic regime become eventually available for information import (green color).

Next, we increase the coupling strength η gradually from zero to a very large value of 20, at which the random input already dominates the system dynamics. For this numerical experiment, we keep the balance and density parameters fixed at b = −0.5, d = 0.5 (oscillatory regime), b = 0, d = 0.5 (chaotic regime), and b = 0.5, d = 0.5 (fixed point regime), respectively.

When in the fixed point regime (Figure 6f), we find that the dependence of the state-to-state correlation C(s_t, s_t+1) on the coupling strength η has the shape of a “resonance peak.” Since η effectively controls the amplitude of “noise” (used by us as pseudo input) added to the system, this corresponds to the phenomenon of “Recurrence Resonance” (RR), which we have previously found in three-neuron motifs (Krauss et al., 2019a): At small noise levels η, the system is stuck in the fixed point attractor, but adding an optimal amount of noise (so that C(s_t, s_t+1) becomes maximal) is freeing the system from this attractor and thus makes recurrent information “flux” possible, even in the fixed point regime. Adding too much noise is however counter-productive and leads to a decrease of C(s_t, s_t+1), as the system dynamics then becomes dominated by noise. We do not observe recurrence resonance in the other two dynamic regimes (Figures 6b,d).

Interestingly, we find very pronounced resonance-like curves also in the dependence of the input-to-state correlation C(x_t, s_t+1) on the coupling strength η, for all dynamical regimes (Figures 6a,c,e). Since C(x_t, s_t+1) is a measure of information import, we call this novel phenomenon “Import Resonance” (IR).

Next, we investigate the ability of the system to import more regular input signals with built-in temporal correlations, as well as inputs that are identical for all neurons. For this purpose, we feed all neurons with the same sinusoidal input signal, using an amplitude of a_sin = 1, an oscillation period of T_sin = 25 time steps, and a coupling strength of η = 2 (Figure 7). The density parameter is again fixed at d = 0.5, while the balance increases from b = −0.6 to b = +0.6 in five steps. We find that the input signal does not affect the evolution of neural states when the system is too far in the oscillatory phase or too far in the fixed point phase (c,g). Only systems where excitatory and inhibitory connections are approximately balanced are capable of information import (d-f). For b = −0.3 (d), most of the neurons are still part of the periodic attractor, but a small sub-population of neurons is taking up the external input signal (d). Interestingly, the system state is reflecting the periodic input signal even in the middle of the chaotic phase (e).

So far, we have analyzed input-to-state and state-to-state correlations mainly for a lag-time Δt = 1. We finally extend this analysis to larger lag-times up to 50 time steps (Figure 8), however only for three selected RNNs in the oscillatory, chaotic and fixpoint regime, using again our standard parameters (w = 0.5, N = 100, η = 0.5). Since our correlation measures C(x_t, s_t+1) (left column) and C(s_t, s_t+1) (right column) are defined as RMS averages, these values never fall below a certain noise level, which is in our case about 0.034. Another consequence of the RMS-average is that perfectly oscillatory RNN states with a period of two show up as C(s_t, s_t+1) = 1 (Figure 8b).

In the oscillatory regime, we find that input-to state correlations (as a measure of information import) remain at the noise level for all lag-times (a), while the system states are bound in a perfectly periodic attractor (b). Also in the fixpoint regime, both types of correlation are negligible for all non-zero lag-times. But remarkably, information can be imported (c) and stored (d) to a small but significant extent even in the middle of the chaotic regime, although the correlations decay back to noise level after about 20 time steps for this specific point in phase space (b = 0, d = 0.5). Future work will analyze how this correlation decay time depends on the statistical system parameters b and d.

We consider a system of N_neu neurons without biases, which are mutually connected according to a weight matrix {w_mk}, where w_mk denotes the connection strength from neuron k to neuron m. The weight matrix is random, but controlled by three statistical parameters, namely the “density” d of non-zero connections, the excitatory/inhibitory “balance” b, and the “width” w of the Gaussian distribution of weight magnitudes. The density ranges from d = 0 (isolated neurons) to d = 1 (fully connected network), and the balance from b = −1 (purely inhibitory connections) to b = +1 (purely excitatory connections). The value of b = 0 corresponds to a perfectly balanced system.

In order to construct a weight matrix with given parameters (b, d, w), we first generate a matrix Mmn(magn) of weight magnitudes, by drawing the Nneu2 matrix elements independently from a zero-mean normal distribution with standard deviation w and then taking the absolute value. Next, we generate a random binary matrix Bmn(nonz)∈{0,1}, where the probability of a matrix element being 1 is given by the density d, i.e., p₁ = d. Next, we generate another random binary matrix Bmn(sign)∈{-1,+1}, where the probability of a matrix element being +1 is given by p₊₁ = (1 + b)/2 where b is the balance. Finally, the weight matrix is constructed by elementwise multiplication wmn=Mmn(magn)· Bmn(nonz)· Bmn(sign). Note that throughout this paper, the width parameter is set to w = 0.5.

The momentary state of the RNN is given by the vector s(t) = {s_m(t)}, where the component s_m(t) ∈ [−1, +1] is the activation of neuron m at time t. The initial state s(t = 0) is set by assigning to the neurons statistically independent, normally distributed random numbers with zero mean and a standard deviation of σ_ini = 1.

We then compute the next state vector by simultaneously updating each neuron m according to (1)sm(t+1)=2πarctan(η xm(t)+∑k=1Nwmksk(t)). Here, x_m(t) are the external inputs of the RNN and η is a global “coupling strength”. Note that the input time series x_m(t) can, but must not be different for each neuron. In one type of experiment, we set the x_m(t) to independent, normally distributed random signals with zero mean and unit variance. In another experiment, we set all x_m(t) to the same oscillatory signal x(t) = a_sin · sin(2πt/T_sin).

After simulating the sequence of system states for N_stp = 1000 time steps, we analyze the properties of the state sequence (see below). For this evaluation, we disregard the first N_tra = 100 time steps, in which the system may still be in a transitory state that depends strongly on the initial condition. The simulations are repeated N_run = 10 times for each set of control parameters (b, d, η).

Consider a vector u(t) in M dimensions and a vector v(t) in N dimensions, both defined at discrete time steps t. The components of the vectors are denoted as u_m(t) and v_n(t). In order to characterize the correlations between the two time-dependent vectors by a single scalar quantity C(u_t, v_t+1), we proceed as follows:

First, we compute for each vector component m the temporal mean, (2)μum=〈um(t)〉t and the corresponding standard deviation (3)σum=〈(um(t)-μum)2〉t. Based on this, we compute the M × N pairwise (Pearson) correlation matrix, (4)Cmn(uv)=〈[um(t)-μum]·[vn(t+1)-μvn]〉tσumσvn, defining Cmn(uv)=0 whenever σ_um = 0 or σ_vn = 0.

Finally we compute the root-mean-squared average of this matrix, (5)C(ut,vt+1)=RMS{Cmn(uv)}mn=1MN∑m=1M∑n=1N|Cmn(uv)|2 This measure is applied in the present paper to quantify the correlations C(s_t, s_t+1) between subsequent RNN states, as well as the correlations C(x_t, s_t+1) between the momentary input and the subsequent RNN state.

In addition to the linear correlations, we consider the mutual information between the two vectors u(t) and v(t), in order to capture also possible non-linear dependencies. However, since the full computation of this quantity is computationally extremely demanding, we binarize the continuous vector components and then consider only the pairwise mutual information between these binarized components.

For the binarization, we first subtract the mean values from each of the components, (6)um(t)→Δum(t)=um(t)-μum. We then map the continuous signals Δu_m(t) ∈ [−∞, +∞] onto two-valued bits u^_m(t) ∈ {0, 1} by defining u^_m(t) = 0 if Δu_m(t) < 0 and u^_m(t) = 1 if Δu_m(t) > 0. In the case of a tie, Δu_m(t) = 0, we set u^_m(t) = 0 with a probability of 1/2.

We next compute the pairwise joint probabilities P(u^m,v^n) by counting how often each of the four possible bit combinations occurs during all available time steps. From that we also obtain the marginal probabilities P(u^_m) and P(v^n).

The matrix of pairwise mutual information is then defined as (7)Imn(uv)=∑u^m=0,1∑v^n=0,1P(u^m,v^n)log[P(u^m,v^n)P(u^m)·P(v^n)], defining all terms as zero where P(u^_m) = 0 or P(v^n) = 0.

Finally we compute the mean over all matrix elements (each ranging between 0 and 1 bit), (8)I(ut,vt+1)=MEAN{Imn(uv)}mn=1MN∑m=1M∑n=1NImn(uv) This measure is applied in the present paper to quantify the mutual information I(s_t, s_t+1) between subsequent RNN states, as well as the mutual information I(x_t, s_t+1) between the momentary input and the subsequent RNN state.