γ And β Band Oscillation in Working Memory Given Sequential or Concurrent Multiple Items: A Spiking Network Model

Visual Abstract


Introduction
Working memory (WM), the ability to actively maintain and manipulate information in the absence of stimuli, plays a crucial role in cognitive function and executive control of behavior (Baddeley, 2003).The information maintained in WM can enter the brain concurrently or sequentially, meaning that WM cannot only hold multiple pieces of information arriving simultaneously as in visuospatial WM tasks (Bays and Husain, 2008;Zhang and Luck, 2008;Gorgoraptis et al., 2011), but also information presented sequentially, as in visual (Gorgoraptis et al., 2011) or speech processing (Cowan, 2010).The neural implementation of information maintenance is not well understood.The oscillatory model proposes that the information of one item is represented in WM by the reactivation of neurons in the g cycle within nested g -theta oscillations, mediated by a slow after-depolarization (ADP) with a time constant that should match the theta oscillation (Lisman and Idiart, 1995;Jensen and Lisman, 1998).The dynamic model posits that memory is maintained by item-specific patterns of synaptic plasticity and that neurons exhibit a nonstationary and short-lived attractor activity, in which only one memory representation can be active at a time, but successive reactivations of neuronal pools memorize different items (Mongillo et al., 2008;Lundqvist et al., 2016Lundqvist et al., , 2018;;Mi et al., 2017).The persistent activity model asserts that persistent activity of neurons induced by Nmethyl-D-aspartate receptor (NMDAR)-mediated recurrent synaptic current encodes the corresponding items during the delay period (Compte et al., 2000;Edin et al., 2009;Wei et al., 2012).The ability of WM to maintain concurrent and sequential information challenges these models.On the one hand, the oscillatory model and dynamic models can use oscillatory activity to maintain sequential information, but it is difficult for them to represent concurrent information (Constantinidis et al., 2018); on the other hand, the persistent models can manipulate concurrent information but cannot address sequential stimuli and often focus on the persistent activity without oscillation.Experiments have shown that g band oscillations are involved in WM maintenance, in particular, that WM load enhances the g band oscillations but suppresses the alpha band oscillations.Here, we aim to construct a biophysically plausible network model to implement the storage of concurrent and sequential information and to investigate how WM load alters oscillations in the network.
The biological brain consists of microcircuits with three types of neurons that play key roles in generating oscillations: pyramidal cells, low-threshold nonfast-spiking (nFS), and high-threshold neuronal FS (Chen et al., 2017;Veit et al., 2017;Domhof and Tiesinga, 2021;Hahn et al., 2022).The pyramidal cells are excitatory, while the FS and nFS are GABAergic inhibitory.Pyramidal cells synapse on nFS and FS.In turn, nFS neurons synapse back onto the dendrites of the pyramidal cell, while FS neurons synapse back onto the soma of the pyramidal cell, creating a dynamic feedback loop that regulates excitatory and inhibitory activities.FS and nFS also inhibit each other, resulting in a competitive relationship.At rest, the nFS is more active, inhibiting the FS and forming a b oscillation with the pyramidal cells.When a perceptual stimulus is presented, the FS cells are activated, inhibiting the nFS cells and creating a g oscillation with the pyramidal cells (Chen et al., 2017;Hahn et al., 2022).Given the importance of oscillations in higher cognitive functions, we aim to explore the role of oscillations and microcircuitry in working memory.
We proposed a spiking neural network to implement WM with concurrent and sequential directional information (Gorgoraptis et al., 2011).The network consists of two-compartment pyramidal cells, nFS, and FS.These cells and compartments are interconnected as in the biological brain.The cells are uniformly arranged in a ring according to their preferred direction (Compte et al., 2000).When activated, pyramidal cells activate pyramidal cells with similar preferences as well as nearby FS cells and nFS cells on the ring.FS and nFS inhibit neighboring neurons on the ring.We found that, regardless of whether the directions are presented concurrently or sequentially, the network elicits corresponding localized activities that persist throughout the delay period.The persistent activity suggests that the model can successfully maintain multiple stimuli presented simultaneously or sequentially in working memory.Furthermore, before cue presentation, the interaction between pyramidal cells and nFS dominates the activity of the network and causes a low-band oscillation (10-30 Hz); the cue presentation induces strong excitation and recruits FS into the network, and the interaction between pyramidal cells and FS dominates the activity of the network, enhancing the g oscillation (35-100 Hz), which persists throughout the delay period.

Model architecture
Our model has 4096 excitatory pyramidal cells, 512 FS cells, and 512 nFS cells.We have two reasons to choose the number of neurons.First, the ratio of the number of excitatory neurons over that of inhibitory neurons is ;4:1 (Braitenberg and Schutz, 1991).Thus, we chose 4906 pyramidal neurons and 1024 inhibitory interneurons (FS 1 nFS, 1024).Second, a network size of 2 n neurons is desired for the fast Fourier transform (FFT) which recursively divides the input data into smaller subsets, and conducts Fourier transform computations on subsets.These neurons are evenly distributed in a ring configuration and are connected by AMPA, NMDA, and GABA synapses, forming an interconnected spiking network.According to our model, the strength of connections between neurons and the angular differences in their distribution across the ring configuration follows a Gaussian distribution.Neurons that are closer in angle on the ring have stronger connections.The footprint of the connections can be described as Equation 1: (1) Equation 1.The standard deviations of the Gaussian distributions for connectivity footprint in our model were as follows: -From pyramidal cells soma to pyramidal cells soma (s s!s ): 12.76°-From pyramidal cells soma to FS cell (s s!fs ): 7.05°-From pyramidal cells soma to nFS cell (s s!nfs ): 1.41°-From FS cell to pyramidal cells soma (s fs!s ): 8.46°-From FS cell to nFS cell (s fs!nfs ): 1.41°-From nFS cell to pyramidal cells dendrite (s nfs!d ): 5.02°T he value of parameters J 1 s!s and s s!s used in our study were adapted from Compte et al. (2000) and Wei et al. (2012).We further introduced structured excitatory projections from pyramidal cells to FS and nFS and inhibitory projections from FS to pyramidal soma and nFS cells or from nFS to pyramidal dendrite, referencing local circuit properties described by Wang et al. (2004) and Wang and Yang (2018).FS belongs to the large basket cell that projects to broader region, while nFS has narrow dendritic and axonal arbors and serves local inhibition.Generally, excitatory neurons project to a larger radius (around 200mmphysiologically; as shown by Fitzpatrick et al., 1985) compared with inhibitory interneurons (around 100mm; as demonstrated by Spiro et al., 1999) based on biological experiments.Nevertheless, the lateral innervation width of interneurons may be similar to that of excitatory neurons.It is because excitatory neurons project to a larger radius, activating inhibitory neurons within that radius to locally innervate inhibition (as assumed by Compte et al).The synaptic weights were calculated by the production of footprint and the maximum conductance as follows: The G syn parameters for these connections are specified in Table 1 (Compte et al., 2000;Wei et al., 2012).To simplify the calculations, we have assumed that the remaining inhibitory connections in our model (e.g., inhibitory connection from nFS to FS or recurrent connection from FS to FS) follow a uniform distribution (a simple diagram is shown in Fig. 1).

Neuron models
Pyramidal cells and inhibitory cells follow the leaky integrate-and-fire model (Tuckwell, 1988).The neurons generate spikes when their membrane potentials reach threshold (V th ¼ À50mV.)And their membrane potential is reset to V reset ¼ À60mV immediately after a spike for t s ¼ 2ms; t d ¼ 2ms; t fs ¼1ms; t nfs ¼1ms as refractory periods.The Pyramidal cells have two compartments (Wang, 1998): soma(s) and dendrite(d).The membrane potential of the somatic compartment (Vs) and dendrite compartment (Vd) of pyramidal cells obey Equations 3 and 4: (3) V 0 ¼ À70 mV, C m ¼ 0:5 nF, g sL ¼ 0:025nS, g dL ¼ 0:025nS, gc 1 ¼ 0:25nS (Troyer and Miller, 1997;Wei et al., 2012).The parameters p1 and p2 characterize the difference of current intensity resulting from the identical current in soma and dendrite because of the difference in membrane surface area.Specifically, p1 represents the ratio of somatic area to the total neuronal area, and p2 represents the ratio of dendritic area to the total area (Wang, 1998).In the current model, p1 is set to 0.6, and p2 is set to 1-p1, which is 0.4.FS and nFS cells contain single compartments, and their membrane potentials obey Equations 5 and 6.Please note that when the neuron (soma) fires, the membrane potential of the dendrite component is not reset to À60 mV.The dendrite component generates its own spikes and resets the membrane potential after each spike.

Synapse models
Synaptic currents are mediated by AMPA, NMDA, and GABA transmissions.Three types of synaptic currents follow the equations: The conductance "g" in our calculations primarily follows the Equation 2 ] in Equation 10 equals 1 mM (Jahr and Stevens, 1990).The gating variable s for AMPA or GABA (s j;AMPA and s j;GABA in Eqs. 9, 11) follows the dynamic in Equation 12,and the dynamic of s (s j;NMDA ) for NMDA was modelled as Equations 13 and 14: t k is the spike sequence of the neuron.The rising time constant for NMDA (t x;NMDA ) is 2 ms, and the decay time constant t s is 2 ms, 10 ms, and 100 ms for AMPA, GABA, and NMDA, respectively (Wei, 2012).
Synapses from Pyrs to nFS mediated by NMDAR were facilitated in short term and follow the equation: U is set to 0.4 in our simulation (Mongillo et al., 2008).

Simulation protocol
We simulated two types of WM tasks: concurrent-cue tasks (CCTs) and sequential-cue tasks (SCTs), similar to the experiments in (Gorgoraptis et al., 2011).In CCTs, cues are presented to neurons concurrently; each pyramidal cell indexed by u received a current I ext u ð Þ during the concurrent cue presentation period from 500 to 700 ms (Fig. 1B).In concurrent-cue tasks, an array of cues is presented to the network from 500 to 700 ms, followed by a delay period of 1300 ms (upper panel).In sequential-cue tasks, cues were presented one after another with an interval of 100 ms, and each cue lasted for 200 ms.
where u in;k is the kth direction and n is the number of directions in the cue array.I 0 ¼ 0:5 nA; s s ¼ 2 degree.In SCTs, pyramidal cells received multiple currents sequentially.Each current can be formulated as: The kth current is applied at 5001300 1B) during the presentation of the kth cue (Wei et al., 2012).

Noise implementation
Noise is applied to pyramidal cells and inhibitory cells via AMPA synapses: g noiseÀpyr ¼ 0:0039 nS; g noiseÀpv ¼ 0:0019 nS; g noiseÀcb ¼ 0:0022 nS; and s noiseÀi follows: where l is the arrival rate of noise.We set l ¼ 1 kHz, which is equivalent to 1000 synaptic input at 1 Hz (Wei et al., 2012).

Starting parameters
To establish a more realistic and stable initial condition, we conducted a recording procedure to obtain a consistent value of multiple variables for each type of neurons.The following variables were used: the membrane potential ("V m "), the gating variables ("s"  and S3) were sequentially presented to the network.The horizontal axis is the time in milliseconds, and the vertical axis is the preferred direction of the pyramidal cell.s1-s3 represent the indices of the cue, and the dotted lines denote the cue presented interval(s).Please note that these symbols have the same interpretations in the following figures.
and "x") for synapses, the time of the last firing activity ("LastTimeEachFired" minus 250 ms), and the AMPA synaptic current ("I_AMPA") for the pyramidal soma (over a period of 10 ms), and these recorded values were then used as initial conditions for the algorithms.Initially, we performed a model run with constant initialization values of À51 mV for the membrane potential and 1e-34 for the "s" parameter without these recordings.

Numeric integration and code accessibility
The model was programmed in MATLAB code and the gating variables were integrated using a second-order Runge-Kutta (RK2) algorithm.
The membrane potential is integrated with Euler scheme as Equations 3-6.This different choice in the integration method is because changes in gating variables during neuronal firing are more complex and usually nonlinear, whereas changes in the subthreshold membrane potential are approximately linear.We conducted a convergence test on two neurons with AMPA connections and provided a 2-mV constant stimulus to the presynaptic neuron for 2000 ms to examine the difference of spike time of postsynaptic neuron between Euler-RK2 and RK2 methods.We found that the spike time errors generated by the two methods are almost negligible for sufficiently small dt (such as 0.02, as in the experiments).By variation of the timesteps, we found the convergence order of integration method R ¼ log 2 jje new jj Àlog 2 jje old jj log 2 jjDt new jj Àlog 2 jjDt old jj equals 2.95, indicating that the error decays to 0.129 of its original error if the timestep is decreased by half.This approach reduces the computational cost without compromising accuracy.The code/software described in the paper is freely available online at https://github.com/scientific-lab/Gamma-and-Beta-Band-Oscillation-in-Working-Memory.We performed our research using GPU(K80) with Ubuntu 18.04.6LTS, Microsoft Surface Pro 5 with Windows 10 operating system, and MacBook Pro 11,1 with Windows 10 operating system.

Firing rate calculation
We recorded the firing activities of neurons (see Fig. 2) and counted the spikes of a certain type of neuron along the time in a single trial and approximated the firing rate by applying a one-dimension filter with a window size of 4 ms to the spike train.

Local field potential and spectrograms
We approximated the local field potential (LFP) using the total synaptic currents onto pyramidal cells following the method proposed by Mazzoni et al. (2015): where t AMPA ¼ 6ms.We applied a one-dimension filter with a window size of 4 ms to the LFP and calculated the power spectrogram of LFP using continuous wavelet transform with Morlet wavelet.The background energy of our model demonstrates the 1/f characteristic for high frequency oscillation (refer to Miller and Wang, 2006;Wang, 2010).We first computed the power-frequency correlation of our model without external input shown in Figure 3A.We observed that the power for low-frequency oscillation remains almost constant, but decays in accordance with a power-law trend.The power-frequency data were smoothed using an average of 250 points and represented by the red line to signify the background power of oscillation.We subtracted the value of the red line from the simulated power.Figure 3B illustrates the power of LFP, which has been normalized by the red curve in Figure 3A.This normalization is applied in subsequent spectral analysis to obtain the normalized power.

Burst rate calculation
A burst in a particular frequency is defined as an interval when spectral power exceeds twice the standard deviation above the average value of that frequency and lasts for at least three cycles (Lundqvist et al., 2016).After we located the bursts in a single trail, we applied the same process to multiple trails and calculated the frequency of the bursts in the b (10-30 Hz) or g (35-100 Hz) band over time.Similar to how we extracted the firing rate, we filtered the calculated burst frequency with a one-dimensional Gaussian filter, with a window size of 60 ms, to approximate the burst rate for a certain frequency range.

Statistical comparison
In our study, we employed independent sample t tests and ANOVA to compare between groups after conducting tests for homogeneity of variance and making necessary corrections.We also provided the effect size and confidence interval for each significant test result.When necessary, we conducted multiple comparison analyses to offer comprehensive results.All statistical analyses were performed using MATLAB.Our model generated substantial data from multiple simulations with high power and narrow confidence intervals.Based on the results from the online power and sample size calculator (https://www.gigacalculator.com/calculators/),the power for all statistical tests were close to 100%.

Results
The continuous spiking network can hold concurrent and sequential cues We considered a spiking neural network consisting of one population of 4098 pyramidal neurons, one population of 512 FS, and one population of 512 nFS.Pyramidal cells are uniformly distributed on a ring according to their preferred cue angles (Fig. 1A; Compte et al., 2000), mimicking the columnar organization of the monkey PFC (Goldman-Rakic, 1995;Rao et al., 1999;Constantinidis et al., 2001).Using this continuous spiking network model, we investigated how WM manipulates sequential and concurrent cues.We first presented one direction to the network for a brief period and then withdraw the cue to explore the spatiotemporal pattern of neural activity elicited by the cue presentation.As shown in Figure 2A, pyramidal cells showed relatively sparse and regular discharges in the absence of stimulus from 0 to 500 ms.At 500-700 ms, the stimulation elicited an intense firing in the cued direction.When the cue was withdrawn, this intense discharge in the preferred direction persisted and showed oscillatory activity throughout the 700-to 2000-ms delay period.We call this sustained oscillatory activity in the cued direction "oscillatory activity bump" and believe it maintains information about the cue.Pyramidal cells in other directions were unaffected during and after the cue presentation and continued to discharge sparsely and regularly.The inhibitory interneurons nFS (Fig. 2B, blue) and FS (Fig. 2B, red) showed opposite trends.In the prestimulus period (0-500 ms), the nFS discharged as sparsely and regularly as the pyramidal cell, whereas the FS barely discharged.Within 500-700 ms of stimulus presentation, the nFS stopped firing in the preferred direction, whereas the FS began to fire strongly and intensely in the preferred direction.This firing pattern continued throughout the 700-to 2000-ms delay period after stimulus withdrawal.The stimulus had little effect on the nonpreferred nFS and FS, and their spikes remained sparse or absent.Figure 4B provides additional information about the average membrane potential of the cells.It also shows examples of pyramidal cells, FS, and nFS in the preferred and nonpreferred directions.We then presented two or three directional cues to the network sequentially or concurrently.Given concurrent cues, the network simultaneously elicits two or three distinct oscillatory activity bumps that persist throughout the delay period (Fig. 2C,D), and each distinct activity bump maintains the information of one cued direction.For sequential cues, the network elicits distinct oscillatory activity bumps one after another (Fig. 2E,F).The later presented cue evokes a new persistent oscillatory activity bump and does not disrupt the persistent oscillatory activity bump elicited by the previous cue.These results show that our continuous network can store not only concurrent information but also sequential information of WM.

Labor division of interneurons in c and b oscillation
To investigate the oscillatory behavior in WM, we calculated the firing rate of pyramidal cells and the LFP given one directional cue in a single trial (Fig. 3).The population firing rate and LFP oscillate before and after the cue presentation and throughout the delay period (Fig. 3C,E).The spectrogram of the LFP shows that the network oscillates in the b band before the cue presentation but oscillates in the b and g bands during the delay period (Fig. 3D,F).Figures 4 and 5 provide an overview of the average membrane potential (Fig. 4A) and gating variables (Fig. 5A) during a single trial for different neurons or compartments.In addition, these figures provide examples of individual neurons in both preferred and nonpreferred directions (Figs.4B, 5B).
c Power increases in the cued direction Experiments (Pesaran et al., 2002) showed that the increase in g band power was sustained during the delay period in the preferred direction, while g band power did not change in the anti-preferred direction.We calculated the LFP based on the pyramidal cells close to the cued direction, which is the preferred direction, and the pyramidal cells far away from the cued direction, which is the anti-preferred direction.We found that in the preferred direction, the LFP mainly consists of b band oscillation before the cue presentation, and a strong g band oscillation appears after the cue presentation in a single trial (Fig. 7A) or averaged over trials (Fig. 7C).The power of g band increases significantly (Mpre ¼ 23.942, Maft ¼ 26.616, t ¼ À80.835, p , 0.001, d ¼ 0.784, cl ¼ [À2.739,À2.609]),but the power of b band decreases (Mpre ¼ 37.157,Maft ¼ 35.381,t ¼ 76.696,p ,0.001,d ¼ 0.737,cl ¼ [1.731,1.822];Fig. 7E).There is little increase in g band oscillation (Mpre ¼ 26.370, Maft ¼ 26.741, t ¼ À54.777, p , 0.001, d ¼ 0.479, cl ¼ [À0.385,À0.358];Fig. 7B) for the anti-preferred direction, and the b band oscillation is not decreased (Mpre ¼ 37.752, Maft ¼ 37.450, t ¼ 47.799, p , 0.001, d ¼ 0.348, cl ¼ [0.289,0.314];Fig. 7F).Note that there was an early increase in g and b band power in Figure 7E,F.We think this was an artifact of the cwtft algorithm.CWTFT calculates the power of a signal by locally integrating the multiplication of signals with shiftable symmetric wavelets of different window sizes and center frequencies.In the same wavelet window, wavelet analysis is isotropic, thus unable to distinguish components before and after stimulus presentation.When the signal frequency undergoes significant changes within a certain window, the window that contains the changing signal cannot differentiate between the prechange and postchange components, showing an earlier increase in power (similar to Kaiju et al., 2017).In short, cue presentation increases the power of the g oscillation and decreases b oscillation in the preferred direction (Fig. 7E,F), which are consistent with the experimental observations (Pesaran et al., 2002).

WM load enhances the c power during the delay period
We further analyzed the spectrogram of LFP given concurrent or sequential cues.The average spectrogram of LFP over 500 trials exhibits several characteristics.First, as previously shown, the network exhibits b band oscillation before the cue presentation, and the b power tend to decreases after cue presentation (two concurrent cues:  [À12.207,À11.963])and g band oscillation maintains during the delay period.Third, the g power remains almost constant during the delay period for concurrent cues (Fig. 8A,B).In contrast, the g power increases along with the arrival of new cues in sequential situations and the power is maintained throughout the delay period (two sequential cues: F ¼ 13,057.29,p , 0.001; h 2 ¼ 0.225, three sequential cues: F ¼ 24,175.79,p , 0.001, h 2 ¼ 0.446, with multiple comparisons p , 0.001; Fig. 8C,D).We show the average g power in sequential trials in Figure 8E,F.The trajectory of g power is similar to the results shown by Lundqvist et al. (2016;their Figure 6), albeit with a slight distinction: the g band exhibited an increase following the presentation of the second stimulus before the presentation of the third stimulus.
We performed a Fourier analysis on the LFP given concurrent cues, and the results are consistent with the spectrogram analysis.g power increases with the WM load (F ¼ 71.06, p , 0.001; Fig. 9).Three cues lead to the highest g power and one cue results in the lowest g power.The mean of the g power are M1cue ¼ 0.493, M2cue ¼ 1.060, and M3cue ¼ 2.357 for 1, 2, and 3 cues, respectively (multiple comparisons p12 ¼ 0.002, p13 , 0.001, p23 , 0.001; Fig. 9).The dependence of g power on WM load can be explained by network activity.Strong excitation induced by cue presentation activates high-threshold FS and the interaction between pyramidal cells and FS results in g oscillation.One cue elicits one distinct activity bump with g oscillation and more cues elicit more distinct activity bumps with g oscillation (Fig. 2), suggesting that more neurons participate in g oscillation and result in higher g power.

Brief bursts of narrow-band c and b oscillation in a single trial
Although the raster plots of spiking in Figure 2 show oscillatory activity, the single spectrogram of LFP (Figs. 3B, 7A,B) shows a different scenario: g and b oscillations occurred irregularly in the form of short bursts that are narrow and variable in frequency.This observation is similar to a previous monkey experiment (Lundqvist et al., 2016).Here, we identified the g and b bursts using the algorithm proposed in the reference (Lundqvist et al., 2016).We identified each burst in 500 trials with sequential cues and plotted each burst as a point in Figure 10A.We also calculated the b and g burst rates (Fig. 10B).The b burst rate decreases slightly after cue presentation (Mpre ¼ 0.0207, Maft ¼ 0.0196, t ¼ 21.08, p , 0.001; Fig. 7F) and remains almost constant during the delay period (Fig. 10B).The g burst rate increases when the cue is presented (Mpre ¼ 0.0570, Maft ¼ 0.0685, t ¼ À48.93, p , 0.001; Fig. 7F) and increases with increasing WM load (F ¼ 229.69, p , 0.001, multiple comparisons p , 0.001; Fig. 10B).During the delay period, the g burst shows a slightly increasing trend, which is consistent with the experimental observation by Lundqvist et al. (2016; their Figure 8B).

Discussion
We have developed a working memory model that utilizes a continuous attractor structure based on the bistable excitatory-inhibitory microcircuit consisting of FS, nFS, and pyramidal cells in PFC.The main contributions of our work are threefold.First, the proposed model can maintain the information of sequentially presented cues and concurrently presented cues.Previous oscillatory and dynamic models collapse in the presence of concurrent information because these models rely on the oscillatory phase or timing to encode the information (Constantinidis et al., 2018).The persistent activity model has not been used to manipulate sequential information (Edin et al., 2009;Wei et al., 2012).Second, we identified the mechanism underlying the enhancement of g power in the WM task.Cue presentation elicits strong excitatory recurrent synaptic currents that activate high-threshold FS.The interaction between pyramidal cells and FS leads to g oscillations and enhances g power during the delay period.Third, we demonstrated the division of labor between high-threshold FS and nFS for oscillation in the WM.The interaction between Figure 9.The normalized power of local field potential.g power induced by two concurrent cues (red line) is higher than that induced by one cue (blue line) but lower than that resulting from three concurrent cues (black line).nFS and pyramidal cells leads to b oscillation, while the interaction between FS and pyramidal cells leads to g oscillation.
We showed that the transient and sparkling g and b burst in a single spectrogram of LFP can be approximated from the persistent oscillatory activity of the network (Fig. 3B), and the g burst rate increases with the WM load (Fig. 8B).The discrete and sparkling g and b burst in a single LFP spectrogram has been considered as evidence that WM is manifested through transient or discrete oscillatory dynamics, rather than sustained firing (Lundqvist et al., 2016).However, our results indicate that transient and sparkling g and b bursts in a single LFP spectrogram may indeed result from persistent oscillatory firing.Therefore, our results favor the theory of persistent activity underlying WM.
In this work, we have demonstrated the division of labor between FS and nFS in generating g and b oscillation.The interplay between high-threshold FS and pyramidal cells results in the g oscillation, whereas the interaction between nFS and pyramidal cells leads to the b oscillation.Before the cue presentation in WM tasks, the weak recurrent excitation from spontaneous discharges of pyramidal cells can activate low-threshold nFS but not the high-threshold FS, resulting in b oscillation.Cue presentation evokes strong activity of pyramidal cells and strong recurrent excitation to activate FS.As a result, the network exhibits g oscillation.This observation is consistent with experimental results that optogenetic activation of FS induces g rhythm (Cardin et al.,2009) and the spiking of somatostatin (SOM) and parvalbumin (PV) cells differentially correlates with b and g oscillations and activation of PV cells enhances g oscillation in V1 area (Chen et al., 2017).
We computationally identified that g power can be enhanced by recruiting more high-threshold FS because of the stronger excitation evoked by the cue presentation in WM.This mechanism may be common to other cognitive processes, such as attention and sensory information processing in perceptual decision-making.It has been extensively observed in cognitive processes such as WM (Jensen et al., 2007;Yamamoto et al., 2014), attention (Engel et al., 2001;Fries et al., 2001;Jensen et al., 2007;Fries, 2009;Gregoriou et al., 2009;Kim et al., 2016), and sensory information processing in perceptual decisionmaking (Frien et al., 2000;Siegel and König, 2003;Kayser and König, 2004;Henrie and Shapley, 2005;Liu and Newsome, 2006;Berens et al., 2008).Studies had observed that the selected neuronal population increased their firing rate when attention was directed to their receptive field (Luck et al., 1997;McAdams and Maunsell, 1999;Bichot et al., 2015;Thiele et al., 2016), implying strong excitation and potentially recruiting more FS neurons.g oscillations in the sensory cortex are often considered as a proxy for the encoding of sensory evidence during perceptual decision-making because of visual g band activity and its dependence on stimulus strength and various stimulus features (Frien et al., 2000;Siegel and König, 2003;Kayser and König, 2004;Henrie and Shapley, 2005;Liu and Newsome, 2006;Berens et al., 2008).Based on these, we hypothesize that sensory input would evoke a stronger excitation in the sensory cortex, which could recruit more FS and increase the g power.
Worth noting that the synapses from pyramidal cells to the nFS are dynamic synapses with short-term potentiation (STP).STP plays an essential role in the modulation of g and b bursts.In the biological brain, there is shortterm synaptic facilitation between pyramidal cells and nFS and short-term synaptic depression between pyramidal cells and FS.These short-term plasticities allow the brain to respond differently to afferents of different durations.Brief currents cause g oscillations, while prolonged currents cause b oscillations.Based on this mechanism, these observations are reproduced in previously published work (Feng et al., 2021).In our simulation, we only briefly present the stimulus to the model, resulting in oscillations at the g frequency, and short-term synaptic facilitation between pyramidal cells and nFS helps control excitability in our model.The synapses equipped with STP can temporarily increase inhibition by increasing the effective excitation from the pyramidal cells to the nFS and expanding the range of parameters to keep excitation and inhibition in balance.Too much excitation can lead to the merging of different information (Wei et al., 2012) or false memories (Edin et al., 2009), while too much inhibition can prevent the network from maintaining information.Suppose the synapses from pyramidal cells to the nFS are fixed without STP; the network is prone to produce spurious bursts of activity because of too much excitation or to forget the information because of too much inhibition.At present, we have not considered the effects of prolonged stimulus.It would be interesting to investigate how the brain adapts and responds to prolonged stimulus in the future.

Figure 1 .
Figure1.Model structure and simulation protocol.A, The network consists of two-compartment pyramidal cells, fast-spiking (FS), and nonfast-spiking (nFS) interneurons.Pyramidal and inhibitory cells are uniformly placed on a ring and labeled by their preferred directions.B, Simulation protocol.In concurrent-cue tasks, an array of cues is presented to the network from 500 to 700 ms, followed by a delay period of 1300 ms (upper panel).In sequential-cue tasks, cues were presented one after another with an interval of 100 ms, and each cue lasted for 200 ms.

Figure 2 .
Figure2.Raster plots of spiking of neurons in the network.A, Raster plots of spiking of pyramidal cells given one cue.B, Raster plots of spiking of low-threshold nFS (blue) and high-threshold FS (red) given one cue.C, Two direction cues S11S2 were concurrently presented to the network during the cue period (500-700 ms).D, Three directional cues S11S21S3 were concurrently presented to the network during the cue period.E, Two directional cues (S1 and S2) were sequentially presented to the network.F, Three directional cues (S1, S2, and S3) were sequentially presented to the network.The horizontal axis is the time in milliseconds, and the vertical axis is the preferred direction of the pyramidal cell.s1-s3 represent the indices of the cue, and the dotted lines denote the cue presented interval(s).Please note that these symbols have the same interpretations in the following figures.

Figure 3 .
Figure 3.The oscillation in the network.A, The power of raw simulated LFP of the network.The black dots plot the power against frequency on a logarithmic scale for LFP and has a power-law distribution.The red curve represents the 250-point average smooth curve of the log-scaled power and the log-scaled frequency.B, The power of LFP normalized by the red curve in panel A. C, Example of the population firing rate of pyramidal cells in a single trial.D, Example of the spectrogram of normalized LFP shown in panel C based on Morlet wavelet.E, Approximated LFP in single trial as in panel C. F, Zoom-in on the approximated local field potential around the time of g and b burst as seen in the spectrogram in panel D. Blue curve shows the local field potential filtered at 42 Hz, and red curve shows local field potential filtered at 17 Hz.

Figure 4 .
Figure 4. Voltage traces from neurons.A, Average membrane potential for pyramidal dendrites (top left), somas (top right), FS cells (bottom left), and nFS cells (bottom right) in response to one cue.Dotted lines indicate the presence of the stimulus.B, Example of single-cell membrane potential for pyramidal dendrites (top left), somas (top right), FS cells (bottom left), and nFS cells (bottom right) in the preferred and nonpreferred direction in response to a cue.The voltage trace of neurons in the preferred direction is represented by the black line and those in the nonpreferred direction by the blue line.

Figure 5 .
Figure 5. Gating variables form neurons. A, Average AMPA, NMDA, and GABA gating variables associated with presynaptic pyramidal dendrites, somas, FS cells and nFS cells in response to one cue.Dashed lines indicate the presence of the stimulus.B, Example of single cell gating variable in response to one cue.The gating variables of individual neurons are shown at the same positions as the corresponding average gating variables in A. The gating variables of neurons in the preferred direction are represented by the black line and those in the nonpreferred direction by the blue line.

Figure 6 .
Figure 6.Labor division of FS and nFS interneuron.A, The total synaptic input to FS close to the cued direction.B, The total synaptic input to FS far away from the cued direction.C, The total synaptic input to nFS close to the cued direction.D, The total synaptic input to nFS far away from the cued direction.E, The schematic diagram of labor division of interneurons in the g band and b band oscillation.Interaction between pyramidal cell and nFS results in b band oscillation, while interaction between pyramidal cell and FS leads to g oscillation.

Figure 7 .
Figure 7. Oscillation power in the preferred direction and anti-preferred direction.A, The spectrogram for LFP approximated from the pyramidal cells preferring the cued direction in a single trial.B, The spectrogram for LFP approximated from the pyramidal cells far away from the cued direction in a single trial.C, The average of spectrogram for preferred direction over 900 trials.D, The average of the spectrogram for anti-preferred direction over 900 trials.E, The average normalized power (b band:10-30 Hz, g band: 35-100 Hz) for the preferred direction.The red line denotes the g power, and the blue line denotes the b power.F, The average normalized power for nonpreferred direction.The red line denotes the g power, and the blue line denotes the b power.

Figure 8 .
Figure 8.The average spectrogram of LFP given concurrent or sequential cues over 500 trials (A-D) and the average normalized g power for sequential cue conditions (E, F).A, The average spectrogram for two concurrent cues.B, The average spectrogram for three concurrent cues.C, The average spectrogram for two sequential cues.D, The average spectrogram for three sequential cues.E, The average normalized g power for two sequential cues.F, The average normalized g power for three sequential cues.

Figure 10 .
Figure 10.The oscillatory bursts.A, The raster plot shows the occurrence of oscillatory bursts in 500 trials of simulation.B, The rate of b and g bursts.