Suppression of chaos in a partially driven recurrent neural network
The dynamics of recurrent neural networks (RNNs), and particularly their
response to inputs, play a critical role in information processing. In many
applications of RNNs, only a specific subset of the neurons generally receive
inputs. However, it remains to be theoretically clarified how the restriction
of the input to a specific subset of neurons affects the network dynamics.
Considering RNNs with such restricted input, we investigate how the proportion,
$p$, of the neurons receiving inputs (the "inputs neurons") and the strength of
the input signals affect the dynamics by analytically deriving the conditional
maximum Lyapunov exponent. Our results show that for sufficiently large $p$,
the maximum Lyapunov exponent decreases monotonically as a function of the
input strength, indicating the suppression of chaos, but if $p$ is smaller than
a critical threshold, $p_c$, even significantly amplified inputs cannot
suppress spontaneous chaotic dynamics. Furthermore, although the value of $p_c$
is seemingly dependent on several model parameters, such as the sparseness and
strength of recurrent connections, it is proved to be intrinsically determined
solely by the strength of chaos in spontaneous activity of the RNN. This is to
say, despite changes in these model parameters, it is possible to represent the
value of $p_c$ as a common invariant function by appropriately scaling these
parameters to yield the same strength of spontaneous chaos. Our study suggests
that if $p$ is above $p_c$, we can bring the neural network to the edge of
chaos, thereby maximizing its information processing capacity, by amplifying
inputs.
Shotaro Takasu, Toshio Aoyagi
# Suppression of chaos in a partially driven recurrent neural network

###### Abstract

The dynamics of recurrent neural networks (RNNs), and particularly their response to inputs, play a critical role in information processing. In many applications of RNNs, only a specific subset of the neurons generally receive inputs. However, it remains to be theoretically clarified how the restriction of the input to a specific subset of neurons affects the network dynamics. Considering recurrent neural networks with such restricted input, we investigate how the proportion, \(p\), of the neurons receiving inputs (the "inputs neurons") and a quantity, \(\xi\), representing the strength of the input signals affect the dynamics by analytically deriving the conditional maximum Lyapunov exponent. Our results show that for sufficiently large \(p\), the maximum Lyapunov exponent decreases monotonically as a function of \(\xi\), indicating the suppression of chaos, but if \(p\) is smaller than a critical threshold, \(p_{c}\), even significantly amplified inputs cannot suppress spontaneous chaotic dynamics. Furthermore, although the value of \(p_{c}\) is seemingly dependent on several model parameters, such as the sparseness and strength of recurrent connections, it is proved to be intrinsically determined solely by the strength of chaos in spontaneous activity of the RNN. This is to say, despite changes in these model parameters, it is possible to represent the value of \(p_{c}\) as a common invariant function by appropriately scaling these parameters to yield the same strength of spontaneous chaos. Our study suggests that if \(p\) is above \(p_{c}\), we can bring the neural network to the edge of chaos, thereby maximizing its information processing capacity, by adjusting \(\xi\).

Large-scale recurrent neural networks (RNNs) exhibit various dynamical patterns, including limit cycles and chaos. They have been used to model brain functions such as working memory [1], motor control [2], and context-based learning [3]. The rich dynamics exhibited by an RNN can also be used for information processing. Reservoir computing (RC) [4; 5] is a machine learning framework that utilizes large RNNs, called "reservoirs", to reproduce the time series data of interest. RC enables quicker learning in RNNs than conventional backpropagation-based methods, because with this method, only the output weights are updated, with the other recurrent weights left unchanged. RC is not restricted to RNNs, and indeed a wide range of dynamical systems can serve as the reservoir under appropriate conditions. Physical reservoir computing, in which a real physical system is used as the reservoir, has been an area of active research in recent years [6].

In general, RNNs must possess rich dynamics in order to assimilate a diverse set of signals to be learned. For this reason, it is advantageous for RNNs to exhibit chaotic spontaneous activity. On the other hand, for an RNN to successfully reproduce a target time series, it must converge to the same state each time it receives a particular set of input signals, regardless of its initial internal state. This property is known as the "echo state property (ESP)" [4] in the context of RC. Hence, it is hypothesized that an RNN that displays varied spontaneous activity while maintaining consistency in response to inputs will exhibit superior computational performance. Such an RNN is commonly referred to as being at the "edge of chaos", and it is known empirically that reservoirs in this regime have the highest computational capacity [7]. In fact, there is experimental evidence suggesting that mammalian neuronal networks operate in this critical regime [8; 9].

Lyapunov spectrum analysis [10] allows us to study the dynamics of RNNs in a quantitative manner. The maximum Lyapunov exponent (MLE) characterizes the exponential rate of separation of infinitesimally close trajectories. A dynamical system with a positive MLE exhibits chaotic behavior. In the case that a dynamical system is driven by given input signals, the MLE is called the maximum conditional Lyapunov exponent (MCLE). Two identical RNNs with slightly different initial states will converge to the same state under the same inputs if and only if the MCLE is negative. Therefore, a negative MCLE is necessary for an RNN to exhibit consistency with respect to inputs.

The MLE and MCLE of an RNN with random weights can be analytically computed in the limit of a large network size, using the dynamic mean-field approach [11]. It is known that the MCLE of a random RNN decreases with the strength of the input signal and eventually becomes negative in both the cases that the signal is white noise [12; 13; 14; 15] and deterministic [16]. In other words, sufficiently amplified driving input signals can suppress chaotic dynamics. These findings suggest that it is possible to shift the state of an RNN exhibiting chaotic spontaneous dynamics toward the edge of chaos by appropriately amplifying the input, and then use its shifted state as an efficient reservoir. It is known that FORCE learning can efficiently train a chaotic RNN with feedback connections, with the feedback input suppressing chaos [17]. There is also experimental evidence indicating similar suppression of chaos in the brain. Specifically, extensive neural recordings taken with monkeys and cats show that visual stimuli can attenuate spontaneous fluc
tuations in cortical neurons [18].

Previous studies of the Lyapunov exponents of large random RNNs assume a model in which every unit in the RNN connects to the input layer and receives driving signals. Hereafter, we refer to RNNs of this type as "full-input RNNs". However, such a model is biologically implausible because biological synapses are sparse [19]. Additionally, in physical reservoir computing, it is often unfeasible to connect the input to all reservoir units. Therefore, it is important to investigate the dynamics of RNNs in the case that only a subset of the neurons receive input signals. We refer to RNNs of this type as "partial-input RNNs." It remains to be elucidated whether sufficiently amplified inputs always suppress the chaotic activity in a partial-input RNN, as in the case of a full-input RNN. In this work, we address this question by analitically calculating the MCLE of a partial-input random RNN.

We investigate the discrete-time dynamics of a random sparse RNN with \(N\gg 1\) neurons of which only \(pN\) receive inputs (Fig.1). We define the parameter \(p\in[0,1]\), called the "input partiality," which determines the fraction of neurons coupling to the input unit. We consider the case in which the input signal \(s(t)\) at time \(t\) is a scalar for simplicity, although our theory can be straightforwardly extended to multi-dimensional inputs. The state of a neuron that receives inputs is represented by a dynamical variable \(x_{i}(t)\in\mathbb{R}\) (\(i=1,\cdots,pN\)) whose evolution obeys the equation

\[x_{i}(t+1)=\sum_{j=1}^{pN}J_{ij}\phi(x_{j}(t))+\sum_{k=pN+1}^{N} J_{ik}\phi(y_{k}(t))+u_{i}s(t), \tag{1}\]

while the state of a neuron that does not receive inputs is represented by a dynamical variable \(y_{i}(t)\in\mathbb{R}\) (\(i=pN+1,\cdots,N\)) whose evolution obeys the equation

\[y_{i}(t+1) =\sum_{j=1}^{pN}J_{ij}\phi(x_{j,t})+\sum_{k=pN+1}^{(1-p)N}J_{ik} \phi(y_{k,t}). \tag{2}\]

Here, \(u_{i}\) (\(i=1,\cdots,pN\)) is the coupling weight connecting the input signal to the \(i\)th neuron, \(J_{ij}\) is a recurrent weight matrix determining the coupling from the \(j\)th to the \(i\)th neuron, and \(\phi\) is the activation function. The value of \(u_{i}\) is drawn randomly from a Gaussian distribution with zero mean and unit variance. We define \(J_{ij}\) to be a random sparse matrix whose elements are chosen as follows. The value of each element is assigned independently. For a given element, first its value is randomly chosen to be zero or non-zero with respective probabilities \(1-\alpha\) and \(\alpha\) (where \(\alpha\in(0,1]\)), and then, in the latter case, it is assigned a value drawn from a Gaussian distribution with zero mean and variance \(g^{2}/N\). The gain parameter \(g\) represents the recurrent coupling strength. The activation function is chosen as \(\phi(x)=\mathrm{erf}(\frac{\sqrt{\pi}}{2}x)\) for analytic tractability. It has been found that in the absence of inputs, this RNN exhibits a transition from fixed-point dynamics to chaotic dynamics at \(g=1\) in the limit of a large network [11]. In previous studies, zero mean white noise has typically been used for the input signal \(s(t)\)[12; 13; 14; 15], but we do not limit the choice of \(s(t)\) in this way, and instead regard it to be an arbitrary time series, except when stated otherwise.

In the investigation of the RNN defined by the above equations, the model parameters to be varied are the input partiality, \(p\), the recurrent connection sparsity, \(\alpha\), and the recurrent coupling strength, \(g\).

We can obtain the statistical properties of \(x_{i}(t)\) and \(y_{i}(t)\) using a mean-field approach [13; 14] in the limit \(N\to\infty\). As seen from Eqs.(1) and (2), \(x_{i}(t+1)\) and \(y_{i}(t+1)\) are obtained as sums of large numbers of identically distributed independent variables, \(\{J_{ij}\phi(x_{j}(t))\}_{j=1}^{pN}\) and \(\{J_{ik}\phi(y_{k}(t))\}_{k=pN+1}^{N}\). Thus, according to the central limit theorem, we can consider \(x_{i}(t)\) and \(y_{i}(t)\) to follow Gaussian distributions. It thus suffices to determine their averages, \(\langle x_{i}(t)\rangle\) and \(\langle y_{i}(t)\rangle\), and variances, \(\langle x_{i}^{2}(t)\rangle\) and \(\langle y_{i}^{2}(t)\rangle\), where \(\langle\cdots\rangle\) denotes the average over realizations of the quenched weights \(J_{ij}\) and \(u_{i}\). Introducing the notation \(K(t):=\langle y_{i}(t)^{2}\rangle\), the Gaussian distributions exhibited by \(x_{i}(t)\) and \(y_{i}(t)\) can be expressed as follows (see Supplemental Material for derivation [20]:

\[x_{i}(t+1) \sim\mathcal{N}(0,K(t+1)+s(t)^{2}), \tag{3a}\] \[y_{i}(t+1) \sim\mathcal{N}(0,K(t+1)). \tag{3b}\]

The value of \(K(t)\) obeys the recurrence relation,

\[K(t+1) =-\alpha g^{2}+\frac{4}{\pi}\alpha g^{2}\bigg{(}(1-p)\arctan \sqrt{1+\pi K(t)} \tag{4}\] \[+p\arctan\sqrt{1+\pi(K(t)+s(t-1)^{2})}\bigg{)}.\]

Sequentially substituting the input time series \(\{s(t)\}_{t}\) into Eq.(4), we can obtain \(\{K(t)\}_{t}\).

The time series \(\{K(t)\}_{t}\) allows us to derive the MCLE, \(\lambda\), of the RNN [12]. The MCLE is defined as the asymptotic growth rate of the distance between two replicated

Figure 1: A schematic depiction of the partial-input RNN studied in this work. The shaded region represents the neurons that receive input signals through input connectivity.
RNNs,

\[\lambda :=\lim_{T\rightarrow\infty,\mathbf{\delta}(0)\to 0}\frac{1}{T}\log \frac{\|\mathbf{\delta}(T)\|}{\|\mathbf{\delta}(0)\|}, \tag{5}\]

where \(\|\mathbf{\delta}(t)\|\) denotes the distance at time \(t\) between two replicas receiving the same input signals. We can calculate \(\lambda\) with random matrix theory [21] (see Supplemental Material for derivation [20]), obtaining

\[\lambda =\lim_{T\rightarrow\infty}\frac{1}{T}\sum_{t=0}^{T}\frac{1}{2} \log\alpha g^{2}\bigg{(}\frac{1-p}{\sqrt{1+\pi K(t)}} \tag{6}\] \[+\frac{p}{\sqrt{1+\pi(K(t)+s(t-1)^{2})}}\bigg{)}.\]

Because \(\{K(t)\}_{t}\) has been determined by Eq.(4), we can obtain \(\lambda\) by substituting \(\{K(t)\}_{t}\) and \(\{s(t)\}_{t}\) into Eq.(6).

Figure 2(a) displays the MCLE, \(\lambda\), as a function of the input intensity \(\sigma\) defined below, for \(p=0.4\) and \(p=0.6\), with \(g=3.0\) and \(\alpha=1.0\). Here we assume the input to be Gaussian white noise with zero mean and standard deviation \(\sigma\). The analytic results obtained from Eq.(6) are consistent with the results of the numerical simulations. For both values of \(p\), \(\lambda\) decreases monotonically as a function of \(\sigma\). However, while the MCLE for \(p=0.6\) falls below \(0\) around \(\sigma=20\), the MCLE for \(p=0.4\) remains positive throughout the range of \(\sigma\) shown in the figure.

In this work, our main interest is to determine whether the MCLE falls below \(0\) with sufficiently amplified inputs. However, it cannot be answered by observing the numerical simulations shown in Figure.2(a) even if simulations are performed for quite large \(\sigma\), which motivates us to derive the analytic value of the MCLE conditioned by infinitely large inputs, \(\lambda_{\infty}\). To answer the above question, it is sufficient to determine the sign of \(\lambda_{\infty}\). We introduce a scaling parameter of input signals, \(\xi\in\mathbb{R}\), and then represent the amplified input signals by \(\{\xi s(t)\}_{t}\). The MCLE conditioned by infinitely large inputs is denoted by \(\lambda_{\infty}:=\lim_{|\xi|\rightarrow\infty}\lambda\).

The quantity \(\lambda_{\infty}\) can be derived analytically as follows. As the value of \(\xi\) becomes sufficiently large, \(K(t)\) approaches a limiting constant value \(K_{\infty}\) given by

\[K_{\infty}=-\alpha g^{2}+\frac{4}{\pi}\alpha g^{2}\bigg{(}\frac {\pi}{2}p+(1-p)\arctan\sqrt{1+\pi K_{\infty}}\bigg{)}. \tag{7}\]

This expression for \(K_{\infty}\) is obtained by replacing \(s(t)\) with \(\xi s(t)\) in Eq.(4) and considering the limit \(|\xi|\rightarrow\infty\). Taking this limit and substituting \(K_{\infty}\) for \(K(t)\) in Eq.(6), we obtain

\[\lambda_{\infty}=\frac{1}{2}\log\alpha g^{2}\bigg{(}\frac{1-p}{ \sqrt{1+\pi K_{\infty}}}\bigg{)}. \tag{8}\]

We have confirmed that the value of \(\lambda_{\infty}\) obtained from Eqs.(7) and (8) is consistent with the results of numerical simulations (Fig.2(b)).

Figure 2(b) depicts the relationship between \(p\) and \(\lambda_{\infty}\) for various values of the recurrent weight intensity \(g\). We define the value of \(p\) at \(\lambda_{\infty}=0\) as the "critical input partiality", \(p_{c}\). Because the condition \(\lambda_{\infty}>0\) always holds for \(p<p_{c}\), we conclude that even sufficiently amplified input signals cannot suppress the chaotic activity of the RNN if \(p<p_{c}\).

Figure 3(a) depicts the dependence of \(p_{c}\) on \(g\) with fixed \(\alpha\). We obtain the curve by solving Eq.(8) for \(K_{\infty}\) with \(\lambda_{\infty}=0\) and substituting this \(K_{\infty}\) into Eq.(7), yield

Figure 2: (a) The maximum conditional Lyapunov exponent (MCLE), \(\lambda\), calculated for various values of the input partiality, \(p\). The theoretical form given in Eq.(6) is plotted with a solid curve for \(p=0.4\) and a dashed curve for \(p=0.6\), where the time average in Eq.(6) is computed up to \(T=10^{5}\). Error bars represent \(\pm\)std of direct numerical simulations based on the definition of \(\lambda\) in Eq.(5). For each plot, the coupling strength is set to \(g=3.0\), and the sparsity is set to \(\alpha=1.0\). The input signal \(s(t)\) is Gaussian white noise with mean \(0\) and variance \(\sigma^{2}\). For the numerical simulations, the value of the MCLE was obtained by directly calculating Eq.(5) for a sufficiently long time, i.e., \(T=10^{4}\). (b) Analytic results (solid, dashed, dot-dashed curves) and numerical results (error bars indicating \(\pm\)std) for \(\lambda_{\infty}\) with sparsity \(\alpha=1.0\). In numerical simulations, the value of \(\lambda_{\infty}\) was obtained by directly calculating Eq.(5) for a sufficiently large input magnitude (\(\sigma=10^{3}\)). In the cases of both (a) and (b), the values of the MCLE (\(\lambda\) and \(\lambda_{\infty}\)) were calculated for \(10\) different network realizations with a network size of \(N=1000\).
ing

\[p_{c} =1-\frac{1}{\alpha g^{2}}\biggl{[}1-\pi\alpha g^{2}\] \[+4\alpha g^{2}\biggl{(}\frac{\pi}{2}p_{c}+(1-p_{c})\arctan\left((1- p_{c})\alpha g^{2}\right)\biggr{)}\biggr{]}^{1/2}.\]

Below the curve, \(\lambda_{\infty}\) is positive, and thus, in this region, chaos is not suppressed no matter how strong the input because \(\lambda\) is a decreasing function of the input intensity as shown in Figure 2(a). As seen in Figure 3(a), \(p_{c}\) is an increasing function of \(g\). This implies that as the values of the coupling strength increases, larger values of \(p\) are needed to suppress chaos.

We next investigate the effect of sparsity \(\alpha\) on \(p_{c}\). Plotting the \(g\)-\(p_{c}\) curves with several values of \(\alpha\), we find that a sparser RNN results in a smaller value of \(p_{c}\) (Fig.3(a)). This is intuitively understandable, because the dynamics of a sparser RNN are less chaotic, and thus a smaller value of the input partiality is sufficient to control the chaos. To take account of this relationship, we introduce the MLE of the RNN with no input, denoted by \(\lambda_{0}\). Clearly, \(\lambda_{0}\) quantifies the strength of chaos in spontaneous activity of the RNN, and it can be determined analytically by substituting \(s(t)\equiv 0\) into Eq.(4) and Eq.(6), yielding

\[K =\alpha g^{2}\left(-1+\frac{4}{\pi}\arctan\sqrt{1+\pi K}\right), \tag{10}\] \[\lambda_{0} =\frac{1}{2}\log\frac{\alpha g^{2}}{\sqrt{1+\pi K}}. \tag{11}\]

Interestingly, we find that when \(p_{c}\) is plotted with respect to \(\lambda_{0}\), the resulting curves for all values of \(\alpha\) coincide, as seen in Figure 3(b). The reason for this coincidence is easily understood by considering Eqs.(9)-(11). From Eqs.(10) and (11), we see that \(\lambda_{0}\) is a function of \(\alpha g^{2}\). Writing the corresponding inverse function as \(\alpha g^{2}=f(\lambda_{0})\), and substituting this into Eq.(9), we obtain \(p_{c}\) expressed as a function of \(\lambda_{0}\) alone. This finding implies that \(p_{c}\) depends primarily on the strength of spontaneous chaos, independently of how sparse the recurrent connection is.

Finally, we study the information processing capability of a partial-input RNN employed as a reservoir for RC. Memory capacity [22] is a commonly used benchmarks for RC. It is a measure of the ability of a reservoir to perform short-term memory tasks through the reconstruction of its past input signals. We define the memory capacity as follows. From the \(N\) reservoir units, \(K\) (\(1\leq K\leq N\)) lead-out units are randomly chosen, and represented by a vector \(\mathbf{\tilde{x}}(t)\in\mathbb{R}^{K}\). The reservoir's output is defined as \(\hat{z}(t):=\mathbf{w}^{\top}\mathbf{\tilde{x}}(t)\), where the vector \(\mathbf{w}\in\mathbb{R}^{K}\) represents the output weights. In a \(\tau\)-delay memory task, the reservoir at time \(t\) is required to output the previous input signal \(s(t-\tau)\), and the output weights are trained to minimize the mean squared error between \(\hat{z}(t)\) and the desired output, \(s(t-\tau)\). This is accomplished with a least-squares method, and the trained output weights are determined as \(\mathbf{\hat{w}}=(\mathbf{XX}^{\top})^{-1}\mathbf{X}\mathbf{s}\), where \(\mathbf{X}:=(\mathbf{\tilde{x}}(1)\cdots\mathbf{\tilde{x}}(T))\) and \(\mathbf{s}:=(s(1)\cdots s(T))^{\top}\) (\(T\) being the length of the simulation). After training, we evaluate the task performance \(M_{\tau}\) defined as

\[M_{\tau}:=1-\frac{\langle(\hat{z}(t)-s(t-\tau))^{2}\rangle}{\langle s(t)^{2} \rangle}, \tag{12}\]

where the brackets represent the time average. Because the numerator of the second term in Eq.(12) is the mean squared error, \(M_{\tau}\) approaches \(1\) as the reservoir learns to accurately reconstruct its past input \(s(t-\tau)\). The memory capacity \(MC\) is defined as the sum of the \(M_{\tau}\), \(MC:=\sum_{\tau=1}^{\infty}M_{\tau}\). It has been mathematically proved that \(MC\) satisfies the inequality \(0\leq MC\leq K\)[22; 23].

Assuming that the input signal \(s(t)\) is Gaussian white noise with zero mean and variance \(\sigma^{2}\), we calculated both the MCLE, \(\lambda\), and memory capacity, \(MC\). The results are plotted as functions of \(\sigma\) in Figure 4(a) and (b). As previously noted, an increase in the input magnitude leads to a decrease in the MCLE (as seen in Fig.4(a)),

Figure 3: (a) Relation between the coupling strength, \(g\), and the critical input partiality, \(p_{c}\), given by Eq.(9) for several values of the sparsity, \(\alpha\). (b) The same data as in (a) plotted with respect to \(\lambda_{0}\) rather than \(g\). It is seen that when plotted in this manner, the data for \(p_{c}\) fall along the same curve for each value of \(\alpha\) considered.
with the result that the plot in Fig.4(c) shifts leftward as \(\sigma\) increases. From Eq.(9), \(p_{c}\) is found to be approximately \(0.074\) under the conditions employed in Fig.4. When \(p=0.15\) and \(p=0.50\), the plots intersect the vertical line \(\lambda=0\) in Fig.4(c), as our theory predicts, and the memory capacity reaches its maximum value near \(\lambda=0\). Contrastingly, the plots with \(p=0.05\) and \(p=0.07\) remain in the chaotic domain (\(\lambda>0\)), and the memory capacity remains relatively low. It is thus seen that once input connections have been built such that \(p\) exceeds \(p_{c}\), optimal computational capability can be realized only by amplifying the input signals appropriately. This finding should be helpful for the physical reservoir computing paradigm, because amplifying input signals is generally easier and more cost effective than adding new input connections.

In the present work, we have examined a partial-input RNN with rate neurons and have analytically shown the existence of a critical input partiality \(p_{c}\) that determines whether the chaotic activity can be suppressed by input signals. Our theory can be applied to realistic situations in which an RNN receives input signals with non-Gaussian statistics or temporal correlations, because we do not assume any particular statistics of the input signals, in contrast to previous theoretical works [12; 13; 14; 15]. In a future work, we will investigate whether there exists critical input partiality in other types of RNNs, such as those with heavy-tailed recurrent weights [24] or unsaturated activation functions [25]. In addition, we have confirmed that memory capacity is maximized near the critical value of MCLE, \(\lambda=0\), which corresponds to the "edge of chaos." Our theory suggests that we can readily construct a partial-input RNN at the edge of chaos by tuning the input signals, as long as the input partiality exceeds \(p_{c}\). The present study provides a possible novel approach to designing reservoir computing.

S.T. was supported by JSPS KAKENHI Grant No. JP22J21559. T.A. was supported by JSPS KAKENHI Grant No.JP20K21810, No.JP20H04144, and No.JP20K20520.

## References

* Rajan _et al._ [2016]K. Rajan, C. Harvey, and D. Tank, Neuron **90**, 128 (2016).
* Laje and Buonomano [2013]R. Laje and D. Buonomano, Nat. neurosci. **16**, 925 (2013).
* Enel _et al._ [2016]P. Enel, E. Procyk, R. Quilodran, and P. F. Dominey, PLOS Comput. Biol. **12**, 1 (2016).
* Jaeger [2001]H. Jaeger, _The "echo state" approach to analysing and training recurrent neural networks_, Tech. Rep. GMD Report 148 (German National Research Center for Information Technology, 2001).
* Maass _et al._ [2002]W. Maass, T. Natschlager, and H. Markram, Neural Comput. **14**, 2531 (2002).
* Nakajima and Fischer [2021]K. Nakajima and I. Fischer, _Reservoir Computing: Theory, Physical Implementations, and Applications_ (Springer, 2021).
* Bertschinger and Natschlager [2004]N. Bertschinger and T. Natschlager, Neural Comput. **16**, 1413 (2004).
* Morales _et al._ [2023]G. B. Morales, S. di Santo, and M. A. Munoz, Proceedings of the National Academy of Sciences **120**, e2208998120 (2023).
* Dahmen _et al._ [2019]D. Dahmen, S. Grun, M. Diesmann, and M. Helias, Proceedings of the National Academy of Sciences **116**, 13051 (2019).
* Pikovsky and Politi [2016]A. Pikovsky and A. Politi, _Lyapunov Exponents: A Tool to Explore Complex Dynamics_ (Cambridge University Press, 2016).
* Sompolinsky _et al._ [1988]H. Sompolinsky, A. Crisanti, and H. J. Sommers, Phys. Rev. Lett. **61**, 259 (1988).
* Molegedy _et al._ [1992]L. Molegedy, J. Schuchhardt, and H. G. Schuster, Phys. Rev. Letters **69**, 3717 (1992).
* Massar and Serge [2013]M. Massar and M. Serge, Phys. Rev. E **87**, 042809 (2013).
* Haruna and Nakajima [2019]T. Haruna and K. Nakajima, Phys. Rev. E **100**, 062312 (2019).
* Schuecker _et al._ [2018]J. Schuecker, S. Goedeke, and M. Helias, Phys. Rev. X **8**, 041029 (2018).
* Rajan _et al._ [2010]K. Rajan, L. F. Abbott, and H. Sompolinsky, Phys. Rev. E **82**, 011903 (2010).
* Sussilo and Abbott [2009]D. Sussilo and L. F. Abbott, Neuron **63**, 544 (2009).

Figure 4: Relation between the memory capacity, \(MC\), and the maximum conditional Lyapunov exponent (MCLE), \(\lambda\), for network size \(N=1000\), coupling strength \(g=1.5\), sparsity \(\alpha=1.0\), and number of lead-out nodes \(K=10\). The values of (a) \(\lambda\) and (b) \(MC\) are respectively plotted as functions of the standard deviation of the input signals, \(\sigma\). (c) Each plot represents \((\lambda,MC/K)\) with various values of \(\sigma\) (\(0.01\leq\sigma\leq 20\)). The sum of the \(M_{\tau}\) is calculated up to \(\tau=500\).
* Churchland _et al._ [2010]M. M. Churchland _et al._, Nat. Neurosci. **13**, 369 (2010).
* Wildenberg _et al._ [2021]G. A. Wildenberg, M. R. Rosen, J. Lundell, D. Paukner, D. J. Freedman, and N. Kasthuri, Cell Rep. **36**, 109709 (2021).
* [20]Supplemental Material.
* Ahmadian _et al._ [2015]Y. Ahmadian, F. Fumarola, and K. D. Miller, Phys. Rev. E **91**, 012820 (2015).
* Jaeger [2001]H. Jaeger, _Short term memory in echo state networks_, Tech. Rep. GMD Report 152 (German National Research Center for Information Technology, 2001).
* Dambre _et al._ [2012]J. Dambre, D. Verstraeten, B. Schrauwen, and S. Massar, Sci. Rep. **2**, 514 (2012).
* Kusmierz _et al._ [2020]L. Kusmierz, S. Ogawa, and T. Toyoizumi, Phys. Rev. Lett. **125**, 028101 (2020).
* Ahmadian _et al._ [2013]Y. Ahmadian, D. B. Rubin, and K. D. Miller, Neural Computation **25**, 1994 (2013).
**Supplemental Material: Suppression of chaos in a partially driven recurrent neural network**

Shotaro Takasu1 and Toshio Aoyagi

Footnote 1: shotaro.takasu.63x@st.kyoto-u.ac.jp

_Graduate School of Informatics, Kyoto University, Kyoto 606-8501, Japan_

Footnote 2: footnotemark:

Footnote 3: footnotemark:

Footnote 4: footnotemark:

November 3, 2021

## I Mean-field theory for a partially driven recurrent neural network

Here we derive statistical properties of a partially driven recurrent neural network (Eqs.(3)-(4) in the main text) applying a mean-field approach [1; 2]. As mentioned in the main text, the values of \(x_{i}(t)\) and \(y_{k}(t)\) obey Gaussian distributions in the limit of a large network size \(N\rightarrow\infty\) because of the central limit theorem. Thus, it is sufficient to determine their mean and variance. Taking the average of the evolution equation of \(x_{i}(t)\) and \(x_{i}(t)^{2}\) over realizations of weights \(J_{ij}\) and \(u_{i}\), we obtain

\[\langle x_{i}(t+1)\rangle =\sum_{j=1}^{pN}\langle J_{ij}\rangle\langle\phi(x_{j}(t))\rangle +\sum_{k=pN+1}^{N}\langle J_{ik}\rangle\langle\phi(y_{k}(t))\rangle+\langle u _{i}\rangle s(t)=0,\] (S.1) \[\langle x_{i}(t+1)^{2}\rangle =\sum_{j,j^{\prime}=1}^{pN}\langle J_{ij}J_{ij^{\prime}}\rangle \langle\phi(x_{j}(t))\phi(x_{j}^{\prime}(t))\rangle+\sum_{k,k^{\prime}=pN+1}^ {N}\langle J_{ik}J_{ik^{\prime}}\rangle\langle\phi(y_{k}(t))\phi(y_{k}^{ \prime}(t))\rangle+\langle u_{i}^{2}\rangle s(t)^{2}\] \[=p\alpha g^{2}\big{\langle}\phi(x_{i}(t))^{2}\big{\rangle}+(1-p) \alpha g^{2}\big{\langle}\phi(y_{k}(t))^{2}\big{\rangle}+s(t)^{2}.\] (S.2)

Similarly, \(\langle y_{i}(t+1)\rangle\) and \(\langle y_{i}(t+1)^{2}\rangle\) are given by

\[\langle y_{i}(t+1)\rangle =0,\] (S.3) \[\langle y_{i}(t+1)^{2}\rangle =p\alpha g^{2}\big{\langle}\phi(x_{i}(t))^{2}\big{\rangle}+(1-p) \alpha g^{2}\big{\langle}\phi(y_{k}(t))^{2}\big{\rangle}.\] (S.4)

Note that we have used the assumption that \(\phi(x_{j})\) and \(\phi(y_{k})\) is independent of its incoming weight \(J_{ij}\) and \(J_{ik}\). This assumption is justified in the limit \(N\rightarrow\infty\) as can be shown using the generating-function formalism [3; 4].

Introducing the notation \(K(t):=\langle y_{i}(t)^{2}\rangle\), the Gaussian distributions exhibited by \(x_{i}(t)\) and \(y_{i}(t)\) can be expressed as follows:

\[x_{i}(t+1) \sim\mathcal{N}(0,K(t+1)+s(t)^{2}),\] (S.5 a) \[y_{i}(t+1) \sim\mathcal{N}(0,K(t+1)).\] (S.5 b)

From Eqs.(S.4 ) and (S.5 ), we can directly calculate \(K(t+1)\), obtaining

\[K(t+1) =p\alpha g^{2}\int Dx\phi^{2}\left(\sqrt{K(t)+s(t-1)^{2}}x\right) +(1-p)\alpha g^{2}\int Dy\phi^{2}\left(\sqrt{K(t)}y\right)\] \[=-\alpha g^{2}+\frac{4}{\pi}\alpha g^{2}\left((1-p)\arctan\sqrt{ 1+\pi K(t)}+p\arctan\sqrt{1+\pi(K(t)+s(t-1)^{2})}\right),\] (S.6)

where \(\int Dx:=\int dx\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}\).

## II Derivation of the maximum conditioned Lyapunov exponent

Here we derive the maximum conditioned Lyapunov exponent \(\lambda\) analytically. Linearizing the evolution equation of the RNN, we obtain the variational equation describing the evolution of infinitesimal perturbation \(\mathbf{\delta}(t)\),

\[\mathbf{\delta}(t+1)=\mathbf{J}\mathbf{\phi^{\prime}}(t)\mathbf{\delta}(t),\] (S.7)
where the matrix \(\mathbf{\phi}^{\prime}(t)\) is the diagonal matrix whose \(i\)th diagonal entry is \(\phi^{\prime}(x_{i}(t))\) for \(1\leq i\leq pN\) or \(\phi^{\prime}(y_{i}(t))\) for \(pN+1\leq i\leq N\). The typical growth rate, \(\|\mathbf{\delta}(t+1)\|/\|\mathbf{\delta}(t)\|\), is determined by the spectral radius of the Jacobian \(\mathbf{J}\mathbf{\phi}^{\prime}(t)\). According to random matrix theory [5], the spectral radius, \(\rho(t)\), is given by

\[\rho(t)^{2}=p\alpha g^{2}\left\langle\phi^{\prime}(x(t))^{2}\right\rangle+(1- p)\alpha g^{2}\left\langle\phi^{\prime}(y(t))^{2}\right\rangle,\] (S.8)

in the limit of large network size. Applying the results of the mean-field theory (Eq.(S.5 )), we can calculate \(\lambda\), yielding

\[\lambda :=\lim_{T\to\infty,\mathbf{\delta}(0)\to 0}\frac{1}{T}\log\frac{\|\mathbf{ \delta}(T)\|}{\|\mathbf{\delta}(0)\|}=\lim_{T\to\infty,\mathbf{\delta}(0)\to 0}\frac{1}{ 2T}\sum_{t=0}^{T-1}\log\frac{\|\mathbf{\delta}(t+1)\|^{2}}{\|\mathbf{\delta}(t)\|^{2}}\] \[\approx\lim_{T\to\infty}\frac{1}{2T}\sum_{t=0}^{T-1}\log\rho(t)^ {2}\] \[=\lim_{T\to\infty}\frac{1}{2T}\sum_{t=0}^{T-1}\log\left(p\alpha g ^{2}\int Dx\phi^{\prime 2}\left(\sqrt{K(t)+s(t-1)^{2}}x\right)+(1-p)\alpha g^{2} \int Dy\phi^{\prime 2}\left(\sqrt{K(t)}y\right)\right)\] \[=\lim_{T\to\infty}\frac{1}{T}\sum_{t=0}^{T}\frac{1}{2}\log\alpha g ^{2}\left(\frac{p}{\sqrt{1+\pi(K(t)+s(t-1)^{2})}}+\frac{1-p}{\sqrt{1+\pi K(t) }}\right),\] (S.9)

where we used \(\phi(x)=\text{erf}(\frac{\sqrt{\pi}}{2}x)\) in the last equation.
