Learning the Synaptic and Intrinsic Membrane Dynamics Underlying Working Memory in Spiking Neural Network Models

2.1 . Units with Diverse Parameter Values Emerge after Training

We applied our method to train spiking networks to perform the context-dependent input integration task previously employed by Mante et al. (2013). Briefly, Mante et al. trained rhesus monkeys to flexibly integrate sensory inputs (color and motion of randomly moving dots presented on a screen). A contextual cue was given to instruct the monkeys which sensory modality (color or motion) they should attend to. The monkeys were required to employ flexible computations as the same modality could be either relevant or irrelevant depending on the contextual cue. Several previous modeling studies have successfully implemented a simplified version of the task and reproduced the neural dynamics present in the experimental data with both continuous-rate RNNs and spiking RNNs converted from rate RNNs (Song et al., 2016a; Miconi, 2017; Kim et al., 2019). With our method, we were able to directly train the first, to our knowledge, spiking RNNs with heterogeneous units whose parameters were within biologically plausible limits.

In order to train spiking RNNs to perform the input integration task, we employed a task paradigm similar to the one used by previous computational studies (Mante et al., 2013; Song et al., 2016a; Miconi, 2017; Kim et al., 2019). A recurrently connected network received two streams of noisy input signals along with a constant-valued signal that encoded the contextual cue (see Figure 1A). The input signals were sampled from a standard gaussian distribution (i.e., with zero mean and unit variance) and then shifted by a positive or negative “offset” value to simulate the evidence presented in the input modalities. The network was trained to produce an output signal approaching either $+$1 or $-$1 depending on the cue and the evidence present in the input signal: if the cued input had a positive mean, the output signal approached $+$1, and vice versa (see Figure 1B, top). The input signal, 150 ms in duration, was given after a fixation period (300 ms), and the network was trained to produce an output signal immediately after the offset of the input signal.

We trained 20 spiking RNNs to perform the context-based input integration task. All of the trainable parameters were initialized with random numbers drawn from a standard gaussian distribution and rescaled to the biologically plausible ranges (see section 4 and Table 1). Each network was trained until the training termination criteria were satisfied (see section 4). On average, $508\pm 46$ training trials were needed for a network to meet the training termination conditions. After training, a wide distribution of the parameters emerged for both excitatory and inhibitory populations (see Figure 1C, top).

Consistent with the previous experimental recordings from cortical neurons, the inhibitory units in our trained RNNs fired at a higher rate compared to the excitatory units (Peyrache et al., 2012). The higher average firing rates of the inhibitory units were largely due to the intrinsic properties that resulted from training. Compared to the excitatory population, the inhibitory units in the trained RNNs had significantly larger input resistance, smaller membrane time constants, and more depolarized resting potential (see Figure 1C; $P<0.0001$, two-sided Wilcoxon rank-sum test). The action potential thresholds and the reset potentials were significantly more depolarized for the inhibitory group. Furthermore, the time constants of the inhibitory synaptic current variable were significantly larger than the excitatory synaptic decay time constants (see Figure 1C).

2.2 . Working Memory Requires Distinct Parameter Distributions

The context-dependent input integration task considered in the previous section did not require complex cognitive skills such as working memory (WM) computations. In order to explore what parameter values are essential for WM tasks, we modified the paradigm to incorporate a WM component by adding a delay period after the delivery of the input signals. The RNN model was trained to integrate the noisy input signals, sustain the integrated information throughout the 300 ms delay period, and produce an output signal (see Figure 1B, bottom). We again trained 20 models for the modified integration task with the same training termination criteria (see section 4). This task required more training trials (on average $1618\pm 346$), but all the models were successfully trained within 2000 training trials.

Overall, the distributions of the trained parameters were similar to those observed from the RNNs trained on the non-WM version of the task (see Figure 1D). The parameters that were significantly different between the two RNN models were the membrane time constant and the synaptic decay time constant. The inhibitory units from the WM model displayed much faster membrane dynamics and slower synaptic decay compared to the inhibitory population of the non-WM model ($P<0.0001$, two-sided Wilcoxon rank-sum test).

To ensure that the patterns of the trained parameters and the distinct distributions of the two parameters (${\tau }_m$ and $\tau$) observed from the delayed integration model were indeed associated with WM computations, we trained RNNs on two additional WM-related tasks: delayed matched-to-sample (DMS) and delayed discrimination (DIS) tasks. For each task, we again trained 20 RNNs. Both task paradigms included two sequential stimuli separated by a brief delay period. For the DMS task, the two input stimuli were either $+$1 or $-$1; if the two sequential had the same sign (i.e., $+$1/$+$1 or $-$1/$-$1), the network was trained to have an output signal approaching $+$1, while if the two stimuli had different signs (i.e., $+$1/$-$1 or $-$1/$+$1), the output signal approached $-$1 (see Figure 2A; see section 4). The two input stimuli for the DIS task were sinusoidal waves with different frequencies, modeled after the task used by Romo, Brody, Hernández, and Lemus (1999) where monkeys were trained to discriminate two vibratory stimuli. If the first stimulus had a higher (lower) frequency, our RNN model was trained to produce a positive (negative) output signal (see Figure 2B; see section 4).

It took longer to train our model on these two tasks compared to the delayed integration task ($7104\pm 3739$ trials for the DMS task and $6985\pm 2112$ trials for the DIS task). The distributions of the tuned parameters from the two WM tasks were similar to the distributions obtained from the delayed integration task (see Figures 2C and 2D). More important, we again observed significantly faster membrane voltage dynamics and slower synaptic decay from the inhibitory units in the DMS and DIS models compared to the inhibitory units from the non-WM task. These findings strongly suggest that the two parameters (${\tau }_m$ and $\tau$) of the inhibitory group contribute to important dynamics associated with WM.

To ensure that the random initialization did not contribute to the heterogeneous distributions after training, we trained another group of RNNs with fixed initial values for the trainable parameters. The trained parameters from the RNNs trained with fixed initial parameter values still displayed heterogeneous values for each trained parameter, indicating that our initialization did not significantly contribute to the heterogeneity (see Figure 3).

2.3 . Shared Intrinsic Properties across Different Working Memory Tasks

Prefrontal cortex and other higher-order cortical areas have been shown to integrate information in a flexible manner and switch between tasks seamlessly (Mante et al., 2013). Along this line of thought, we hypothesized that the intrinsic properties optimized for one WM task should be generalizable to other tasks that also require WM. In order to test this hypothesis, we retrained all the RNNs that were trained in the previous sections to perform the DMS task without tuning the intrinsic parameters. For example, given a network trained on the non-WM integration task, we froze its intrinsic ($R$, ${\tau }_m$, $v_{\text{rest}}$, $v_{\text{reset}}$$\vartheta$) along with the synaptic decay time constant ($\tau$) and optimized the recurrent connections ($W$) only using BPTT (see section 4). Therefore, each of the 20 RNNs trained for each of the four tasks (non-WM integration, delayed integration, DMS, and DIS tasks) was retrained to perform the DMS task. As expected, the average number of trials required to successfully retrain the RNNs previously trained for the DMS task was low at 4409 $\pm$ 3596 (see Figure 4A). The number of trials required to retrain the RNNs from the DIS task was also low at 4180 $\pm$ 2693. The RNNs trained for the delayed integration task took longer to retrain, at 5392 $\pm$ 2198. The non-WM RNNs required the most number of training trials to perform the DMS task (9648 $\pm$ 2933). These findings indicate that the intrinsic properties from one WM model are transferable to other WM models.

Based on these previous results, the membrane time constant (${\tau }_m$) and the synaptic decay ($\tau$) variables appeared to be the two most important parameters for the transferability of WM. To test this, we repeated the retraining procedure with both ${\tau }_m$ and $\tau$ either fixed (“frozen”) or optimized (“tuned”) for the non-WM RNNs (see section 4). For the frozen condition (i.e., ${\tau }_m$ and $\tau$ frozen while the other parameters optimized), the number of trials required to retrain the non-WM RNNs to perform the DMS task was high and not significantly different from the number of trials it took with the intrinsic parameters fixed (see Figure 4B). On the other hand, retuning only ${\tau }_m$ and $\tau$ with the other parameters fixed (i.e., tuned condition) resulted in a significant reduction in training time (see Figure 4B), suggesting that these two parameters are indeed critical for performing WM. Optimizing both ${\tau }_m$ and $\tau$ resulted in a significant decrease in ${\tau }_m$ for both excitatory and inhibitory populations (see Figure 4C). The synaptic decay values decreased for the excitatory units after retuning (see Figure 4D left). For the inhibitory population, $\tau$ was significantly increased (see Figure 4D right).

2.4 . Membrane and Synaptic Decay Time Constants Critical for WM Maintenance

Pyramidal excitatory neurons and parvalbumin (PV) interneurons make up the majority of the neuronal cell population in the cortex, and they have been shown to specialize in fast and reliable encoding of information with high temporal precision (Tremblay, Lee, & Rudy, 2016). To further investigate if the fast membrane and slow synaptic dynamics of the units from our WM RNNs are aligned with previous experimental findings and to probe how they contribute to WM maintenance, we manipulated ${\tau }_m$ and $\tau$ during different epochs of the DMS task paradigm.

For each of the RNNs trained from the DMS task, we first divided the population into two subgroups based on their ${\tau }_m$ values (see section 4). The short ${\tau }_m$ group contained units whose ${\tau }_m$ was smaller than the lower quartile value, while the long ${\tau }_m$ group contained units whose ${\tau }_m$ was greater than the upper quartile. During each of the four epochs (fixation, first stimulus, delay, and second stimulus), we then inhibited the two ${\tau }_m$ subgroups separately by hyperpolarizing them and assessed the task performance (see section 4). As shown in Figure 5, inhibiting the short ${\tau }_m$ subgroup during the two stimulus windows significantly impaired task performance (see Figures 5B and 5D), while disrupting the long ${\tau }_m$ group did not result in significant changes in task performance in all four task epochs.

We repeated the above analysis with two subgroups derived from a quartile split of the synaptic decay time constant ($\tau$; see section 4). Suppressing the synaptic connections in the long $\tau$ subgroup during the first stimulus window and the delay period significantly impaired task performance (see Figures 5B and 5C). Inhibiting the short $\tau$ group at any of the four epochs did not affect the task performance.

Therefore, the units with the fast membrane voltage dynamics (${\tau }_m$) were important for encoding of stimuli, while the slow synaptic dynamics ($\tau$) were critical for maintaining the first stimulus information throughout the period spanning from the first stimulus window to the end of the delay window.

4.1 . Spiking Network Structure and Discretization

Our spiking RNN model consisted of $N$ integrate-and-fire (LIF) units is governed by

where ${\tau }_{m,i}$ is the membrane time constant of unit $i$, $v_i\left(t\right)$ is the membrane voltage of unit $i$ at time $t$, $v_{\text{rest},i}$ is the resting potential of unit $i$, and $R_i$ is the input resistance of unit $i$, and $\xi$ is the membrane voltage spontaneous fluctuation. $I_i\left(t\right)$ represents the current input to unit $i$ at time $t$, which is given by

where $N$ is the total number of units in the network, $s_{ij}\left(t\right)$ is the filtered spike train of unit $j$ to unit $i$ at time $t$, and $I_{{\text{ext}}_i}\left(t\right)$ is the external current source into unit $i$ at time $t$. For this study, $N=400$ for all tasks and networks trained.

The external current ${\mathbf{I}}_{\text{ext}}\left(t\right)$ encodes the task-specific input at time $t$,

where the time-varying stimulus signals $bmu\left(t\right)\in {\mathbb{R}}^{N_{\text{in}}\times 1}$ are fed into the network via $W_{\text{in}}\in {\mathbb{R}}^{N\times N_{\text{in}}}$, which can be viewed as presynaptic connections to the network that convert analog input into firing rates. $N_{\text{in}}$ corresponds to the number of channels in the input signal.

We used a single exponential synaptic filter

where ${\tau }_{ij}$ is the synaptic decay time constant from unit $j$ to unit $i$, $w_{ij}$ is the synaptic strength from unit $j$ to unit $i$, $t_j^{\left(k\right)}$ denotes the time of the $k$th action potential of unit $j$, and $\delta \left(x\right)$ is the Dirac delta function. Once the membrane voltage of the unit $i$ crosses its action potential threshold (${\vartheta }_i$), its membrane voltage is brought back down to its reset voltage ($v_{\text{reset},i}$).

The output of our spiking model at time $t$ is given by

where $W_{\text{out}}\in {\mathbb{R}}^{1\times N}$ are the readout weights, and $\mathbf{r}\left(t\right)\in {\mathbb{R}}^{N\times 1}$, which can be interpreted as the firing rate of units, are given by

where ${\tau }_{r,i}$ is the synaptic decay time constant of firing rate estimate for unit $i$.

We converted the continuous-time differential equations to discrete-time iterative equations and used numerical integration (Euler's method) to solve the equations. The membrane voltage $\mathbf{v}\in {\mathbb{R}}^{1\times N}$ at step $n+1$ is given by

where $\Delta t$ is the sampling rate (or step size), which was set $\Delta t=1$ ms for this study, ${\mathbf{\tau }}_m\in {\mathbb{R}}^{1\times N}$ is the membrane time constant, ${\mathbf{v}}_{\text{rest}}\in {\mathbb{R}}^{1\times N}$ is the resting potential, $\odot$ refers to Hadamard operation (element-wise multiplication), $\frac{·}{·}$ refers to the element-wise division, and $\mathbf{R}\in {\mathbb{R}}^{1\times N}$ is the input resistance. The term $c\mathcal{N}(0,\Delta t)$ injects spontaneous membrane fluctuations, where $\mathcal{N}(0,\Delta t)\in {\mathbb{R}}^{1\times N}$ is a gaussian random vector consisting of $N$ independent gaussian random variables with mean 0 and variance $\Delta t$, and $c$ is the scaling constant for the amplitude of fluctuations, set as $c=5$ throughout the study.

There are two time-varying terms in equation 2.7, the membrane voltage after reset (${\tilde{\mathbf{v}}}^{\left(n\right)}$) and input current (${\mathbf{I}}^{(n+1)}$). The voltage reset in the LIF model after action potentials at step $n$ is formulated as

where ${\mathbf{v}}_{\text{reset}}\in {\mathbb{R}}^{1\times N}$ is the reset potential, $\mathbf{\vartheta }\in {\mathbb{R}}^{1\times N}$ is the action potential thresholds, and $H\left(x\right)$ is the element-wise Heaviside step function. The term $H\left({\mathbf{v}}^{(n+1)}-\mathbf{\vartheta }\right)$ represents the spiking output activities at step $n+1$. The input current at step $n+1$ is given by

where $\mathbf{1}\in {\mathbb{R}}^{1\times N}$ is the column vector with all ones and $S^{\left(n\right)}$ is the filtered spike train matrix at step $n$, which follows the iteration

where $T\in {\mathbb{R}}^{N\times N}$ is the matrix of synaptic decay time constants and $W\in {\mathbb{R}}^{N\times N}$ is the matrix of synaptic strengths. Here, $W\in {\mathbb{R}}^{N\times N}$ is a matrix, and $H\left({\mathbf{v}}^{\left(n\right)}-\mathbf{\vartheta }\right)\in {\mathbb{R}}^{1\times N}$ is a row vector. The notation $A\odot \mathbf{v}$ refers to element-wise multiplication of matrix $A$ row by row with the row vector $\mathbf{v}$.

The output at step $n+1$ is computed by

in which

where ${\mathbf{\tau }}_r\in {\mathbb{R}}^{1\times N}$ is the synaptic decay time constants of firing rate estimate.

4.2 . Training Details

In this study, we used only the supervised backpropagation of errors learning algorithm. The loss function ($\mathcal{L}$) is defined in terms of the root mean square error (RMSE) with respect to a task-specific target signal ($\mathbf{z}$) and the network output signal ($\mathbf{o}$),

where $M$ is the total time steps.

We used adaptive moment estimation (ADAM) stochastic gradient descent algorithm (Kingma & Ba, 2014) with mini-batch training. The mollifier gradient approximations were employed to address non-differentiability problem associated with the spiking process (see section 4.3). The learning rate was set to 0.01, the batch size was set to 10, and the first and second moment decay rates were 0.9 and 0.999, respectively. The trainable parameters include input weights ($W_{\text{in}})$, synaptic strengths ($W$), readout weights ($W_{\text{out}}$), synaptic decay time constants ($T$), membrane time constants (${\mathbf{\tau }}_m$), input resistances ($\mathbf{R}$), resting potentials (${\mathbf{v}}_{\text{rest}}$), reset voltages (${\mathbf{v}}_{\text{reset}}$), action potential thresholds ($\mathbf{\vartheta }$), and synaptic decay time constants for firing rate estimates (${\mathbf{\tau }}_r$).

A nonlinear projected gradient method was used to constrain parameters within the biologically realistic ranges described in Table 1. A linear projection map forces some solutions to be projected on the boundary. That is, there are always some units whose parameters take the min and max values of the constraint. On the other hand, a nonlinear projection guarantees that no values are on the boundary almost surely, a more realistic situation to consider. Specifically, to bound a parameter $p$ at iteration $i+1$ into the range $[p_{min},p_{max}]$, we have

where ${\tilde{p}}_{i+1}$ is the projected solution of parameter $p$ at iteration $i+1$, $p_{i+1}$ is the unconstrained solution given by the gradient descent algorithm at iteration $i+1$, $p_{\text{max}}$ and $p_{\text{min}}$ are the maximum and minimum values of parameter $p$, and $\sigma \left(x\right)$ is the sigmoid function, defined as

We initialized all parameters, except the input weights ($W_{\text{in}}$), as samples from the standard gaussian distribution with zero mean and unit variance, whereas the input weights were drawn from gaussian distribution with zero mean and variance 400. This is because our input signals were bounded within the range $[-1,1]$, insufficient to bring the membrane voltage from the resting potential above the action potential threshold. Hence, to accelerate training, it was necessary to make sure units were excited by the input signals in the first place. The synaptic strength matrix ($W$) was also initialized sparse, with the percentage of connectivity being only 20%. We say the network successfully did the task if the output signal hits above $+0.8$ (or below $-0.8$) if the target output is $+$1 (or $-$1). We stopped training when the loss ($\mathcal{L}$) is less than 15 and the accuracy over 100 trials is above 95%.

The method proposed by Song et al. (2016a) was used to impose Dale's principle with separate excitatory and inhibitory populations. The synaptic connectivity matrix ($W$) in the model was parameterized by

where ${\tilde{W}}_{i+1}$ is the resulting matrix that encoded separate populations at update step $i+1$, $W_{i+1}$ is the solution given by the gradient descent algorithm at step $i+1$, and ${[·]}_{+}$ is the rectified linear unit (ReLU) operation applied at the end of each update step. The ReLU operation is to ensure that entries of the matrix are always nonnegative before being multiplied by the matrix $D$, as the negative weight connections updated from gradient descent are pruned by the end of each update. The diagonal matrix ($D\in {\mathbb{R}}^{N\times N}$) encodes $+$1 for excitatory units and $-$1 for inhibitory units. The value of matrix ($D$) was randomly assigned before training according to a preset proportion between inhibitory and excitatory units, and the value $D$ was fixed through the whole training process. The I/E units proportion in this study was 20% to 80%.

In order to capture the biologically realistic dynamics of SNNs, the temporal resolution ($\Delta t$) was set to be no longer than the duration of absolute refractory period to ensure that the spiking activities are not affected by the numerical integration process. Therefore, we set $\Delta t=1$ ms during training. Due to the vanishing gradient problem occurring in training RNNs (Hochreiter, Bengio, Frasconi, & Schmidhuber, 2001), with $\Delta =1$ ms, it is impossible to train tasks with duration longer than 1 second (i.e., $M>1000$). It is notable that in the above formulation, only membrane time constant (${\tau }_m$) and synaptic time decay ($\tau$) are dependent on the sampling rate ($\Delta t$; see equations 4.7 and 4.10). Hence, after the models are trained, we can make sampling rate ($\Delta t$) smaller (i.e., having finer temporal resolution) while still keeping the same dynamics of the trained networks. Increasing $\Delta t$ by a factor is equivalent to decreasing $\tau$ and ${\tau }_m$ altogether by the same factor, as $\tau$ and ${\tau }_m$ are inversely proportional to $\Delta t$ in equations 4.7 and 4.10. Hence, to train a network performing tasks with duration longer than 1 second, we need to make the temporal resolution coarser (i.e., increasing $\Delta t$ by a factor $s$) so that with the same trainable range of time steps (i.e., a fixed $M\le 1000$), the duration of task becomes longer by the same factor $s$. This “decrease in temporal resolution” can be interpreted as shortening $\tau$ and ${\tau }_m$ instead of an actual decrease in temporal resolution. Applying this trick enables us to train tasks with arbitrary duration by rescaling the ranges of $\tau$ and ${\tau }_m$ into a smaller one while still making the spiking activities biologically realistic. In practice, we simply scaled down $\tau$ and ${\tau }_m$ by a factor $s=3$ with a fixed number of time steps ($M$), and later during the testing stage, we rescaled $M,\tau$ and ${\tau }_m$ up by the same factor $s$.

4.3 . Mollifier Gradient Approximations

In the above formulation, the Heaviside step function $H\left(x\right)$ is not continuous. As a result, the loss function $\mathcal{L}$ is not differentiable. This poses a major problem when applying the traditional backpropagation algorithm for training neural networks because the backpropagation algorithm uses gradient descent methods that require the function being minimized to be differentiable, or at least to be continuous. However, the derivative of Heaviside step function $H\left(x\right)$ is Dirac delta function $\delta \left(x\right)$, which is 0 everywhere except at 0, where the function value is $\infty$. It is difficult to use this derivative for the gradient descent methods because the value of the gradients is 0 almost everywhere.

To address the discontinuity problem, we employed mollifier gradient method proposed by Ermoliev et al. (1995). The method can be applied to any strongly lower semicontinuous functions to find local minima following an iterative gradient descent in which the gradients change over iterations based on averaged functions derived from the original objective function. The family of averaged functions $f_{\varepsilon }$ of function $f$ is defined by convolution of $f$ with a mollifier: ${\psi }_{\varepsilon }$

where ${\psi }_{\varepsilon }\in \{{\psi }_{\varepsilon }:{\mathbb{R}}^n\to {\mathbb{R}}_{+},\varepsilon >0\}$, a family of compactly supported (generalized) functions named mollifiers that satisfy

It was shown that for any strongly lower semicontinuous functions $f$, the averaged functions $f_{\varepsilon }$ epi-converge to $f$ as $\varepsilon \to 0$, a type of convergence that preserves the local minima and minimizers. Therefore, it is possible to use the gradients of averaged functions to minimize the original lower semicontinuous functions and find the local minima. We used the conventional family of mollifiers obtained by normalizing a probability density function $\psi$:

In our case, $n=1$ as the domain of $H\left(x\right)$ is the real line:

For any $\varepsilon >0$, the gradient of $H_{\varepsilon }\left(g\left(x\right)\right)$ with respect to parameter $p$ is given by

where $\psi$ is some symmetric density function and $g\left(x\right)$ is any function with $\mathbb{R}$ as its codomain. Since our goal was not to find a local minimum $x^{*}$ that satisfies the optimality condition $\underset{\varepsilon \to 0}{lim}\parallel \nabla f_{\varepsilon }\left(x^{*}\right)\parallel =0$ as defined by Ermoliev et al. (1995), but rather to minimize the loss function for its value to be sufficiently small so that the network can perform the task correctly, we did not vary the gradients during the minimization process. Instead, we fixed an approximation of the gradient and used the approximation throughout the training process. We chose the normalized box function, that is, the density function of uniform distribution $\mathcal{U}(-\varepsilon /2,\varepsilon /2)$, as the kernel,

and fixed $\varepsilon =5$.

We found no difference in the trained models with different choices of $\varepsilon$ as long as the value was large enough to keep the gradients active so that they did not vanish through time steps. There was also no difference between models trained with fixed $\varepsilon$ and those trained with the original scheme in Ermoliev et al. (1995) where $\varepsilon \to 0$ as the number of iterations increases. The purpose for fixing the value of $\varepsilon$ was to compare the training epochs (iterations) among the retraining paradigms (see Figure 3) with the same gradient.

4.4 . Retraining Models for DMS Task

To test whether intrinsic properties optimized for one WM task are generalizable to other tasks that also require WM, we retrained our models to perform the DMS task with all intrinsic properties fixed. In contrast to the training paradigm described in the previous sections, the trainable parameters for retraining only include input weights ($W_{\text{in}})$, synaptic strengths ($W$), and readout weights ($W_{\text{out}}$). Each of the 20 RNNs trained for each of the four tasks (non-WM integration, delayed integration, DMS, and DIS tasks) used in this study was retrained to perform the DMS task.

To test whether synaptic decay time constants ($\tau$) and membrane time constants (${\tau }_m$) are the most crucial parameters for transferability of WM tasks, we repeated the retraining procedure with both ${\tau }_m$ and $\tau$ either fixed or optimized for the non-WM RNNs. The RNNs optimized to perform the context-based input integration task were used for retraining under two schemes: the tuned scheme and the frozen scheme. For the tuned scheme, the trainable parameters include input weights ($W_{\text{in}})$, synaptic strengths ($W$), readout weights ($W_{\text{out}}$), synaptic decay time constants ($T$), membrane time constants (${\mathbf{\tau }}_m$), and synaptic decay time constants for firing rate estimates (${\mathbf{\tau }}_r$). For the frozen scheme, the trainable parameters include input weights ($W_{\text{in}})$, synaptic strengths ($W$), readout weights ($W_{\text{out}}$), input resistances ($\mathbf{R}$), resting potentials (${\mathbf{v}}_{\text{rest}}$), reset voltages (${\mathbf{v}}_{\text{reset}}$), and action potential thresholds ($\mathbf{\vartheta }$).

4.5 . Units Function Analysis

For Figure 4, we manipulated ${\tau }_m$ and $\tau$ during different epochs of the DMS task paradigm to investigate if fast membrane and slow synaptic dynamics are responsible for WM maintenance. For each of the RNNs trained from the DMS task, we first divided the population into two subgroups based on their ${\tau }_m$ values. The short ${\tau }_m$ group contained units whose ${\tau }_m$ was smaller than the median value of ${\tau }_m$ of all units in the RNN, while the long ${\tau }_m$ group contained units whose ${\tau }_m$ was greater than the median value. The average median value of ${\tau }_m$ across all 20 models was $19.64\pm 2.45$ ms. During each of the four epochs (fixation, first stimulus, delay, and second stimulus), we inhibited the two ${\tau }_m$ subgroups separately by hyperpolarizing them and then assessed the task performance. The hyperpolarization was done by setting the membrane voltage $v=-100$ mV for the intended subgroup of units. Similar to the training stage, we say that the network successfully did the task if the output signal hits above $+0.8$ (or below $-0.8$) if the target output is $+$1 (or $-$1). If the target output is between $-0.8$ and $+0.8$, the network is considered having no response. If the output signal is above $+0.8$ (or below $-0.8$) while the target output is $-$1 (or $+$1), we say that the network gives an incorrect response.

We conducted a similar analysis based on two subgroups of synapses derived from a quartile split of synaptic decay time constant ($\tau$). The short $\tau$ group contained synapses whose $\tau$ was smaller than the 25th percentile of all $\tau$ in the RNN, while the long $\tau$ group contained synapses whose $\tau$ was greater than the 75th percentile. The average 25th percentile across all 20 models was $25.36\pm 2.40$ ms, and the average 75th percentile was $66.18\pm 1.17$ ms. The targeted subgroup of synapses was suppressed by setting the connection strength $w=0$ during each of the four epochs of DMS task.