# Bayesian inference is facilitated by modular neural networks with different
time scales

Kohei Ichikawa
Graduate School of Arts and Sciences
The University of Tokyo
Meguro-ku, Tokyo 153-8902, Japan
& Kunihiko Kaneko
The Niels Bohr Institute
University of Copenhagen
Blegdamsvej 17, Copenhagen, 2100-DK, Denmark

###### Abstract

Various animals, including humans, have been suggested to perform Bayesian
inferences to handle noisy, time-varying external information. In performing
Bayesian inference, the prior distribution must be shaped by sampling noisy
external inputs. However, the mechanism by which neural activities represent
such distributions has not yet been elucidated. In this study, we demonstrated
that the neural networks with modular structures including fast and slow
modules effectively represented the prior distribution in performing accurate
Bayesian inferences. Using a recurrent neural network consisting of a main
module connected with input and output layers and a sub-module connected only
with the main module and having slower neural activity, we demonstrated that
the modular network with distinct time scales performed more accurate Bayesian
inference compared with the neural networks with uniform time scales. Prior
information was represented selectively by the slow sub-module, which could
integrate observed signals over an appropriate period and represent input
means and variances. Accordingly, the network could effectively predict the
time-varying inputs. Furthermore, by training the time scales of neurons
starting from networks with uniform time scales and without modular structure,
the above slow-fast modular network structure spontaneously emerged as a
result of learning wherein prior information was selectively represented in
the slower sub-module. These results explain how the prior distribution for
Bayesian inference is represented in the brain, provide insight into the
relevance of modular structure with time scale hierarchy to information
processing, and elucidate the significance of brain areas with slower time
scales.

_Keywords_ Recurrent Neural Network $\cdot$ Bayesian inference $\cdot$ Neural
dynamics

##  1 Introduction

In the human and various animal brain, information processing involves
inference based on inputs from the external world through the sensory systems,
which obtains information with uncertainty due to noise. Previous studies
suggested that animals such as humans and monkeys process inputs according to
a Bayesian inference framework to deal with such uncertainty(Knill and Pouget,
2004; Angelaki et al., 2009; Haefner et al., 2016; Ernst and Banks, 2002;
Friston, 2012; Merfeld et al., 1999; Doya et al., 2007; Pouget et al., 2013;
Beck et al., 2011; Geisler and Kersten, 2002; Honig et al., 2020).

Bayesian inference is performed by calculating the posterior from the prior,
which refers to the information possessed in advance about the signal, and the
likelihood estimated by observing the input signal. Hence, it is believed that
the prior must first be represented in the brain, but how prior information is
shaped in the brain remains unclear. In previous studies, the prior has often
been treated as a given value(Echeveste et al., 2020), and the mechanism for
shaping the prior by learning has not been considered. Evolutionary
acquisition of the prior has been proposed(Campbell, 2016; Lo and Zhang,
2021), whereas it is naturally expected that such information should be shaped
within one generation through observing and learning time-dependent signals.
Experimental results suggest that the prior and the likelihood for Bayesian
inference are encoded in different brain areas(Vilares et al., 2012; Chan et
al., 2016; d’Acremont et al., 2013). Still, the validity and the mechanisms
underlying the results remain controversial, and how area differentiation is
relevant to the accuracy of Bayesian inference is not well understood. A
simulation(Quax et al., 2021) suggested that a gain of the activation function
encodes the prior. However, because the prior was fixed in this study, how
shaping occurs when the prior varies over time was not considered.

In general, to obtain the prior, it is necessary to estimate the prior
distribution based on previous observations, and the population of neurons
that represents the prior must integrate observed inputs over time. One
possible mechanism for achieving such integration may be two neural modules
functioning at distinct time scales: a downstream neuron population with
slower activity changes separated from an upstream neuron population that
processes input information. In this structure, the slow module that does not
directly receive inputs may facilitate integration. Some experimental reports
have suggested that the time scale of neural activities in downstream areas of
the brain that do not directly receive external input is slow(Murray et al.,
2014; Cavanagh et al., 2020; Golesorkhi et al., 2021). On this basis, we
evaluated recurrent neural networks (RNNs) with two modules; a main module
with direct connection to the input-output layer and a sub-module with a
direct connection to the main module and no connections to the input-output
layer (i.e., a hierarchical structure)(Fig.1). Then, we examined the role of
modular structure and the relevance of the time scale difference between the
main and sub-modules for the prior representation for Bayesian inference.

We found that RNNs with a modular structure shape the prior more appropriately
than regular RNNs. Further, Bayesian inference is more accurate when the time
scale of the sub-module is appropriately slow. When the time scale is uniform,
prior information is maintained in both the main module and sub-module. In
contrast, when the time scales are different, prior information is represented
by the slow sub-module. Comparing these two cases revealed that the coded
variance of prior on the neural manifold was easier to decode in the time
scale difference model, which facilitated the distinction of the average input
change from noise.

In addition, we examined if the modular structure with distinct time scales
would emerge from a homogeneous neural network. We trained the network in a
Bayesian inference task where the time scale of each neuron varied in time. As
the training progressed, we observed that the time scales of neurons
differentiated into slower and faster scales. A modular structure arose in
which slow neurons were separated from the input/output layers, which were
predominantly connected to the fast neurons, and a sub-module with slow
neurons represented the prior information.

These results are crucial for understanding the prior representation mechanism
in Bayesian inference and provide insight into the relationships between
neural network structure, neural dynamics(Amunts et al., 2022; Mastrogiuseppe
and Ostojic, 2018; Vyas et al., 2020; Beiran et al., 2021), and time
scales(Papo, 2013) underlying information processing in the brain, which is
considered the central issue of computational neuroscience.

##  2 Model

###  2.1 Recurrent Neural Networks with/without modular structure

To investigate the effect of structure and time scale on Bayesian inference,
we considered the following RNNs(Barak, 2017).

First, we established a regular RNN consisting of an input layer, a
recurrent(hidden) layer, and an output layer, as shown in Fig.1(a). The
following equation represents the dynamics of the recurrent layer:

${\bf x}(t+1)=({\bf I}-{\boldsymbol{\alpha}}){\bf
x}(t)+{\boldsymbol{\alpha}}{\rm ReLU}(W_{in}{\bf u}(t)+W{\bf
x}(t))+\sqrt{\boldsymbol{\alpha}}{\boldsymbol{\xi}},$ (1)

where
${\boldsymbol{\alpha}}=(\alpha_{1},\alpha_{2},...,\alpha_{200})^{\mathsf{T}}$
represents a vector to introduce the time scale of the neurons as

Figure 1: Schematic of RNN. (a) Standard RNN without modular structure (b) RNN
with modular structure

$\alpha_{i}=\left\\{\begin{array}[]{ll}\alpha_{m}&(1\leq i\leq 150)\\\
\alpha_{s}&(150\leq i\leq 200),\\\ \end{array}\right.$ (2)

where the standard homogeneous network is given by $\alpha_{s}=\alpha_{m}$;
the case with $\alpha_{s}<\alpha_{m}$ was also studied to investigate the
effect of time scale difference. Although we mainly studied the systems with
150 fast, and 50 slow neurons, the results to be discussed are not altered, as
long as both the numbers are sufficient (say 100 vs 50, 150 vs 150 for fast
and slow neurons). Here, ${\bf u}(t)$ is the input signal, and ${\bf x}$ is
the state of the neurons in the recurrent layer. We adopted the activation
function ReLU(${\rm RELU}(z)=0$ for $z\leq 0$ and $=0$ for $z>0$)(Nair and
Hinton, 2010). Then, the output of the RNN was determined by the linear
combination of the internal states as follows.

${\bf y}(t)=W_{out}{\bf x}(t)$ (3)

In Eq.1, ${\boldsymbol{\xi}}$ was used to account for noise in dynamics given
by a random variable that follows a normal distribution with mean $0$ and
standard deviation $0.05$.

Next, we introduced a modular structure to the above RNN to ensure the
distinction of main and sub-modules(Fig.1(b)). Only the main module was
connected to the input/output layers. Thus, the dynamics of the recurrent
layer are given by

${\bf x}_{m}(t+1)=(1-\alpha_{m}){\bf x}_{m}(t)+\alpha_{m}{\rm ReLU}(W_{in}{\bf
u}(t)+W_{main}{\bf x}_{m}(t)+W_{s\rightarrow m}{\bf
x}_{s}(t))+\sqrt{\alpha_{m}}{\boldsymbol{\xi}}_{m}$ (4) ${\bf
x}_{s}(t+1)=(1-\alpha_{s}){\bf x}_{s}(t)+\alpha_{s}{\rm ReLU}(W_{sub}{\bf
x}_{s}(t)+W_{m\rightarrow s}{\bf
x}_{m}(t))+\sqrt{\alpha_{s}}{\boldsymbol{\xi}}_{s},$ (5)

where ${\bf x}_{m}and{\bf x}_{s}$ represent the firing rate of neurons in the
main and sub-modules, respectively. Here, $\alpha_{m}and\alpha_{s}$ represent
the time scale of the main and the sub-module, respectively. $\alpha_{m}$ is
fixed at $1$, while we varied $\alpha_{s}$ from $1$ to $0.01$ to examine the
effect of the time scale difference. The RNN output was determined by the
linear combination of internal states of the main module.

${\bf y}(t)=W_{out}{\bf x}_{m}(t)$ (6)

###  2.2 Task

In this study, we considered a task in which Bayesian inference improves
estimation accuracy. Specifically, the RNN was tasked with estimating the true
value from an observed signal with noise. We generated the external input as
follows: First, the true value $y_{true}$ was randomly sampled from a
generator(cause) distribution, given by the normal distribution with mean
$\mu_{g}$ and variance $\sigma_{g}^{2}$. Next, the observed signal $s$ was
generated from $y_{true}$ by adding noise so that the input is given by the
normal distribution with mean $y_{true}$ and variance $\sigma_{s}^{2}$. The
generator did not remain constant: It changed with probability $p_{t}$ over
time. When the generator changed, $\mu_{g},\sigma_{g}$ were sampled uniformly
from $\mu_{g}\in[-0.5,0.5],\sigma_{g}\in[0,0.8]$ respectively.

As mentioned in the Introduction, the prior distribution needed for Bayesian
estimation must be estimated from the observed signal so that it is close to
the generator distribution. Then, ${\bf u}(t)$ for Eq.1 (or 4,5) is given by
using the Probabilistic Population Code (PPC), which has been proposed as the
neural basis for Bayesian inference(Ma et al., 2006). PPC assumes that the
information in a signal is encoded by a population of neurons with a position-
based preferred stimulus that fires probabilistically according to a Poisson
distribution. It has been shown that neural networks with a population of
neurons following PPC as the input layer can learn probabilistic inference
effectively(Orhan and Ma, 2017). Therefore, in this study, we also assumed
that the activity ${\bf u}$ of the input-layer neurons encoding the observed
signal followed the PPC model. ${\bf u}$ was sampled from the following
Poisson distribution(Ichikawa and Kataoka, 2022):

$p({\bf u}|s)=\prod_{i}\frac{e^{-f_{i}(s)}f_{i}(s)^{u_{i}}}{u_{i}!}$ (7)

Here, $s$ is the observed signal generated from $y_{true}$ by adding noise,
and $f_{i}$ is the tuning curve of the neurons. This selective firing occurs
in proportion to the gain when the observed signal is generated. This gain is
inversely proportional to the noise variance as $g=1/\sigma_{l}^{2}$, and
corresponds to signal clarity. Namely, the gain decreases and noise increases
due to uncertainty in observations(Tolhurst et al., 1983). Considering the
gain, we obtain:

$f_{i}(s)=g\exp\biggl{(}\frac{-(s-\phi_{i})^{2}}{2\sigma_{\rm
PPC}^{2}}\biggr{)},$ (8)

where $\phi_{i}$ represents the preferred stimuli of neurons in the input
layer. It was assumed that $\phi_{i}$ follows an arithmetic sequence for $i$
($\phi_{i}=-1/2+i/m$ when the number of neurons in the input layer is
$m$)(Swindale, 1998). Also, $\sigma_{\rm PPC}^{2}$ is a constant that
represents the ease of firing and was set as $\sigma_{\rm PPC}^{2}=1/2$ in
this study.

In this task, the true value $y_{true}$ was to be estimated based on the input
signal ${\bf u}$. Therefore, training was performed to minimize the mean
squared error (MSE) between the neural network output $y(t)$ and the true
value $y_{true}(t)$. Note that the loss function was not based on the Bayesian
optimal value calculated from the generator distribution and the noise in the
observed signal but only calculated based on the true value.

$L=\frac{1}{T}\sum_{t}(y(t)-y_{true}(t))^{2}$ (9)

Training was performed by the backpropagation method (Rumelhart et al., 1986;
Werbos, 1990). An efficient Stochastic Gradient Descent method, Adam(Kingma
and Ba, 2014), was used for optimization. The batch size of training samples
was set to 50, and the weight decay rate was set to 0.0001; training was
performed for 6000 iterations(See Table.1 for the hyperparameters used in the
experiment).

Table 1: Hyperparameters Attribute | Value
---|---
Range of $\mu_{p}$ | $-0.5\leq\mu_{p}\leq 0.5$
Range of $\sigma_{p}$ | $0\leq\sigma_{p}\leq 0.8$
Range of $\sigma_{l}$ | $\sqrt{1/5}\leq\sigma_{l}\leq 1$
Switching probability of prior | $p_{t}=0.03$
Length of ${\bf u}(t)$ | $100$
$\sigma_{\rm PPC}$ | 0.5
Lasting time of ${\bf u(t)}$ | $T=120$
#Neurons in the main module | 150
#Neurons in the submodule | 50
$\alpha_{m}$ | 1
$\alpha_{s}$ | 1, 0.5, 0.2, 0.1, 0.05, 0.01
Batch size | 50
Optimization algorithm | Adam
Learning rate | 0.001
Iteration | 6000
Weight decay | 0.0001

## Results1: Fixed structure and time scales

### Bayesian optimality

Because the generated signal $s$ was observed under noise, the neural network
was required to estimate the true value sampled from the generator. If the
information from the generator was known, $y_{true}$ would be estimated by
minimizing the long-term MSE, which reveals the optimal $y$ value as follows
(maximum a posteriori(MAP) estimation (Bishop, 2006)).

$y_{opt}=\frac{\sigma_{g}^{2}}{\sigma_{g}^{2}+\sigma_{s}^{2}}s+\frac{\sigma_{s}^{2}}{\sigma_{g}^{2}+\sigma_{g}^{2}}\mu_{g}$
(10)

However, as described in the “Task” section, the information from the
generator was not explicitly given to the neural network, so it must be
estimated from observed signals as a prior distribution. First, we examined
whether the neural network could achieve this prior-based estimation.

The output $y$ of RNN with modular structure trained with $\alpha_{s}=0.1$,
when given an observed signal $s$, is shown Fig.2. $s$ was sampled from the
prior with $\mu_{g}=0.5,\sigma_{g}=0.5$, and $\sigma_{s}=\sqrt{1/5}$ of noise
was added. The green points represent the estimation based on the maximum
likelihood estimation $y_{ML}$, which is that with the highest accuracy when
no prior information is available. Here, this estimation is equal to the
observed signal $s$. The blue points represent $y_{opt}$ when estimated
according to the MAP estimation, and the orange points represent the actual
neural network output $y$. Fig.2 shows that the output of RNN is closer to the
blue points $y_{opt}$ rather than to the green points, indicating that
approximate Bayesian inference (Near-optimal Bayesian inference) with a well-
estimated prior is achieved (the mean squared error between $y$ and $y_{ML}$
is $0.15$, and the mean squared error between $y$ and $y_{opt}$ is $0.019$,
the latter being smaller).

Figure 2: The output $y$ of RNN against the observed signal value $s$. Before
$s$ is input, the time series signal, which is sampled from the normal
distribution with the mean $\mu_{g}=0.5$ and the standard deviation
$\sigma_{g}=0.5$ and then the noise with the standard deviation
$\sigma_{s}=\sqrt{1/5}$ is added in the input. The accuracy can be increased
by estimating prior based on the signal input before $s$ and performing
Bayesian inference. Blue points represent
$y_{opt}=\frac{\sigma_{g}^{2}}{\sigma_{g}^{2}+\sigma_{s}^{2}}s+\frac{\sigma_{s}^{2}}{\sigma_{g}^{2}+\sigma_{g}^{2}}\mu_{g}$
value, orange points represent the output of RNN $y$, and green points
represent estimation based on maximum likelihood estimation $y_{ML}=s$. The
result is for a model with $\alpha_{s}=0.1$.

Next, we examined the optimality of the Bayesian estimation for networks with
and without modular structures and time scale differences. Fig.3(a) shows the
MSE between $y$ and $y_{opt}$ by the RNN trained under each condition. This
result shows that the modular structure improved the accuracy of Bayesian
estimation, which was further increased when $\alpha_{s}$ decreased to an
appropriate degree. In fact, we found the optimal time scale
$\alpha_{s}=0.1\sim 0.2$, at which maximum accuracy was achieved. Even without
modular structure, the time scale difference contributed to inference
accuracy, but the accuracy increased significantly with both the modular
structure and time scale difference.

Figure 3:  (a) MSE between the optimal value $y_{opt}(t)$ and the output of
RNN $y(t)$, plotted against the time scale $\alpha_{s}$. $\cdot$ with modular
and $\times$ without modular structure. RNNs with a modular structure is more
accurate. In addition, those with $\alpha_{s}\sim(0.1\sim 0.2)$ have optimal
error. (b) MSE between the true value $y_{true}(t)$ and the output of RNN
$y(t)$ for the network with $\alpha_{s}=1(\cdot)$ and
$\alpha_{s}=0.1(\times)$. The value increases as $p_{t}$ increases, but the
model with $\alpha_{s}=0.1$ is always more accurate.

### Adjustability to rapid generator switching

So far, we studied the performance of Bayesian inference models under a fixed
generator to compare the accuracy of Bayesian inference itself. Next, we
examined their performance when the generator changes in time. To perform
Bayesian inference for a rapidly changing input, it was necessary for the
model to quickly approach the new optimal value $y_{opt}$ to yield a good
estimation. To verify the accuracy of the RNN in this case, we compared the
MSE between $y_{true}(t)$ generated by the generator and the output $y(t)$ of
RNN under various $p_{t}$(Fig.2(b)). The model with $\alpha_{s}=0.1$ was found
to be more accurate for all values of $p_{t}$.

As a special case, we considered a setting where the input moves back and
forth between two generators, A and B. Then we examined whether the prior
distribution estimated by the RNN was closer to the distribution of either
generator. Specifically, we adopted the generator A with
$(\mu_{g},\sigma_{g}^{2})=(\mu_{A},\sigma_{A}^{2})$ and the generator B with
$(\mu_{g},\sigma_{g}^{2})=(\mu_{B},\sigma_{B}^{2})$ and computed the following
values when the Bayesian optimal estimates under each generators were
$y_{opt}^{A},y_{opt}^{B}$.

$a(t)=\frac{y_{opt}^{B}-y(t)}{y_{opt}^{B}-y_{opt}^{A}}$ (11)

When $a(t)$ is close to 1, the model’s prior is closer to generator B, and
when $a(t)$ is close to -1, it is closer to prior A.

Comparing the change in $a(t)$ between the model with $\alpha_{s}=0.1$ and the
model with $\alpha_{s}=1$, we found that the model with $\alpha_{s}=0.1$ was
more adjustable to the generator change as shown in Fig.4(a). This result
shows that the model with $\alpha_{s}=0.1$ was more responsive to the changes
of the generators and recognized the generator change more quickly in all
runs. The difference between the two models was especially pronounced in the
extreme case in which the two generators switched every time(Fig.4(b)).
Intuitively, having a population of slow neurons would seem to be a
disadvantage in responding to rapid environmental changes, but the results
showed the opposite. The network with $\alpha_{s}=1$ could not follow rapid
input changes, whereas that with $\alpha_{s}=0.1$ could estimate the input
prior effectively. We discuss the importance of slow neurons in responding to
rapid changes below.

Figure 4:  Adjustability to rapid generator change. (a) $a(t)$ for the case
where generator A($(\mu_{A},\sigma_{A}^{2})=(-0.5,0.04)$) and generator
B($(\mu_{A},\sigma_{A}^{2})=(0.5,0.04)$) switch alternately with probability
$p_{t}=0.2$. The model with $\alpha_{s}=0.1$ adjusted more quickly to the
generator change. The thin line represents $a(t)$ when the output was fully
adjusted to generator switching. (b) $a(t)$ for the case where generator
A($(\mu_{A},\sigma_{A}^{2})=(-0.5,0.04)$) and generator
B($(\mu_{A},\sigma_{A}^{2})=(0.5,0.04)$) switch every time(periodic switch).

### Representation of the prior

We investigated how the slow sub-module facilitated improved prior
representation for Bayesian inference. Beginning from the hypothesis that a
group of downstream slow neurons represent the prior by integrating the
observed signal over time, we investigated which side of the main/sub-module
was responsible for the prior information in the modular RNN.

Here, by using the prior information, the estimated value was shifted from the
observed signal $s$ to an appropriate value $y_{opt}$(Eq.10). In other words,
even given the same signal input $s$, the output varied depending on which
time series signal was input before $s$ (because the prior estimation
changed). Even if one module returned to its original state, the output
shifted from $s$ because the prior information remained in the other module.
The scale of this change is considered to represent the degree to which the
module utilizes the estimated prior information. Therefore, it is possible to
estimate the extent to which each module plays a role in prior information
processing by examining the change in the output $y(t)$ when the internal
state of each main and sub-modules is changed to the value corresponding to a
different prior.

First, let ${\bf x}_{m}(\mu_{g},\sigma_{g}),{\bf x}_{s}(\mu_{g},\sigma_{g})$
be the internal states of the main and sub-module, respectively when the input
signal $s$ from a generator $(\mu_{g},\sigma_{g})$ is applied for a certain
period. Because the output $y$ is determined by the internal states of two
modules and the input signal, it can be written as $y({\bf
x}_{m}(\mu_{g},\sigma_{g});{\bf x}_{s}(\mu_{g},\sigma_{g}),s,\sigma_{s})$.
From this, the change in output $y$ is computed by fixing one of the two
modules and varying the other to a different internal state ${\bf
x}_{i}(\mu_{g},\sigma_{g})\rightarrow{\bf
x}_{i}(\mu_{g}^{\prime},\sigma_{g}^{\prime})$. The degree of change in $y$
represents the impact on the output of each module reflecting the prior
information. Hence, by comparing the above variances of $y$ by ${\bf
x}_{m}$(or ${\bf x}_{s}$) with fixed ${\bf x}_{s}$(or ${\bf x}_{m}$)
respectively, it is possible to estimate how much each module is responsible
for the prior representation. Specifically, we fixed one of the modules at
$\mu_{g}=0$, $\sigma_{g}=0.4$ (These values are set to the median of the range
of values $-0.5\leq\mu_{g}\leq 0.5$, $0\leq\sigma_{p}\leq 0.8$), i.e. ${\bf
x}_{i}(0,0.4)$, while for the other module $\mu_{g}$ and $\sigma_{g}$ are
changed as ${\bf x}_{g}(\mu_{g},\sigma_{g})$. Then, we calculated the variance
of $y$ as

$V_{m}=\langle{\rm Var}[y({\bf x}_{m}(\mu_{p},\sigma_{p}),{\bf
x}_{s}(0,0.4),s,\sigma_{l})]_{(\mu_{p},\sigma_{p})}\rangle_{(s,\sigma_{l})}$
(12) $V_{s}=\langle{\rm Var}[y({\bf x}_{m}(0,0.4),{\bf
x}_{s}(\mu_{p},\sigma_{p}),s,\sigma_{l})]_{(\mu_{p},\sigma_{p})}\rangle_{(s,\sigma_{l})},$
(13)

where ${\rm Var}[\ ]_{(\mu_{p},\sigma_{p})}$ denotes the variance over the
changes of $(\mu_{p},\sigma_{p})$, and $\langle\rangle_{(s,\sigma_{l})}$
denotes the average over the changes of $(s,\sigma_{l})$. The larger $V_{s}$
or $V_{m}$ indicates that the sub-module or main module strongly reflects the
difference in the prior distribution to the difference in output,
respectively.

Dependencies of $V_{s}$ and $V_{m}$ on different $\alpha_{s}$ are shown in
Fig.5. This result shows that when $\alpha_{s}=1$ (i.e., the time scale is
uniform), both the main and sub-modules contribute to the representation of
prior distribution to the same degree. Conversely, when $\alpha_{s}=0.1\sim
0.5$, $V_{s}$ is much larger than $V_{m}$, meaning that the sub-module
selectively contributes to the representation of the prior. In particular,
when $\alpha_{s}=0.1$ and $0.2$, the differentiation of representation between
the main and sub-modules is more pronounced. Note that the contribution of the
main module is large when $\alpha_{s}=0.01$, probably because the time scale
of the sub-module is too slow to code the information of the prior. Comparing
of Fig.5 and Fig.3 shows that the highly accurate Bayesian inference is
achieved when the prior distribution information is localized in the sub-
module.

Figure 5:  Division of roles for representing prior distribution.
$V_{s},V_{m}$ defined in the text Eqs. (12,13) plotted for different values of
$\alpha_{s}$ computed over $1000$ samples of data. $V_{s}$ and $V_{m}$
represent the degree to which the sub-module and the main module are
responsible for prior-based information processing. When $\alpha_{s}=0.2,0.1$,
the sub-module selectively contributes to the representation of the prior.

Next, we investigated how the prior is represented by the main and sub-modules
by visualizing the neural activity by principal component analysis(PCA)(Mante
and et al., 2013; Ichikawa and Kaneko, 2021). First, ${\bf
x}_{m}(\mu_{g},\sigma_{g})$ and ${\bf x}_{s}(\mu_{g},\sigma_{g})$ were
computed for various $(\mu_{g},\sigma_{g})$ in a model with $\alpha_{s}=0.1$,
and made PCA. The results were projected on a plane using the first and second
principal components and color-coded according to $\mu_{g}$ and $\sigma_{g}$
(Fig.6(a,b)). The neural activity in the main module was loosely distributed
on a one-dimensional manifold, represented by the first principal
component(PC1). This PC1 approximately corresponded to the $\mu_{g}$ value,
although the distinction was not clear. In contrast, the activity in the sub-
module was clearly represented by 2-dimensional manifolds, as in Fig.6(b2),
where PC1 corresponds to $\mu_{g}$ and PC2 corresponds to $\sigma_{g}$, rather
well.

Figure 6:  The neural activities of the main module ${\bf
x}_{m}(\mu_{g},\sigma_{g})$ and sub-module ${\bf x}_{s}(\mu_{g},\sigma_{g})$
were plotted by the first and second principal component spaces. (a,b) is the
result of $\alpha=0.1$, and (c,d) is of $\alpha_{s}=1$. (a1,c1) Main module,
color-coded by $\mu_{g}$. (a2, c2) Main module, color-coded by $\sigma_{g}$.
(b1,d1) Sub-module, color-coded by $\mu_{g}$. (b2, d2) Sub-module, color-coded
by $\sigma_{g}$. 300 data are plotted.

Then, we performed the same analysis on the model with $\alpha_{s}=1$
(Fig.6(c,d)). In this case, the manifolds of neural activities for the main
and sub-modules did not change significantly. Both were represented in a one-
dimensional manifold corresponding to $\mu_{g}$; there was no axis
corresponding to $\sigma_{g}$. The decodability of $\sigma_{g}$ achieved in
the internal states of sub-module with $\alpha_{s}=0.1$ was not observed for
$\alpha_{s}=1$. In fact, the coefficient of determination when $\sigma_{g}$
was calculated by Ridge regression from the internal state of the sub-module
with $\alpha_{s}=0.1$ was $0.68$, while that using the sub-module with
$\alpha_{s}=1$ is $-0.03$. This suggests that the model with $\alpha_{s}\sim
0.1$ can better distinguish the input’s variance from noise to accurately
perform Bayesian inference.

When the generator changed rapidly, the variance of the prior was larger than
the variance of the generator, as shown in the SI for the case with
$\alpha_{s}=0.1$. When $\sigma_{g}$ was large, as seen from Eq.10, the
influence of the observed signal $s$ was larger than that of $\mu_{g}$,
allowing the model to "keep up" with large changes in the observed signal.
This explains the higher adjustability to rapid generator changes as seen in
Fig.4.

### Effects of different time scales

To examine the impact of $\alpha_{s}$ differences on Bayesian inference
accuracy in detail, we considered how each model with $\alpha_{s}=1$ and
$\alpha_{s}=0.1$ represents prior as a function of the input signal. As seen
in Fig.6, when the generator is constant, the internal state of the neural
network corresponds with the state of the generator $(\mu_{g},\sigma_{g})$.
Conversely, when the generator changes, the internal state at a certain time
does not necessarily correspond to the state of the generator at that time
because some time is needed to estimate the state of the prior after the
generator switches. Let $\mu_{p}$ and $\sigma_{p}$ be the mean and standard
deviation of prior used by the neural network to compute $y(t)$. Hence,
$\mu_{p}$ must memorize the input $s(t)$ for a certain time in the form of

$\mu_{p}\simeq\sum a_{k}s(t-k).$ (14)

To examine how many past steps $k$ are memorized, $\mu_{p}$ must be estimated.
This can be achieved by estimating $\mu_{p}$ from the internal state ${\bf
x}(t)$.

First, we calculated the internal state ${\bf x}$ for the observed signal with
a fixed generator instead of a time-varying case. Then, we found the
transformation matrix $W_{\mu_{p}}$ from the internal state ${\bf x}$ to the
recognized prior $\mu_{p}$ by assuming that $\mu_{p}$ can be represented by
linear transformations of the internal state as $\mu_{p}\simeq W_{\mu_{p}}{\bf
x}$. This transformation matrix $W_{\mu_{p}}$ was obtained by a pseudo-inverse
method(Schrauwen et al., 2007)(SI).

Next, we obtained ${\bf x}(t)$ against the time-varying signal with a
probability of $p_{t}=0.03$. By applying the above transformation matrix,
$W_{\mu_{p}}$ to ${\bf x}(t)$ obtained at this time, the prior $\mu_{p}$ was
estimated accordingly. The state of the prior was thus obtained for the time
series of the observed signal $s(t)$.

Then, $a_{k}$ in Eq.14 was obtained to minimize the difference between the two
sides of Eq.14. Because the obtained coefficients correspond to the
contribution of the signal before $k$ time steps, we could estimate the extent
to which the neural network uses past information when estimating the prior.

Figure 7:  $a_{k}$ defined by Eq.14 is plotted against $t$, for the model with
$\alpha_{s}=1$ and $\alpha_{s}=0.1$ using 3000 data points.

The estimated coefficients of Eq.14 were plotted against $k$(Fig.7), revealing
that the model with $\alpha_{s}=0.1$ used more past information in estimating
prior information than the model with $\alpha_{s}=1$. This difference in time
windows leads to a difference in accuracy for prior encoding.

## Results2: Modular structure organization and time-scale separation by
learning

So far, we investigated neural networks with fixed and modular structures
along fixed time scales and demonstrated that those with fast and slow modules
effectively represented the prior distribution. Then, we investigated whether
such a structure would emerge by training a neural network to predict
$y_{true}$. We again used the same neural network model as the normal RNN.

${\bf x}(t+1)=({\bf I}-{\boldsymbol{\alpha}}){\bf
x}(t)+{\boldsymbol{\alpha}}{\rm ReLU}(W_{in}{\bf u}(t)+W_{rec}{\bf
x}(t))+\sqrt{\boldsymbol{\alpha}}{\boldsymbol{\xi}},$ (15)

where ${\boldsymbol{\alpha}}$ represents a vector of time scales of neurons
consisting of $\alpha_{i}$. These ${\boldsymbol{\alpha}}$ values, as well as
elements of $W$, change by training to start from initial values set randomly
according to $\mathcal{N}(0.5,0.1)$. During training, each matrix $W$ and
${\boldsymbol{\alpha}}$ are optimized according to the gradient descent
method(Perez-Nieves et al., 2021) at each step. The number of neurons in the
recurrent layer of the neural network was set to 80.

The change in ${\boldsymbol{\alpha}}$ distribution during the learning task is
shown in Fig.8(a). As shown, ${\boldsymbol{\alpha}}$ split into two groups
over the learning period: one with large values close to 1 and the other with
small values near 0.1.

Next, we measured the contribution of prior representation as examined in the
"Representation of the prior" section for groups of neurons with large values
${\boldsymbol{\alpha}}$ (neurons with $\alpha_{i}>0.8$) and groups of neurons
with small values ${\boldsymbol{\alpha}}$(neurons with $\alpha_{i}<0.2$) for
three epochs in the learning process(Fig.8(b)). We found that after 10000
epochs, the slow neurons were responsible for the representation of prior
distribution, as in the model with $\alpha_{s}=0.1$ in the fixed time scale
setting.

Figure 8:  RNN features obtained by learning when ${\boldsymbol{\alpha}}$ is
variable by learning. (a) Frequency distribution of ${\boldsymbol{\alpha}}$
for all neurons at 200, 1000, and 10000 learning epochs. At 10000 epochs, the
learning process was complete. (b) Division of roles $V_{slow}$ and
$V_{fast}$. See the "Representation of the prior" section for definitions of
$V_{slow}$ and $V_{fast}$; $V_{slow}\ (V_{fast})$ was computed for neurons
with $\alpha<0.2\ (\alpha>0.8)$ respectively. (c) The average degree of RNN
connections of 10000 epochs. Connections between input layer, recurrent
neurons with $\alpha<0.2$, $0.2\leq\alpha\leq 0.8$, $\alpha>0.8$, and output
layer. Each was normalized so that the maximum value was 1.

Finally, we investigated the neural network structure shaped by training. In
Fig.8(a), the recurrent layer neurons of the network of epoch 10000 was split
into the three groups, divided by the magnitude of $\alpha_{i}$, slow neurons
with $\alpha_{i}<0.2$, fast neurons with $\alpha_{i}>0.8$, and
$0.2\leq\alpha_{i}\leq 0.8$ neurons as the others. The average connectivity
between the input layer, each group, and the output layer is shown in
Fig.8(c)(Yang et al., 2019). The connection from the input layer to the group
of fast neurons and that from the fast neurons to the output layer were
distinctively larger than those to or from the slow neurons. Among connections
within the recurrent layer, those between the fast and slow neurons were
larger than others. In summary, a modular structure, shown in Fig.1(b),
emerged through learning alone.

## Discussion

In this study, we demonstrated that neural networks with slow and fast
activity modules play an essential role in the prior representation for
Bayesian inference. We set up a task to predict a time-varying signal under
noise that could be estimated by Bayesian inference and trained RNNs with or
without modular structure and with or without time scale differences.

The RNN could learn to approximate Bayesian inference using
prior(approximating the generator distribution) in all conditions tested.
However, the accuracy was higher in the modular RNN; further, the accuracy was
significantly higher when the time scale of the sub-module was moderately
slower than that of the main module. In addition, the increase in accuracy was
pronounced against a rapidly varying input, for which it was necessary to
generate a prior that changes quickly. To achieve such accuracy with a slow
sub-module, the sub-module was found to dominantly represent the prior,
indicating role differentiation between representation of the prior and
representation of the observed signal (likelihood). Of note, such functional
differentiation is caused by differences in time scales. This result is
consistent with experimental observations in the brain in which areas that
code the prior and likelihood in Bayesian inference are different(Vilares et
al., 2012; Chan et al., 2016; d’Acremont et al., 2013). Finally, it was shown
that a modular structure with distinct time scales was spontaneously organized
in the RNN by learning.

It is important to note that a relatively slow time scale of the neuron
population encoding the prior is required, but the difference between fast and
slow neurons should not be excessive. If the time scale is too small, the
accuracy is decreased (Fig.3) in which case the sub-module is not responsible
for representing the prior (Fig.5). This is because prior construction
requires a larger time span to address changes in external input for a neural
network with such a slow time scale. Therefore, we suggest that there is an
optimal time scale for the slow sub-module. Future research should investigate
how this optimal time scale depends on the time scale of environmental
changes.

It has been suggested that the time scale of neurons slows down hierarchically
from the area where the signal is directly applied to the area where
information is proceed(Murray et al., 2014; Cavanagh et al., 2020; Golesorkhi
et al., 2021). This time scale hierarchy with a modular structure(Yamashita
and Tani, 2008) is suggested to be relevant to information processing(Kurikawa
and Kaneko, 2021; Yamashita and Tani, 2008; Tanaka et al., 2022). Our study
showed that modular structures with two-level time scales could deal with
slowly changing inputs. A deeper modular structure with multiple time scales
may be necessary to deal with further complex changes in environments. With
such a structure, Bayesian inference against complex temporal changes could be
achieved by extrapolating the results of this study. Further research
verifying this finding will elucidate the significance of hierarchical
structuring in the brain. It is noteworthy that the time scale separation was
not only found to be influential for accurate Bayesian inference but also
emerged from learning in our simulation. Considering these findings, a similar
process may be expected in evolution(Yamaguti and Tsuda, 2021).

The modular network with slow/fast time scales could integrate out noise and
distinguish the average change in the inputs from fast noise. In fact, the
network could effectively predict temporal changes in the input, even under
rapidly changing conditions. The brain must adapt to time-varying, noisy
inputs; hence, the performance of Bayesian inference by the network design
reported herein is considered relevant to brain information processing.

We adopted a simple RNN and trained it using backpropagation.
redBackpropagation is often believed to be different from the learning
algorithm implemented in the brain(Bengio et al., 2015; Lillicrap et al.,
2016), so care should be taken when generalizing our results. However,
previous studies have also suggested that neural networks obtained by
backpropagation can show similar behavior to that of the actual brain(Richards
et al., 2019; Yang and Wang, 2020; Mante and et al., 2013; Barak et al., 2013;
Cueva and Wei, 2018; Yamins and DiCarlo, 2016; Haesemeyer et al., 2019). With
these considerations, our findings are considered to be relevant to the
brain’s learning processes despite the potential limitation of
backpropagation.

Unraveling the relationship between the structure of neural networks, neural
dynamics, and the information processing performed by the brain is a primary
goal in computational neuroscience(Mastrogiuseppe and Ostojic, 2018; Dubreuil
et al., 2022; Vyas et al., 2020; Beiran et al., 2021). In this study, the
relevance of modular structure and time scale difference in neural dynamics to
the representation of the prior in Bayesian inference is demonstrated, as well
as their formation by learning(Lorenz et al., 2011; Kashtan and Alon, 2005),
which will support ongoing research in the field.

### Data Availability

Source codes for these models can be found at
https://github.com/tripdancer0916/slow-reservoir

## Appendix

### Appendix A: Adjustment to rapid environment changes

The trajectories of the internal state ${\bf x}_{sub}(t)$ are plotted by the
first and second principal components in Fig. 6 for the cases in which
generators A and B switch every 2 time steps and every 30 time steps.
Generators A and B both have $\sigma_{g}=0.04$. In the case of switching every
30 time steps, they were located in the region taken by the internal state
when $\sigma_{g}$ was small. In the case of switching every 2 time steps, they
were located in the region taken by the internal state when $\sigma_{g}$ was
large (See Fig. S1). This occurred because the generators switched so rapidly
that the RNN recognized that the signal was created by a generator with a
large variance. This made it possible to switch $y(t)$ quickly because the
information of the observed signal $s$ was prioritized over the prior
information when calculating the output $y$.

Figure 9: Trajectory of the internal state ${\bf x}_{sub}(t)$ of the sub-
module when generator A($(\mu_{A},\sigma_{A}^{2})=(-0.5,0.04)$) and generator
B($(\mu_{A},\sigma_{A}^{2})=(0.5,0.04)$) switch alternately.

### Appendix B: Results of the RNN with $\alpha_{m}=\alpha_{s}=0.1$

We argue that the slower time scale of the sub-module relative to the main
module is important for accurate Bayesian inference. Here, to investigate
whether the difference in time scale or the slower time scale itself were more
influential, we trained an RNN with $\alpha_{s}=\alpha_{m}=0.1$ and examined
its accuracy. We found that when $(\alpha_{m},\alpha_{s})=(0.1,0.1)$, the MSE
was larger, and the accuracy was worse than that in the cases with
$(\alpha_{m},\alpha_{s})=(1,1)$ and $(\alpha_{m},\alpha_{s})=(1,0.1)$, as
shown in Fig. S2. Therefore, it is not simply the slower time scale of the
neurons but the time scale difference between the main and sub-modules that
facilitate accurate Bayesian inference.

Figure 10: Mean squared error between the optimal value $y_{opt}(t)$ and the
output of RNN $y(t)$, plotted against the setting
$(\alpha_{m},\alpha_{s})=(1,1),(1,0.1),(0.1,0.1)$.

### Appendix C: Detailed analysis of the time scale difference effect

To obtain the transformation matrix $W_{\mu_{p}}$, we created a data vector
$M_{g}$ that arranges $\mu_{g}$ and a data matrix $X$ that arranges the
internal state $\bf x$ obtained when fixed $\mu_{g}$ was input as follows:

$M_{g}=(\mu_{g}^{1},\mu_{g}^{2},...)^{\mathsf{T}},\ X=({\bf x}^{1},{\bf
x}^{2},...)$

Then, we attempted to find matrix $W_{\mu_{p}}$ such that $M_{g}\simeq
W_{\mu_{p}}X$. Using the Moore-Penrose pseudo inverse, we can find the best-
fit matrix as $W_{\mu_{p}}=M_{g}X^{\dagger}$(Penrose, 1955). Let $\mu_{p}$ be
the result of the transformation by $W_{\mu_{p}}$. As Fig. S3(a) shows,
$\mu_{g}\simeq\mu_{p}$ is valid. To find $a_{k}$, we calculated ${\bf x}(t)$
in the case that $\mu_{g}$ varies randomly using $W_{\mu_{p}}$ and obtained
$\mu_{p}$. Then, we created a data vector $M_{p}$ that arranges $\mu_{p}$ and
a data matrix $S$ that arranges ${\bf
s}=(s(t-1),s(t-2),...,s(t-K))^{\mathsf{T}}$. Using the Moore-Penrose pseudo
inverse, we found ${\bf a}=(a_{1},a_{2},...,a_{K})^{\mathsf{T}}$. As Fig.
S3(a) shows, $\mu_{p}\simeq\sum_{k}^{K}a_{k}s(t-k)$ was valid.

Figure 11: (a)Comparison between the estimated mean of prior $\mu_{p}$ and the
mean of generator $\mu_{g}$. (b)Comparison between the linear weighted sum of
past signals $s(t-k)$ and the estimated mean of prior $\mu_{p}$.

## References

  * Knill and Pouget [2004] David C. Knill and Alexandre Pouget.  The bayesian brain: the role of uncertainty in neural coding and computation.  _Trends in Neurosciences_ , 27(12):712–719, 2004\.  ISSN 0166-2236.  doi:https://doi.org/10.1016/j.tins.2004.10.007.  URL https://www.sciencedirect.com/science/article/pii/S0166223604003352.
  * Angelaki et al. [2009] Dora E Angelaki, Yong Gu, and Gregory C DeAngelis.  Multisensory integration: psychophysics, neurophysiology, and computation.  _Current Opinion in Neurobiology_ , 19(4):452–458, 2009.  ISSN 0959-4388.  doi:https://doi.org/10.1016/j.conb.2009.06.008.  URL https://www.sciencedirect.com/science/article/pii/S0959438809000725.  Sensory systems.
  * Haefner et al. [2016] Ralf M. Haefner, Pietro Berkes, and József Fiser.  Perceptual decision-making as probabilistic inference by neural sampling.  _Neuron_ , 90(3):649–660, 2016.  ISSN 0896-6273.  doi:https://doi.org/10.1016/j.neuron.2016.03.020.  URL https://www.sciencedirect.com/science/article/pii/S0896627316300113.
  * Ernst and Banks [2002] Marc O. Ernst and Martin S. Banks.  Humans integrate visual and haptic information in a statistically optimal fashion.  _Nature_ , 415(6870):429–433, Jan 2002.  ISSN 1476-4687.  doi:10.1038/415429a.  URL https://doi.org/10.1038/415429a.
  * Friston [2012] Karl Friston.  The history of the future of the bayesian brain.  _NeuroImage_ , 62(2):1230–1233, 2012.  ISSN 1053-8119.  doi:https://doi.org/10.1016/j.neuroimage.2011.10.004.  URL https://www.sciencedirect.com/science/article/pii/S1053811911011657.  20 YEARS OF fMRI.
  * Merfeld et al. [1999] Daniel M. Merfeld, Lionel Zupan, and Robert J. Peterka.  Humans use internal models to estimate gravity and linear acceleration.  _Nature_ , 398(6728):615–618, Apr 1999.  ISSN 1476-4687.  doi:10.1038/19303.  URL https://doi.org/10.1038/19303.
  * Doya et al. [2007] Kenji Doya, Shin Ishii, Alexandre Pouget, and Rajesh P N Rao.  _Bayesian Brain: Probabilistic Approaches to Neural Coding_.  MIT Press, 2007.
  * Pouget et al. [2013] Alexandre Pouget, Jeffrey M Beck, Wei Ji Ma, and Peter E Latham.  Probabilistic brains: knowns and unknowns.  _Nature Neuroscience_ , 16(9):1170–1178, 2013\.  doi:10.1038/nn.3495.  URL https://doi.org/10.1038/nn.3495.
  * Beck et al. [2011] Jeffrey M. Beck, Peter E. Latham, and Alexandre Pouget.  Marginalization in neural circuits with divisive normalization.  _Journal of Neuroscience_ , 31(43):15310–15319, 2011.  ISSN 0270-6474.  doi:10.1523/JNEUROSCI.1706-11.2011.  URL https://www.jneurosci.org/content/31/43/15310.
  * Geisler and Kersten [2002] Wilson S. Geisler and Daniel Kersten.  Illusions, perception and bayes.  _Nature Neuroscience_ , 5(6):508–510, 2002.  doi:10.1038/nn0602-508.  URL https://doi.org/10.1038/nn0602-508.
  * Honig et al. [2020] Maija Honig, Wei Ji Ma, and Daryl Fougnie.  Humans incorporate trial-to-trial working memory uncertainty into rewarded decisions.  _Proceedings of the National Academy of Sciences_ , 117(15):8391–8397, 2020.  doi:10.1073/pnas.1918143117.  URL https://www.pnas.org/doi/abs/10.1073/pnas.1918143117.
  * Echeveste et al. [2020] Rodrigo Echeveste, Laurence Aitchison, Guillaume Hennequin, and Máté Lengyel.  Cortical-like dynamics in recurrent circuits optimized for sampling-based probabilistic inference.  _Nature Neuroscience_ , 23(9):1138–1149, 2020\.  doi:10.1038/s41593-020-0671-1.  URL https://doi.org/10.1038/s41593-020-0671-1.
  * Campbell [2016] John O. Campbell.  Universal darwinism as a process of bayesian inference.  _Frontiers in Systems Neuroscience_ , 10, 2016.  ISSN 1662-5137.  doi:10.3389/fnsys.2016.00049.  URL https://www.frontiersin.org/articles/10.3389/fnsys.2016.00049.
  * Lo and Zhang [2021] Andrew W. Lo and Ruixun Zhang.  The evolutionary origin of bayesian heuristics and finite memory.  _iScience_ , 24(8):102853, 2021.  ISSN 2589-0042.  doi:https://doi.org/10.1016/j.isci.2021.102853.  URL https://www.sciencedirect.com/science/article/pii/S258900422100821X.
  * Vilares et al. [2012] Iris Vilares, James D. Howard, Hugo L. Fernandes, Jay A. Gottfried, and Konrad P. Kording.  Differential representations of prior and likelihood uncertainty in the human brain.  _Current Biology_ , 22(18):1641–1648, 2022/07/26 2012.  doi:10.1016/j.cub.2012.07.010.  URL https://doi.org/10.1016/j.cub.2012.07.010.
  * Chan et al. [2016] Stephanie C. Y. Chan, Yael Niv, and Kenneth A. Norman.  A probability distribution over latent causes, in the orbitofrontal cortex.  _Journal of Neuroscience_ , 36(30):7817–7828, 2016.  ISSN 0270-6474.  doi:10.1523/JNEUROSCI.0659-16.2016.  URL https://www.jneurosci.org/content/36/30/7817.
  * d’Acremont et al. [2013] Mathieu d’Acremont, Wolfram Schultz, and Peter Bossaerts.  The human brain encodes event frequencies while forming subjective beliefs.  _Journal of Neuroscience_ , 33(26):10887–10897, 2013.  ISSN 0270-6474.  doi:10.1523/JNEUROSCI.5829-12.2013.  URL https://www.jneurosci.org/content/33/26/10887.
  * Quax et al. [2021] Silvan C. Quax, Sander E. Bosch, Marius V. Peelen, and Marcel A. J. van Gerven.  Population codes of prior knowledge learned through environmental regularities.  _Scientific Reports_ , 11(1):640, 2021.  doi:10.1038/s41598-020-79366-z.  URL https://doi.org/10.1038/s41598-020-79366-z.
  * Murray et al. [2014] John D Murray, Alberto Bernacchia, David J Freedman, Ranulfo Romo, Jonathan D Wallis, Xinying Cai, Camillo Padoa-Schioppa, Tatiana Pasternak, Hyojung Seo, Daeyeol Lee, and Xiao-Jing Wang.  A hierarchy of intrinsic timescales across primate cortex.  _Nature Neuroscience_ , 17(12):1661–1663, 2014\.  doi:10.1038/nn.3862.  URL https://doi.org/10.1038/nn.3862.
  * Cavanagh et al. [2020] Sean E. Cavanagh, Laurence T. Hunt, and Steven W. Kennerley.  A diversity of intrinsic timescales underlie neural computations.  _Frontiers in Neural Circuits_ , 14, 2020.  ISSN 1662-5110.  doi:10.3389/fncir.2020.615626.  URL https://www.frontiersin.org/articles/10.3389/fncir.2020.615626.
  * Golesorkhi et al. [2021] Mehrshad Golesorkhi, Javier Gomez-Pilar, Federico Zilio, Nareg Berberian, Annemarie Wolff, Mustapha C. E. Yagoub, and Georg Northoff.  The brain and its time: intrinsic neural timescales are key for input processing.  _Communications Biology_ , 4(1):970, 2021.  doi:10.1038/s42003-021-02483-6.  URL https://doi.org/10.1038/s42003-021-02483-6.
  * Amunts et al. [2022] Katrin Amunts, Javier DeFelipe, Cyriel Pennartz, Alain Destexhe, Michele Migliore, Philippe Ryvlin, Steve Furber, Alois Knoll, Lise Bitsch, Jan G. Bjaalie, Yannis Ioannidis, Thomas Lippert, Maria V. Sanchez-Vives, Rainer Goebel, and Viktor Jirsa.  Linking brain structure, activity, and cognitive function through computation.  _eNeuro_ , 9(2), 2022.  doi:10.1523/ENEURO.0316-21.2022.  URL https://www.eneuro.org/content/9/2/ENEURO.0316-21.2022.
  * Mastrogiuseppe and Ostojic [2018] Francesca Mastrogiuseppe and Srdjan Ostojic.  Linking connectivity, dynamics, and computations in low-rank recurrent neural networks.  _Neuron_ , 99(3):609–623.e29, 2018.  ISSN 0896-6273.  doi:https://doi.org/10.1016/j.neuron.2018.07.003.  URL https://www.sciencedirect.com/science/article/pii/S0896627318305439.
  * Vyas et al. [2020] Saurabh Vyas, Matthew D. Golub, David Sussillo, and Krishna V. Shenoy.  Computation through neural population dynamics.  _Annual Review of Neuroscience_ , 43(1):249–275, 2020.  doi:10.1146/annurev-neuro-092619-094115.  URL https://doi.org/10.1146/annurev-neuro-092619-094115.  PMID: 32640928.
  * Beiran et al. [2021] Manuel Beiran, Alexis Dubreuil, Adrian Valente, Francesca Mastrogiuseppe, and Srdjan Ostojic.  Shaping Dynamics With Multiple Populations in Low-Rank Recurrent Networks.  _Neural Computation_ , 33(6):1572–1615, 05 2021\.  ISSN 0899-7667.  doi:10.1162/neco_a_01381.  URL https://doi.org/10.1162/neco_a_01381.
  * Papo [2013] David Papo.  Time scales in cognitive neuroscience.  _Frontiers in Physiology_ , 4, 2013.  ISSN 1664-042X.  doi:10.3389/fphys.2013.00086.  URL https://www.frontiersin.org/articles/10.3389/fphys.2013.00086.
  * Barak [2017] Omri Barak.  Recurrent neural networks as versatile tools of neuroscience research.  _Curr Opin Neurobiol._ , 46:1–6, 2017.
  * Nair and Hinton [2010] Vinod Nair and Geoffrey E. Hinton.  Rectified linear units improve restricted boltzmann machines.  In Johannes Fürnkranz and Thorsten Joachims, editors, _ICML_ , pages 807–814. Omnipress, 2010.  URL http://dblp.uni-trier.de/db/conf/icml/icml2010.html#NairH10.
  * Ma et al. [2006] Wei Ji Ma, Jeffrey M. Beck, Peter E. Latham, and Alexandre Pouget.  Bayesian inference with probabilistic population codes.  _Nature Neuroscience_ , 9(11):1432–1438, Nov 2006\.  ISSN 1546-1726.  doi:10.1038/nn1790.  URL https://doi.org/10.1038/nn1790.
  * Orhan and Ma [2017] A. Emin Orhan and Wei Ji Ma.  Efficient probabilistic inference in generic neural networks trained with non-probabilistic feedback.  _Nature Communications_ , 8(1):138, Jul 2017.  ISSN 2041-1723.  doi:10.1038/s41467-017-00181-8.  URL https://doi.org/10.1038/s41467-017-00181-8.
  * Ichikawa and Kataoka [2022] Kohei Ichikawa and Asaki Kataoka.  Dynamical Mechanism of Sampling-Based Probabilistic Inference Under Probabilistic Population Codes.  _Neural Computation_ , 34(3):804–827, 02 2022\.  ISSN 0899-7667.  doi:10.1162/neco_a_01477.  URL https://doi.org/10.1162/neco_a_01477.
  * Tolhurst et al. [1983] D.J. Tolhurst, J.A. Movshon, and A.F. Dean.  The statistical reliability of signals in single neurons in cat and monkey visual cortex.  _Vision Research_ , 23(8):775–785, 1983.  ISSN 0042-6989.  doi:https://doi.org/10.1016/0042-6989(83)90200-6.  URL https://www.sciencedirect.com/science/article/pii/0042698983902006.
  * Swindale [1998] N. V. Swindale.  Orientation tuning curves: empirical description and estimation of parameters.  _Biological Cybernetics_ , 78(1):45–56, 1998\.  doi:10.1007/s004220050411.  URL https://doi.org/10.1007/s004220050411.
  * Rumelhart et al. [1986] D. E. Rumelhart, G. E. Hinton, and R. J. Williams.  _Learning Internal Representations by Error Propagation_ , page 318–362.  MIT Press, Cambridge, MA, USA, 1986.  ISBN 026268053X.
  * Werbos [1990] P. J. Werbos.  Backpropagation through time: what it does and how to do it.  _Proceedings of the IEEE_ , 78(10):1550–1560, 1990.
  * Kingma and Ba [2014] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2014.  URL http://arxiv.org/abs/1412.6980.
  * Bishop [2006] Christopher M. Bishop.  _Pattern Recognition and Machine Learning_.  Springer, 2006.
  * Mante and et al. [2013] V Mante and et al.  Context-dependent computation by recurrent dynamics in prefrontal cortex.  _Nature_ , 503:78–84, 2013.  doi:10.1038/nature12742.
  * Ichikawa and Kaneko [2021] Kohei Ichikawa and Kunihiko Kaneko.  Short-term memory by transient oscillatory dynamics in recurrent neural networks.  _Phys. Rev. Research_ , 3:033193, Aug 2021.  doi:10.1103/PhysRevResearch.3.033193.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.3.033193.
  * Schrauwen et al. [2007] Benjamin Schrauwen, David Verstraeten, and Jan Campenhout.  An overview of reservoir computing: Theory, applications and implementations.  pages 471–482, 01 2007.
  * Perez-Nieves et al. [2021] Nicolas Perez-Nieves, Vincent C. H. Leung, Pier Luigi Dragotti, and Dan F. M. Goodman.  Neural heterogeneity promotes robust learning.  _Nature Communications_ , 12(1):5791, 2021.  doi:10.1038/s41467-021-26022-3.  URL https://doi.org/10.1038/s41467-021-26022-3.
  * Yang et al. [2019] Guangyu Robert Yang, Madhura R. Joglekar, H. Francis Song, William T. Newsome, and Xiao-Jing Wang.  Task representations in neural networks trained to perform many cognitive tasks.  _Nature Neuroscience_ , 22(2):297–306, 2019.  doi:10.1038/s41593-018-0310-2.  URL https://doi.org/10.1038/s41593-018-0310-2.
  * Yamashita and Tani [2008] Yuichi Yamashita and Jun Tani.  Emergence of functional hierarchy in a multiple timescale neural network model: A humanoid robot experiment.  _PLOS Computational Biology_ , 4(11):1–18, 11 2008.  doi:10.1371/journal.pcbi.1000220.  URL https://doi.org/10.1371/journal.pcbi.1000220.
  * Kurikawa and Kaneko [2021] Tomoki Kurikawa and Kunihiko Kaneko.  Multiple-timescale neural networks: Generation of history-dependent sequences and inference through autonomous bifurcations.  _Frontiers in Computational Neuroscience_ , 15, 2021.  ISSN 1662-5188.  doi:10.3389/fncom.2021.743537.  URL https://www.frontiersin.org/articles/10.3389/fncom.2021.743537.
  * Tanaka et al. [2022] Gouhei Tanaka, Tadayoshi Matsumori, Hiroaki Yoshida, and Kazuyuki Aihara.  Reservoir computing with diverse timescales for prediction of multiscale dynamics.  _Phys. Rev. Research_ , 4:L032014, Jul 2022.  doi:10.1103/PhysRevResearch.4.L032014.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.4.L032014.
  * Yamaguti and Tsuda [2021] Yutaka Yamaguti and Ichiro Tsuda.  Functional differentiations in evolutionary reservoir computing networks.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 31(1):013137, 2021.  doi:10.1063/5.0019116.  URL https://doi.org/10.1063/5.0019116.
  * Bengio et al. [2015] Yoshua Bengio, D. Lee, J. Bornschein, and Zhouhan Lin.  Towards biologically plausible deep learning.  _ArXiv_ , abs/1502.04156, 2015.
  * Lillicrap et al. [2016] Timothy P. Lillicrap, Daniel Cownden, Douglas B. Tweed, and Colin J. Akerman.  Random synaptic feedback weights support error backpropagation for deep learning.  _Nature Communications_ , 7(1):13276, 2016.  doi:10.1038/ncomms13276.  URL https://doi.org/10.1038/ncomms13276.
  * Richards et al. [2019] Blake A. Richards, Timothy P. Lillicrap, Philippe Beaudoin, Yoshua Bengio, Rafal Bogacz, Amelia Christensen, Claudia Clopath, Rui Ponte Costa, Archy de Berker, Surya Ganguli, Colleen J. Gillon, Danijar Hafner, Adam Kepecs, Nikolaus Kriegeskorte, Peter Latham, Grace W. Lindsay, Kenneth D. Miller, Richard Naud, Christopher C. Pack, Panayiota Poirazi, Pieter Roelfsema, João Sacramento, Andrew Saxe, Benjamin Scellier, Anna C. Schapiro, Walter Senn, Greg Wayne, Daniel Yamins, Friedemann Zenke, Joel Zylberberg, Denis Therien, and Konrad P. Kording.  A deep learning framework for neuroscience.  _Nature Neuroscience_ , 22(11):1761–1770, 2019\.  doi:10.1038/s41593-019-0520-2.  URL https://doi.org/10.1038/s41593-019-0520-2.
  * Yang and Wang [2020] Guangyu Robert Yang and Xiao-Jing Wang.  Artificial neural networks for neuroscientists: A primer.  _Neuron_ , 107(6):1048–1070, 2020.  ISSN 0896-6273.  doi:https://doi.org/10.1016/j.neuron.2020.09.005.  URL https://www.sciencedirect.com/science/article/pii/S0896627320307054.
  * Barak et al. [2013] Omri Barak, David Sussillo, Ranulfo Romo, Misha Tsodyks, and L.F. Abbott.  From fixed points to chaos: Three models of delayed discrimination.  _Progress in Neurobiology_ , 103:214–222, 2013.  ISSN 0301-0082.  doi:https://doi.org/10.1016/j.pneurobio.2013.02.002.  URL https://www.sciencedirect.com/science/article/pii/S0301008213000129.  Conversion of Sensory Signals into Perceptions, Memories and Decisions.
  * Cueva and Wei [2018] Christopher J. Cueva and Xue-Xin Wei.  Emergence of grid-like representations by training recurrent neural networks to perform spatial localization.  In _International Conference on Learning Representations_ , 2018.  URL https://openreview.net/forum?id=B17JTOe0-.
  * Yamins and DiCarlo [2016] Daniel L K Yamins and James J DiCarlo.  Using goal-driven deep learning models to understand sensory cortex.  _Nature Neuroscience_ , 19(3):356–365, 2016.  doi:10.1038/nn.4244.  URL https://doi.org/10.1038/nn.4244.
  * Haesemeyer et al. [2019] Martin Haesemeyer, Alexander F. Schier, and Florian Engert.  Convergent temperature representations in artificial and biological neural networks.  _Neuron_ , 103(6):1123–1134.e6, 2019.  ISSN 0896-6273.  doi:https://doi.org/10.1016/j.neuron.2019.07.003.  URL https://www.sciencedirect.com/science/article/pii/S0896627319306014.
  * Dubreuil et al. [2022] Alexis Dubreuil, Adrian Valente, Manuel Beiran, Francesca Mastrogiuseppe, and Srdjan Ostojic.  The role of population structure in computations through neural dynamics.  _Nature Neuroscience_ , 25(6):783–794, 2022.  doi:10.1038/s41593-022-01088-4.  URL https://doi.org/10.1038/s41593-022-01088-4.
  * Lorenz et al. [2011] Dirk M. Lorenz, Alice Jeng, and Michael W. Deem.  The emergence of modularity in biological systems.  _Physics of Life Reviews_ , 8(2):129–160, 2011\.  ISSN 1571-0645.  doi:https://doi.org/10.1016/j.plrev.2011.02.003.  URL https://www.sciencedirect.com/science/article/pii/S1571064511000170.
  * Kashtan and Alon [2005] Nadav Kashtan and Uri Alon.  Spontaneous evolution of modularity and network motifs.  _Proceedings of the National Academy of Sciences_ , 102(39):13773–13778, 2005.  doi:10.1073/pnas.0503610102.  URL https://www.pnas.org/doi/abs/10.1073/pnas.0503610102.
  * Penrose [1955] R. Penrose.  A generalized inverse for matrices.  _Mathematical Proceedings of the Cambridge Philosophical Society_ , 51(3):406–413, 1955.  doi:10.1017/S0305004100030401.