Skip to main content

Main menu

  • HOME
  • CONTENT
    • Early Release
    • Featured
    • Current Issue
    • Issue Archive
    • Collections
    • Podcast
  • ALERTS
  • FOR AUTHORS
    • Information for Authors
    • Fees
    • Journal Clubs
    • eLetters
    • Submit
  • EDITORIAL BOARD
  • ABOUT
    • Overview
    • Advertise
    • For the Media
    • Rights and Permissions
    • Privacy Policy
    • Feedback
  • SUBSCRIBE

User menu

  • Log in
  • My Cart

Search

  • Advanced search
Journal of Neuroscience
  • Log in
  • My Cart
Journal of Neuroscience

Advanced Search

Submit a Manuscript
  • HOME
  • CONTENT
    • Early Release
    • Featured
    • Current Issue
    • Issue Archive
    • Collections
    • Podcast
  • ALERTS
  • FOR AUTHORS
    • Information for Authors
    • Fees
    • Journal Clubs
    • eLetters
    • Submit
  • EDITORIAL BOARD
  • ABOUT
    • Overview
    • Advertise
    • For the Media
    • Rights and Permissions
    • Privacy Policy
    • Feedback
  • SUBSCRIBE
PreviousNext
Research Articles, Behavioral/Cognitive

Trial-to-Trial Variability of Spiking Delay Activity in Prefrontal Cortex Constrains Burst-Coding Models of Working Memory

Daming Li, Christos Constantinidis and John D. Murray
Journal of Neuroscience 27 October 2021, 41 (43) 8928-8945; DOI: https://doi.org/10.1523/JNEUROSCI.0167-21.2021
Daming Li
1Department of Physics, Yale University, New Haven, Connecticut 06511
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
Christos Constantinidis
2Department of Biomedical Engineering, Vanderbilt University, Nashville, Tennessee 37235
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
John D. Murray
1Department of Physics, Yale University, New Haven, Connecticut 06511
3Department of Psychiatry, Yale University, New Haven, Connecticut 06511
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
  • ORCID record for John D. Murray
  • Article
  • Figures & Data
  • Info & Metrics
  • eLetters
  • PDF
Loading

Abstract

A hallmark neuronal correlate of working memory (WM) is stimulus-selective spiking activity of neurons in PFC during mnemonic delays. These observations have motivated an influential computational modeling framework in which WM is supported by persistent activity. Recently, this framework has been challenged by arguments that observed persistent activity may be an artifact of trial-averaging, which potentially masks high variability of delay activity at the single-trial level. In an alternative scenario, WM delay activity could be encoded in bursts of selective neuronal firing which occur intermittently across trials. However, this alternative proposal has not been tested on single-neuron spike-train data. Here, we developed a framework for addressing this issue by characterizing the trial-to-trial variability of neuronal spiking quantified by Fano factor (FF). By building a doubly stochastic Poisson spiking model, we first demonstrated that the burst-coding proposal implies a significant increase in FF positively correlated with firing rate, and thus loss of stability across trials during the delay. Simulation of spiking cortical circuit WM models further confirmed that FF is a sensitive measure that can well dissociate distinct WM mechanisms. We then tested these predictions on datasets of single-neuron recordings from macaque PFC during three WM tasks. In sharp contrast to the burst-coding model predictions, we only found a small fraction of neurons showing increased WM-dependent burstiness, and stability across trials during delay was strengthened in empirical data. Therefore, reduced trial-to-trial variability during delay provides strong constraints on the contribution of single-neuron intermittent bursting to WM maintenance.

SIGNIFICANCE STATEMENT There are diverging classes of theoretical models explaining how information is maintained in working memory by cortical circuits. In an influential model class, neurons exhibit persistent elevated memorandum-selective firing, whereas a recently developed class of burst-coding models suggests that persistent activity is an artifact of trial-averaging, and spiking is sparse in each single trial, subserved by brief intermittent bursts. However, this alternative picture has not been characterized or tested on empirical spike-train data. Here we combine mathematical analysis, computational model simulation, and experimental data analysis to test empirically these two classes of models and show that the trial-to-trial variability of empirical spike trains is not consistent with burst coding. These findings provide constraints for theoretical models of working memory.

  • bursting
  • computational model
  • Fano factor
  • prefrontal cortex
  • spike-train statistics
  • working memory

Introduction

Working memory (WM) is a core cognitive function for temporary internal maintenance, and manipulation, of information across seconds-long delays. A key neuronal correlate of WM found by many studies is memorandum-selective spiking activity during mnemonic delays in neurons recorded from PFC of monkeys performing WM tasks (Funahashi et al., 1989; Miller et al., 1996; Romo et al., 1999; Leavitt et al., 2017a). These common observations have motivated an influential class of theoretical models proposing neurobiological mechanisms for how PFC circuits could produce selective spiking activity patterns that can persist over seconds to support WM (Amit and Brunel, 1997; Wang, 1999, 2001; Lim and Goldman, 2013).

Recently, alternative proposals have challenged the role of persistent activity in WM on the basis that stable activity patterns may be inadequate to describe the dynamics of neural activity on single trials. In particular, Lundqvist et al. (2016) contended that “this seemingly continuous delay activity may, however, reflect averaging across trials and/or neurons.” Their analysis of local field potential (LFP) recordings from PFC found intermittent bursts of γ-band power on single trials, and argued that these findings provide evidence in favor of models of WM maintenance based on intermittent bursts of neuronal activity (Lundqvist et al., 2016). Theoretical modeling studies have proposed potential mechanisms for how WM traces can be maintained without persistent spiking activity or stable activations across trials. In particular, short-term facilitation of synapses can maintain WM information, with neuronal spiking exhibiting sparse bursting through synaptic reactivation and adaptation (Mongillo et al., 2008; Lundqvist et al., 2011; Fiebig and Lansner, 2017; Mi et al., 2017).

These divergent proposals for mechanistic theories of WM function have contributed to a debate on the nature of WM delay activity in PFC (Constantinidis et al., 2018; Lundqvist et al., 2018a). Prior studies characterizing the temporal stationarity of WM delay spiking activity were mostly based on trial-averaged firing rates, and therefore did not specifically examine the stability of those activations across trials (Murray et al., 2017; Spaak et al., 2017; Cavanagh et al., 2018). The burst-coding proposal, whose support comes from across-trial instability of LFP power (Lundqvist et al., 2016), has not been directly tested at the level of single-neuron spiking activity. It is therefore necessary to create a testing framework addressing this issue. Here we propose that stability across trials is a key feature differentiating these two classes of models, and characterize the spike count trial-to-trial variability measured by Fano factor (FF). FF has been shown to be associated with the engagement of neurons in behavioral tasks (Churchland et al., 2010; Hussar and Pasternak, 2010; Chang et al., 2012; Purcell et al., 2012) and was shown in modeling studies to be higher in a clustered network with neurons transitioning between high and low activity states (Deco and Hugues, 2012; Litwin-Kumar and Doiron, 2012).

In this study, we analyzed a doubly stochastic Poisson model, which allows different levels of spiking burstiness, and demonstrate that FF is a sensitive measure that can well dissociate bursting and nonbursting regimes. We show that, if selective WM delay spiking activity were mediated by intermittent bursting, then testable implications include significant increase in FF and thus loss of stability across trials during the WM delay. Next, we simulate biophysically grounded spiking neural circuit models of WM maintenance, which can generate nonbursting persistent activity or intermittent bursting, to further confirm these results. Finally, we test the theoretical implications using single-neuron spike-train data recorded in macaque lateral PFC during three WM tasks of different modalities. In clear contrast to burst-coding model predictions, a population-level decrease in FF was observed in all the tasks, implying that across-trial stability is strengthened during delay. The magnitudes of FF were also smaller than burst-coding model predictions. The incompatibility between theoretical predictions and results from empirical spike-train data suggests that the extent to which intermittent burst spiking could contribute to WM delay activity is strongly limited.

Materials and Methods

FF

Single-neuron spike-train data consist of multiple trials of recordings of spiking times under each task condition. Spiking activity is analyzed from the perspectives of intensity and variability across trials during each task epoch, leading to the concepts of mean firing rate and FF. These measures are both calculated using spike trains from multiple trials and can be calculated for a specific time bin during the trial locked to a task event, such as stimulus onset. We divide the time axis into bins of size Δ, and count the number of spikes that fall into each time bin in each single trial. Let N(t,Δ) be the spike count in the time window (t,t + Δ) in a trial. The mean firing rate rmean at a given time bin under a stimulus condition is given by the following: rmean(t,Δ)=〈N(t,Δ)〉(1) where 〈...〉 denotes the average calculated across trials. FF at a given time bin under a stimulus condition is defined by the following: FF(t,Δ)=Var(N(t,Δ))〈N(t,Δ)〉(2) where the average and the variance are also both calculated across trials (Fano, 1947). FF is therefore a second-order statistic quantifying trial-to-trial variability. Low FF implies stability across trials, whereas high FF implies instability across trials. The FF for a Poisson process is 1.

Doubly stochastic Poisson model

Poisson spiking models are commonly used as statistical models of cortical spiking variability, and were motivated based on the observed proportional relationship between the mean and variance of neural responses (Tolhurst et al., 1981, 1983; Softky and Koch, 1993; Shadlen and Newsome, 1998; Oram et al., 1999). An inhomogeneous Poisson spiking process is a point process determined by a varying firing rate without random fluctuations, which has equal mean and variance of spike count. A doubly stochastic Poisson process is an extension to this with the underlying firing rate following another random process (Cox, 1955), which can model the spike count overdispersion in empirical data (Churchland and Abbott, 2012). This paradigm is suitable for the purpose of this study because it allows incorporation of the quiescent and bursty state transitions directly into the firing rate (Miller, 2006). In our study, the underlying firing rate follows a two-level telegraph process, as illustrated in Figure 1A.

In the two-level telegraph process model, the two hidden states, denoted L (low) and H (high), have firing rates rL and rH, respectively. The transition rates leaving the respective states are denoted μH and μL. Equivalently, we introduce the mean dwell time of each state: τH and τL, which are related to the transition rates by μH=1τH and μL=1τL. The dwell time in each state follows exponential distribution, with mean equal to τH and τL, respectively. The four parameters {rL,rH,τL,τH} fully define the spiking process, and are same across trials under one stimulus condition.

A key advantage of this simple model is that it can be solved analytically. With the model parameters defined above, one has immediately the mean firing rate as follows: rmean=τHrH + τLrLτH + τL(3)

If we call the time-dependent firing rate λ(t), and denote N(t,Δ) the spike count in the time window (t,t + Δ), then we have the following: 〈〈N(t,Δ)〉〉=∫tt+Δds〈λ(t)〉=rmeanΔ(4) where the inner average is with respect to the fluctuations of the firing rate, and the outer average is with respect to the Poisson process. By the properties of Poisson process, the second moment of the spike count is given by the following: 〈〈N(t,Δ)2〉〉=〈∫tt+Δdsλ(s)+[∫tt+Δds〈λ(s)〉]2〉=rmeanΔ + ∫tt+Δds∫tt+Δds′〈λ(s)λ(s′)〉(5)

Hence the variance of the spike count is given by the following: 〈〈N(t,Δ)2〉〉−〈〈N(t,Δ)〉〉2=rmeanΔ + ∫tt+Δds∫tt+Δds′〈δλ(s)δλ(s′)〉(6) where δλ(t)=λ(t)−〈λ(t)〉.

The autocovariance of the underlying telegraph firing rate process has the following form: 〈δλ(s)δλ(s′)〉=σ2exp(−(1τH+1τL)|s−s′|)(7) where σ2=(rH−rL)2τHτL(τH+τL)2

Therefore, the characteristic time of the exponential decay is τ=τHτLτH+τL. Direct calculation of the integral gives the following: 〈〈N(t,Δ)2〉〉−〈〈N(t,Δ)〉〉2=rmeanΔ + 2σ2τ2[exp(−Δτ)−(1−Δτ)](8)

Therefore, FF is given by the following: FF(Δ)=1 + 2σ2τ2[exp(−Δτ)−(1−Δτ)]rmeanΔ(9)

FF therefore depends on the four parameters {rL,rH,τL,τH} and on the bin size Δ.

Circuit model simulations

We simulated three biophysically based cortical circuit models to test that FF is a sensitive measure that can differentiate distinct WM mechanisms. Compte et al. (2000) built a recurrent circuit model with synaptic mechanisms for stable activity in spatial WM processes, in which leaky integrate-and-fire neurons were arranged on a ring with the angular location corresponding to the preferred cue position. Nonbursting persistent activity was maintained by strong recurrent connections. Here, we adopt a spatially uniform discrete node attractor network model based on similar mechanism applied to a visual motion discrimination decision-making task (Wang, 2002) (Model 1 below), first introduced by Brunel and Wang (2001). Nonbursting persistent activity realized by stable attractor networks typically fails to reproduce irregular spiking activity (Durstewitz and Gabriel, 2007; Renart et al., 2007; Barbieri and Brunel, 2008). To address this issue, we also included an attractor model with balanced slow excitation and fast inhibition implementing a negative derivative feedback mechanism, introduced by Lim and Goldman (2013), which generates persistent irregular firing patterns during WM delays (Model 2 below). For the generation of bursting in selective WM activity, Lundqvist et al. (2011) proposed that bursting can be realized through the competition between cellular adaptation and synaptic augmentation in a cyclic attractor network. Similarly, we added spike-frequency adaptation and short-term facilitation to Model 1 and obtained intermittent bursting (Model 3 below).

For FF analysis, we simulated each model the same timing of the task epochs as in our experimental datasets (with a foreperiod of 1000 ms, a stimulus onset of 500 ms, and a delay of 3000 ms). We then calculated FF on spiking activity from excitatory neurons in the pool that was stimulated by its preferred cue input. We calculated FF for each neuron using 10 simulated trials, to be in line with the typical number of trials recorded per stimulus condition for neurons in the experimental datasets (see Data preprocessing).

Model 1: stable attractor network

The model for generating nonbursting persistent activity is the same as described by Wang (2002) (see also Brunel and Wang, 2001), with minor modification to the parameter values. There are three excitatory pyramidal cell pools with a total of 1600 pyramidal cells, with two selective to the stimuli (240 neurons for each) and the third one nonselective (1120 neurons), as well as an inhibitory interneuron pool (400 neurons). Both pyramidal cells and interneurons are described by leaky integrate-and-fire neurons and are characterized by a resting potential V L = –70 mV, a firing threshold V th = –50 mV, a reset potential V reset = –55 mV, a membrane capacitance Cm = 0.5 nF for pyramidal cells and 0.2 nF for interneurons, a membrane leak conductance gL = 25 nS for pyramidal cells and 20 nS for interneurons, and a refractory period τref = 2 ms for pyramidal cells and 1 ms for interneurons. The corresponding membrane time constants are τm=Cm/gL=20 ms for excitatory cells and 10 ms for interneurons. Below the threshold, the membrane potential V(t) of each unit follows: CmdV(t)dt=−gL(V(t)−VL)−Isyn(t)(10) where Isyn (t) represents the total synaptic current flowing into the cell.

The network consists of pyramid-to-pyramid, pyramid-to-interneuron, interneuron-to-pyramid, and interneuron-to-interneuron connections. Recurrent EPSCs have two components mediated by AMPARs and NMDARs, respectively. External synaptic inputs send to the network all the information (stimuli) received from the outside world, as well as background noise because of spontaneous activity outside the local network. In simulations, external EPSCs were mediated exclusively by AMPARs. The total synaptic currents are given by the following: Isyn(t)=Iext,AMPA(t) + Irec,AMPA(t) + Irec,NMDA(t) + Irec,GABA(t)(11) in which Iext,AMPA(t)=gext,AMPA(V(t)−VE)sext,AMPA(t)(12) Irec,AMPA(t)=grec,AMPA(V(t)−VE)∑j=1CEwjsjAMPA(t)(13) Irec,NMDA(t)=gNMDA(V(t)−VE)(1 + [Mg2+]exp(−0.062V(t)/3.57))∑j=1CEwjsjNMDA(t)(14) Irec,GABA(t)=gGABA(V(t)−VI)∑j=1CIsjGABA(t)(15) where V E = 0 mV, V I = –70 mV. The dimensionless weights wj represent the structured excitatory recurrent connections (see below), the sum over j represents a sum over the synapses formed by presynaptic neurons j. NMDA currents have a voltage dependence that is controlled by extracellular magnesium concentration, [Mg2+] = 1 mm. The gating variables, or fraction of open channels s, are described as follows. The AMPA (external and recurrent) channels are described by the following: dsjAMPA(t)dt=−sjAMPA(t)τAMPA + ∑kδ(t−tjk)(16) where the decay time of AMPA currents is taken to be τAMPA = 2 ms, and the sum over k represents a sum over spikes emitted by presynaptic neuron j. In the case of external AMPA currents, the spikes are emitted according to a Poisson process with rate νext = 2.3 kHz independently from cell to cell. NMDA channels are described by the following: dsjNMDA(t)dt=−sjNMDA(t)τNMDA,decay + αxj(t)(1−sjNMDA(t))(17) dxj(t)dt=−xj(t)τNMDA,rise + ∑kδ(t−tjk)(18) where the decay time of NMDA currents is taken to be τNMDA,decay = 100 ms, α = 0.5 ms –1, and τNMDA,rise = 2 ms. Last, the GABA synaptic variable obeys the following: dsjGABA(t)dt=−sjGABA(t)τGABA + ∑kδ(t−tjk)(19) where the decay time constant of GABA currents is taken to be τGABA = 5 ms. The very fast rise times (<1 ms) of both AMPA and GABA currents are neglected. All synapses have a latency of 0.5 ms.

We used the following values for the recurrent synaptic conductances (in nS) in the N = 2000 neurons network: for pyramidal cells, gext,AMPA=2.07,grec,AMPA=0.05,gNMDA=0.165, and gGABA = 1.3; for interneurons, gext,AMPA=1.62,grec,AMPA=0.04,gNMDA=0.13, and gGABA = 1.0.

The coupling strength between a pair of neurons is prescribed according to a “Hebbian” rule: the synapse is strong (weak) if in the past the two cells tended to be active in a correlated (anticorrelated) manner. Hence, inside a selective pool, wj=w+, where w+ > 1 is a dimensionless parameter that is equal to the relative strength of potentiated synapses with respect to the baseline. We used w+ = 1.84. Between two different selective pools, and from the nonselective pool to selective ones, wj=w−, where w– < 1 measures the strength of synaptic “depression.” We used w– = 0.852. Other connections have wj = 1. We chose the stimulus strength to be μ = 38 Hz and the task coherence to be 50%.

Model 2: stable attractor network with negative derivative feedback

We use the same spiking network as described by Lim and Goldman (2013), with 1600 pyramidal neurons, 400 interneurons, and 2000 neurons in the external pool (not shown). The random connection probability p is set to 0.5. The remaining parameters are the same as listed in Lim and Goldman (2013). This network implements a negative derivative feedback mechanism through balanced slow excitation and fast inhibition. A brief stimulus is given to the external pool, and the negative derivative feedback provides a corrective signal that maintained delay activity. Simulation was done by running the MATLAB code available at ModelDB (https://senselab.med.yale.edu/ModelDB/ShowModel?model=181010).

Model 3: burst-firing cyclic attractor network (Model 1 with addition of short-term facilitation and spike-frequency adaptation)

To create WM delay activity maintained by intermittent bursts, we use a synaptic model allowing a cyclic attractor. Motivated by the model of Lundqvist et al. (2011), we add spike-frequency adaptation mediated by calcium-activated potassium channels and NMDA short-term facilitation to the pyramidal leaky integrate-and-fire neurons to Model 1. First, we add a spike-frequency adaptation current to Equation 10 as follows: CmdV(t)dt=−gL(V(t)−VL)−Isyn(t)−Iadaptation(t)(20) where Iadaptation=gKCae(V−EKCae)Cae(21) dCaedt=−CaeτCae + wCae∑lδ(t−tl)(22) with l the index of spikes. The calcium ions decay exponentially at a fast timescale τCae=90 ms, and on each firing there is an amount of calcium influx wCae=0.324. Next, we add short-term facilitation to NMDA synapses between pyramidal neurons, so that Equations 17 and 18 now read as follows: dsjNMDA(t)dt=−sjNMDA(t)τNMDA,decay + αxj(t)(1−sjNMDA(t))(23) dxj(t)dt=−xj(t)τNMDA,rise + αxFjNMDA∑kδ(t−tjk)(24) and dFjNMDAdt=−FjNMDAτF + αF(1−FjNMDA)∑kδ(t−tjk)(25) with αx=300,αF=0.008, and τF = 3000 ms, such that FNMDA decays on a much longer timescale compared with adaptation. To obtain strong intermittent bursts, we modify the previously defined parameters such that w+=5.14,w−=0.5,νext=1.9 kHz, and μ = 152 Hz, with the rest being the same. In the FF analysis, for each trial, we shifted the delay activity by a random amount between 0 and 2000 ms, to account for the fact that bursts of spiking may not be phase-locked with the stimulus onset in empirical data.

WM tasks and datasets

We applied our analysis framework on previously collected single-neuron spike-train datasets recorded from monkey PFC during three WM tasks of different modalities: an oculomotor delayed response (ODR, visual) task (Constantinidis et al., 2001); a vibrotacticle delayed discrimination (VDD, somatosensory) task (Romo et al., 1999); and a spatial match/nonmatch (MNM, visual) task (Riley et al., 2018).

The ODR data were collected from 2 macaque monkeys (COD, MAR), recorded from the dorsolateral PFC (areas 8 and 46) of the left hemisphere from both monkeys. The ODR task (see Fig. 4A) has eight stimuli for angular locations (0°, 45°, 90°, 135°, 180°, 225°, 270°, 315°). The subject fixates on a central point for a period of 1 s, and a visuospatial cue of one of these spatial angles is presented for 0.5 s, followed by a 3 s mnemonic delay. After the delay, the subject makes a saccadic eye movement to the remembered location.

The VDD data were collected from 3 macaque monkeys (R13, R14, R15), recorded from the inferior convexity of the PFC of the hemispheres on the contralateral side. The VDD task data (see Fig. 4B) consist of two sets of stimuli for vibrotactile frequencies (10, 14, 18, 22, 26, 30, and 34 Hz) for Monkey R13 and (10, 14, 18, 24, 30, and 34 Hz) for Monkeys R14 and R15. The analysis results of these two sets of VDD stimuli were finally concatenated, as we were only interested in characterizing spiking trial-to-trial variability with respect to tuning, while the exact frequencies were not important. After a 1 s foreperiod, the subject receives a 0.5 s vibrotactile stimulus of variable mechanical frequency (cue 1, f1) to the finger, followed by a 3 s mnemonic delay. After the delay, a second stimulus (cue 2, f2) is presented and the subject reports, by level release, which stimulus had a higher frequency.

The MNM data were collected from 3 male rhesus monkeys (ADR, ELV, NIN), recorded from areas 8a, 8b, 9, 10, 9/46, 45, 46, and 47/12 of the lateral PFC of the right hemispheres of ADR and NIN, and from both hemispheres of ELV. In the MNM task (see Fig. 4C), a 2° square stimulus was presented in one of nine possible locations (labeled by angular locations similar to ODR, plus one at the center) arranged in a 3 × 3 grid with a 10° distance between adjacent stimuli. The subject fixates for a period of 1 s, followed by a stimulus onset of duration 0.5 s. After a mnemonic delay of 1.5 s, a second stimulus is shown for 0.5 s either matching the first in location or of a different location. After a second delay of 1.5 s, the subject makes a saccade to one of the two choice targets representing match and nonmatch shown respectively on an additional display.

In those experiments, typically 5-20 trials (correct and incorrect) for each stimulus were recorded from each neuron. In each dataset, neurons were recorded from arrays of multiple electrodes, and neurons were sampled and recorded in an unbiased fashion without regard to their response or tuning properties (Romo et al., 1999; Constantinidis et al., 2001; Riley et al., 2018). Because neurons in these datasets were not preselected for specific tuning properties, bias in characterizing population activity is minimized.

Data preprocessing

All the data manipulations and statistical tests were done using scripts written in Python. We followed a procedure similar to what was described by Barak et al. (2010) for filtering neurons for analysis. We only included correct trials in the analysis. Some trials, particularly in the VDD and some in the MNM datasets, contain unbiologically short interspike intervals (<1 ms). This may be because of factors, such as malfunctioning of the electrodes, or that some of the trials may contain spikes from multiple units such that they partially overlap. For that reason, as an upstream step, trials were excluded if > 0.1% of the interspike intervals were <1 ms. Neurons for which >10% of the trials met this criterion were removed. At this step, 0 of 820 neurons in the ODR dataset, 756 of 3479 neurons in the VDD dataset, and 12 of 1305 neurons in the MNM dataset were removed. Then, for each dataset, a neuron was removed if its mean firing rate was <1 sp/s (spike per second) in every task epoch (foreperiod, cue, delay) under any stimulus condition. This step only removed neurons that were not responsive to the tasks but did not bias away from potential sparse representations. At this step, 34 of 820 neurons in the ODR dataset, 719 of 2723 neurons in the VDD dataset, and 103 of 1293 neurons in the MNM dataset were removed. The VDD dataset contains many unstable recordings that result in a drifting firing across trials. To exclude such neurons, stability of the recordings was assessed by performing a t test on the spike counts of the first versus the second half of the trials during the forperiod. Neurons for which p < 0.05 were removed; 710 of 2004 neurons in the VDD dataset were removed at this step. The ODR and MNM datasets suffer much less from this issue, and only <5 such neurons were excluded by visual inspection when performing the analysis. At the end, in each dataset, neurons with at least 8 correct trials per stimulus condition were finally kept for analysis. All the above minimal filtering steps were only based on data quality, but not related to tasks or tuning, and should not introduce sampling bias to the study. This procedure eventually yielded 642 neurons for the ODR task, 609 neurons for the VDD task, and 1009 neurons for the MNM task.

Data analysis

When analyzing the experimental datasets, we divided the time axis into bins of Δ = 250 ms, which is relatively long compared with proposed burst durations. If there was no spike in a given time bin in any trial under a certain stimulus condition, then that time bin was not included in the analysis. Such time bins only exist under the least preferred stimulus and is of small fraction; therefore, not much bias is introduced. In the analysis of delay activity, the first time bin immediately after the cuing period was not included to account for the potential nonzero relaxation time. ODR and MNM neurons with stimulus-selective delay activity were selected using one-way ANOVA (p < 0.05) over spike counts across stimulus conditions. VDD neurons with stimulus-selective delay activity were selected by testing a linear dependence between the mean firing rates during delay and the stimulus frequency strengths (p < 0.05). Neurons with at least consecutive 500 ms of selective activity during delay were selected in this procedure, which we also called neurons with well-tuned delay activity, and the analysis was done only on the time bins with stimulus-selective activity. In the Results section, when referring to “delay,” we mean these time bins with stimulus-selective activities of the delay. Under each stimulus condition, rmean and FF were computed during foreperiod and in the stimulus-selective time bins of the delay, respectively, then were averaged across time bins to yield single scalars. For each single neuron, the stimulus under which rmean was highest was called the preferred stimulus, whereas that under which rmean was smallest was called the least preferred stimulus. We compared FF under these two stimulus conditions during delay to the baseline mean FF during the foreperiod. We defined neurons with increased FF under the preferred stimulus those whose FF was higher during delay under the preferred stimulus than both in the foreperiod and during delay under the least preferred stimulus. We also counted neurons with FF and rmean significantly positively correlated across stimulus conditions (Pearson r > 0, p < 0.05), by concatenating all the time bins with well-tuned activity during delay.

To take into account of the temporal dynamics of coding and FF, we also compared FF in each single time bin during delay to the foreperiod baseline, in which case the preferred and the least preferred stimulus were defined according to the firing rate in each time bin.

When analyzing the simulated data of Models 1 and 3, we computed FF during foreperiod and delay, only for the neural pool that was activated by the stimulus, without discarding the time bin immediately after the cuing period. When analyzing the simulated data of Model 2, we computed FF during the last 500 ms of the foreperiod, and the last 3000 ms during the delay, for pyramidal neurons with nonzero firing rates during these periods. The time bin size used was still 250 ms.

Data and code availability

Data and code for simulation and analysis are available on request from the authors.

Results

Types of stability of spiking activity during WM

There are various candidate neural mechanisms for WM delay activity (Barak and Tsodyks, 2014; Chaudhuri and Fiete, 2016; Murray et al., 2017). For instance, stable attractor networks implementing a positive feedback mechanism allow excitation to reverberate in a recurrent network, thereby storing information (Compte et al., 2000; Wang, 2002). Burst-coding models, on the other hand, rely on synaptic mechanisms, such as synaptic facilitation with information encoded in calcium kinetics (Lundqvist et al., 2011). Attractor networks are not limited to discrete or one-dimensional attractors. For instance, chaotic attractor models can generate more complex patterns subserved by highly irregular dynamic attractors (Pereira and Brunel, 2018). There are also transient memory mechanisms, such as the feedforward chains, where the signal decays rapidly at each unit but lifetime becomes longer when passed along the line (Goldman et al., 2008).

These models can generate memorandum-selective neural states that account for different types of stability of WM. We here propose that models of WM delay activity can be qualitatively assessed from two distinct properties (Table 1). One property is stationarity across time. This can be assessed by analysis of the trial-averaged peristimulus time histogram (PSTH), with stationarity defined such that the mean firing rate during delay keeps roughly at a constant level. The other property is stability across trials, which is assessed by the variability of spiking activity across trials. These two properties provide axes that allow us to classify existing WM models, and guide us choosing a measure to dissociate divergent mechanisms.

View this table:
  • View inline
  • View popup
Table 1.

Qualitative characterization of various WM models in terms of stationarity across time and stability across trials

For instance, the stable attractor models, the burst-coding models, and the chaotic attractor models of WM are all featured by stationarity across time (Compte et al., 2000; Lundqvist et al., 2011; Pereira and Brunel, 2018), whereas the feedforward chain model of WM exhibits nonstationarity through its temporal variations of trial-averaged firing rate (Goldman et al., 2008). On the other hand, the burst-firing and the chaotic attractor models are featured by the sparseness of spiking within single trials, which imply that the bursts of spikes are mostly not aligned across trials, resulting in low levels of across-trial stability. In contrast, the stable attractor model and the feedforward chain model exhibit higher levels of stability across trials, as within each single trial, the stable attractor is characterized by sustained activity, and the feedforward chain model is characterized by consistent transient dynamics following stimulus onset.

PSTH time courses with and without strong temporal variations have both been widely observed in empirical PFC recordings (Constantinidis et al., 2001; Brody et al., 2003; Shafi et al., 2007; Tiganj et al., 2018), and these findings can test the within-trial stationarity of neuronal responses in relation to theoretical models (Murray et al., 2017). Because PSTH analyses rely on trial-averaging, they are not able to test the stability of neuronal firing across trials, which can provide important constraints to theoretical proposals for burst-firing. We therefore focus in this study on characterizing across-trial stability in PFC recordings and in computational models.

FF as a measure of trial-to-trial variability

FF is a statistical measure of variability in counts, defined as the ratio of the variance to the mean (Fano, 1947; Shadlen and Newsome, 1998; Miller, 2006). Applied to neuronal spike counts across trials within a trial-defined time bin, FF has been applied in multiple studies to characterize spiking variability, and task-related changes in FF have been found to reflect neural engagement in cognitive tasks (Churchland et al., 2010; Hussar and Pasternak, 2010; Chang et al., 2012; Purcell et al., 2012). Consistent with these prior studies, here we define FF with the mean and variance of spike count taken across trials at each time bin (Eq. 2). Such a definition allows us to assess trial-to-trial variability and thus stability across trials during WM delay, while accounting for temporal dynamics of firing rates within a trial. We note that our more-common use of FF to measure trial-to-trial variability is distinct from using FF to characterize within-trial temporal dynamics. For example, Shafi et al. (2007) defined an intratrial FF with the mean and variance taken across time bins within a trial, to quantify how dynamic activity is within a trial. These different uses of FF highlight the critical distinction made in the previous section, between stationarity across time and stability across trials.

Poisson processes are widely chosen spiking models based on the observed proportional relationship between the mean and variance of spiking responses in cortical neurons (Tolhurst et al., 1981, 1983; Softky and Koch, 1993; Shadlen and Newsome, 1998; Oram et al., 1999). A simple extension for neuronal spiking is an inhomogeneous Poisson model, determined by a time-varying firing rate without randomness, which has equal mean and variance of spike counts within the same time window, and thus FF equal to 1. In a stable attractor network, excitatory input stabilizes a specific attractor, so that spiking variability is quenched. This is likely to result in sub-Poisson statistics and FF <1 (Durstewitz and Gabriel, 2007; Renart et al., 2007; Roudi and Latham, 2007; Barbieri and Brunel, 2008). Burst-coding models, on the other hand, have two sources of variability: the variability of the firing rates and the variability of spikes (Miller, 2006; Churchland et al., 2011). Therefore, the extra variability is likely to result in super-Poisson statistics and FF >1. This motivates us to build a doubly stochastic Poisson spiking model with the firing rate randomly transitioning between high and low activity states to capture the main characteristics of the alternative burst-coding proposal.

A doubly stochastic Poisson spiking model

We first aimed to build a parsimonious model that can be parametrically tuned to reach either the stable or burst-spiking regimes. This allows us to make predictions on how FF behaves differently in the two distinct scenarios. A doubly stochastic Poisson process is an extension to the classic Poisson model, with the underlying firing rate following another random process, which can capture firing rate fluctuations in addition to spiking variability. This paradigm is suitable for the purpose of this study because it allows incorporation of the low-activity quiescent and high-activity bursty state transitions directly into the firing rate. Burst-coding WM models hypothesize that each single trial consists of sparse sharp bursts of spiking, while during most of the time, neurons are in a baseline state of low-rate firing. According to this proposal, the memorandum selectivity of WM-related delay activity arises from modulation of these burst events (Lundqvist et al., 2016, 2018b). These proposals motivated us to translate these features into a doubly stochastic Poisson spiking model, where neuronal spiking is a Poisson process with the underlying firing rate following another two-level telegraph random process, as illustrated in Figure 1A.

Figure 1.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 1.

Two-state telegraph process as a doubly stochastic Poisson model of neuronal burst spiking. A, Neuronal spiking is modeled as a doubly stochastic Poisson process with the underlying firing rate following a random two-level telegraph process. H and L denote the high and low firing states, with respective firing rates rH and rL, transition rates leaving the states μH and μL, τH and τL are the mean dwell time in each state. B, Three sample trials of model-simulated telegraph processes with the corresponding generated spike trains. The model parameters used here are in a biophysically plausible range: rL = 5 sp/s, rH = 100 sp/s, τL = 350 ms, τH = 65 ms. Colored bursts typically occur in the high firing state. C, Any homogeneous Poisson process with constant firing rate has FF equal to 1, whereas the example doubly stochastic Poisson process, defined in B, has FF equal to 6.15. D, The model in general can generate different levels of burstiness at fixed mean firing rate rmean by tuning the four parameters. The intermittent bursting regime is characterized by low rL, high rH, and low τHτL. Less burstiness can be achieved by raising the low firing rate rL, or by increasing the ratio τHτL.

Here the two firing states L (low) and H (high or bursting) have firing rates rL and rH, respectively. In the burst-coding scenario, rL corresponds to the background firing that is homogeneous across all stimulus conditions, and should be relatively low. Bursts are most likely to occur when the neuron is in the high firing state H, with rH ≫ rL. We also introduce the mean dwell time in each state: τH and τL. Since the bursts are claimed to be sparse, we have τL ≫ τH. The four parameters {rL,rH,τL,τH} fully define the spiking process and are the same across trials under one stimulus condition. Sample trials of a two-level telegraph process with their generated spike trains are shown in Figure 1B. Based on the aforementioned considerations, we ignored possible ramping of mean firing rate rmean which occur at a slower timescale, and focused on characterizing the across-trial variability (with FF) of spike counts, which does not depend on slow modulation of the mean firing rate as a function of time from delay period onset.

This model can also produce a nonbursting regime at given mean firing rate by flexibly tuning the parameters, for instance, by raising the low-state background firing rate rL, or by increasing the mean dwell-time ratio τHτL (Fig. 1D). rmean and FF of this model can be analytically solved in the general form (with derivations shown in Materials and Methods).

Theoretical implications of burst-coding models

To investigate what the burst-coding hypothesis implies within this doubly stochastic Poisson framework, we characterized a biologically plausible parameter space subject to the following constraints motivated by the PFC neurophysiology literature. First, we varied the parameter τH between 40 and 80 ms, based on the report from Lundqvist et al. (2016) of an average LFP burst duration of 67 ms. Second, we varied the mean dwell time ratio τHτL between 0.1 and 0.3, based on reports from Lundqvist et al. (2016, 2018b) of a burst rate between slightly lower than 0.1 and roughly 0.25 during the delay epochs in the tasks studied. Third, we varied rmean between 10 and 50 sp/s. This range was selected because at a temporal resolution of 250 ms, the mean firing rates under the preferred stimulus among all the selected neurons in the three WM tasks (for details, see Materials and Methods) range from 1.2 to 98 sp/s (spike per second), with most of the rates between 10 and 50 sp/s after dropping the two 10% quantile tails. The overall mean of rmean under the preferred stimulus for each task is as follows: 20.9 sp/s for ODR, 26.8 sp/s for VDD, and 14.7 sp/s for MNM. Based on these values, we chose 20 sp/s as the default parameter for rmean in our model simulations. Fourth, we varied the parameter rL between 1 and 10 sp/s, considering that rL corresponds to such a default state that does not modulating the cognitive process and is assumed to be homogeneous in the burst-coding proposal. This range is in line with our spike-train data, in that the median of the mean firing rates during the delay under the least preferred stimulus is ∼3 sp/s for the ODR task, 8 sp/s for the VDD task, and 2 sp/s for the MNM task. Fifth, based on Equation 3, we solved for a range of rH between 70 and 200 sp/s, which corresponds to τH = 65 ms, rmean = 20 sp/s, rL = 5 sp/s, and τHτL between 0.1 and 0.3.

If WM delay were modulated by sharp intermittent bursts, there would be the following important implications derived from the doubly stochastic Poisson model in the burst-coding regime, which can be tested in empirical WM neuronal recordings.

Implication 1: burst strength tends to be strong if the burst coding assumptions were valid

If WM delay coding were instantiated as sparse intermittent bursts, we can estimate how strong the burst would need to be. For this purpose, we transformed Equation 3 and expressed the burst strength (firing rate in the high state) rH as a function of the mean firing rate rmean, the mean dwell time ratio τHτL and the firing rate in the low state rL as follows: rH=rmean(τHτL + 1)−rLτHτL(26)

If WM were subserved by sharp intermittent bursts, we would have rH ≫ rL and τHτL≪1. Figure 2A is a visualization of the magnitudes of rH by varying other parameters in a biophysically reasonable range consistent with burst-coding model assumptions. We fixed the firing rate in the low state rL = 5 sp/s, and varied rmean between 10 and 50 sp/s and the mean dwell time ratio τHτL between 0.1 and 0.3. It can be seen that the implied burst strength rH tends to be high, even potentially exceeding the biophysical limit when both the observed rmean is high and the bursts are sparse (low τHτL). This provides a first constraint on the burst sparsity of delay activity: the bursts must be very strong and cannot be very sparse if the observed mean firing rate is high.

Figure 2.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 2.

Theoretical predictions of the intermittent burst-coding model. A, Firing rate in the high state should be very high under the sparse intermittent burst-coding scenario, given fixed firing rate in the low state. Heatmap of rH is shown as a function of the mean firing rate rmean and the dwell-time ratio τHτL according to Equation 26, with firing rate in the low state fixed. Parameters were chosen in a biophysically reasonable range: rL = 5 sp/s is fixed; rmean is varied between 10 and 50 sp/s as typically observed in experiments; τHτL is varied between 0.1 and 0.3 as implied by sparsity. B, FF grows with burst strength rH, with rL = 5 sp/s, τH = 65 ms, τL = 350 ms, and time-bin width Δ = 250 ms fixed. C, FF grows with burst duration τH in the sparse coding regime, with rL = 5 sp/s, rH = 100 sp/s, τHτL=0.15, and Δ = 250 ms fixed. D, FF grows with burst rate τHτL in the sparse coding regime, with rL = 5 sp/s, rH = 100 sp/s, τL = 350 ms, and Δ = 250 ms fixed. E-G, Magnitudes of FF tend to be large within a biophysical range, Δ = 250 ms, varying (rH, τH), (rL, τL), and (τL, τH) respectively, with the mean firing rate rmean = 20 sp/s fixed. H, FF increases monotonically with the mean firing rate. When varying the burst rate and equivalently τH between 30 and 100 ms, we set rL = 5 sp/s, rH = 80 sp/s, τL = 350 ms, and Δ = 250 ms. When varying the burst strength rH between 45 and 120 sp/s, we set rL = 5 sp/s, τH = 65 ms, τL = 350 ms, Δ = 250 ms. I, FF is still positively correlated with the mean firing rate after adding refractoriness into model simulation by allowing only the interspike intervals >2 ms. Varying burst rate: Pearson r = 0.65 (p < 10–120). Varying burst strength: Pearson r = 0.88 (p < 10–300). The parameters used in simulation are the same as in H.

Implication 2: FF tends to be very large if the burst-coding assumptions were valid

It can be shown that FF as a function of the model parameters in general is given by the following: FF(Δ)=1+2σ2τ2[exp(−Δτ)−(1−Δτ)]rmeanΔ(27) with Δ the time bin size, σ2=(rH−rL)2τHτL(τH+τL)2 and τ=τHτLτH+τL (for derivation, see Materials and Methods). Here FF must be >1 because of the extra variability because of the fluctuations of the underlying firing rate, compared with homogeneous Poisson model whose FF is equal to 1 (Fig. 1C). Equation 27 has important neurophysiological implications. In the fast rate fluctuations limit (Δ ≫ τ) that we are interested in, FF saturates to 1+2σ2τrmean. Exhibiting sharp intermittent bursts require rH ≫ rL and τH ≪ τL, so that τ is close to τH. We chose a time bin size of 250 ms for analysis, which is long enough compared with the mean burst duration. Figure 2B-I shows the properties of FF implied by Equation 27, in a biophysical parameter range. If we fix rL = 5 sp/s and rmean to be what is usually observed in empirical data (e.g., 20 sp/s) and vary each pair of the free model parameters, burst-coding assumptions imply large magnitudes of FF in a biophysically plausible regime (Fig. 2E-G). Relatively smaller FF only appears at regions with stronger baseline firing and less sparseness, such that the high and low state discharge patterns are not as sharply distinguishable. Overall, these results suggest that if burst-coding assumptions applied, across-trial stability would be low as reflected in a high FF.

Implication 3: FF should be positively correlated with the mean firing rate if the burst-coding assumptions were valid

Equation 27 also implies that FF and rmean are positively correlated under the burst-coding assumptions. With other model parameters fixed, FF grows with burst strength rH, burst duration τH, and mean dwell time ratio τHτL in the biophysically plausible regime (Fig. 2B-D). Varying burst strength through rH and burst rate through τH with other parameters held constant trace trajectories in the (FF, rmean) plane. It can be seen that FF is positively correlated with rmean (Fig. 2H), which means that FF is larger for the neuron's preferred stimuli. We also tested the robustness of this result under the introduction of refractoriness following each spike, by running simulations of the model allowing only interspike intervals >2 ms, and still found the significant positive correlation between FF and firing rate (Fig. 2I). A natural consequence of these results is that burst-coding models necessarily imply an increase in FF during the delay compared with the foreperiod.

Trial-to-trial variability measured by FF dissociates distinct circuit mechanisms

The doubly stochastic Poisson model is based on the mathematical abstraction of the burst-coding assumptions. Although it captures the most important properties of the burst-coding mechanism, the simplicity of the model does not include underlying biophysical processes included in circuit models of WM maintenance. We therefore sought to validate that our FF measures can be used to distinguish different dynamic regimes of WM delay activity in influential circuit models. We simulated three spiking circuit models to generate both bursting and nonbursting WM delay activities, and confirm that FF is a sensitive measure that can well dissociate WM circuit models instantiating distinct mechanisms.

Attractor networks have been the most popular class of models for information storage (Wang, 2001). To generate nonbursting persistent activity, we used the spatially uniform stable attractor network model developed by Wang (2002) with minor modifications, focusing on the delay activity in the neural pool that was activated under the stimulus (Model 1, Fig. 3A; see Materials and Methods). Information is maintained by strong recurrent connections between neurons within the same selective neuronal pool. We adopted this discrete node model architecture because it allows generation of localized delay activity in both nonbursting and bursting scenarios (Model 3 below), facilitating direct comparison.

Figure 3.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 3.

Numerical simulations confirm that FF is a sensitive measure which can dissociate distinct circuit mechanisms for WM delay activity. A, The architecture of the spatially uniform discrete node spiking network model generating persistent activity (Wang, 2002). There are three pyramidal neural pools with two (A and B) selective to stimulus and one nonselective (not shown), as well as an interneuron pool. The two selective pools have strong recurrent excitation within themselves, mediated by NMDARs and AMPARs. The two pools have lateral inhibition between them mediated by the GABAergic inhibitory population. Pool A is of interest in the study as it has elevated persistent activity under the stimulus. B, Spike raster plot across a subset of neurons (40 from each pool) in a single trial with long-lasting persistent activity in Pool A during delay (0-3 s: foreperiod; 3-4 s: stimulus onset; 4-10 s: delay). C, With nonbursting persistent activity, there is a global reduction in FF during delay compared with foreperiod. D, Spike raster plot across trials during delay (only 3 s displayed) for a single neuron with nonbursting persistent activity. The mean firing rate with chosen simulation parameters is ∼50 sp/s. E, The architecture of the network implementing a negative derivative feedback mechanism (Lim and Goldman, 2013). F, Spike raster plot across a subset of pyramidal neurons (80) in a single trial with long-lasting persistent activity realized by negative derivative feedback (0-3 s: foreperiod; 3-4 s: stimulus onset; 4-10 s: delay). The mean firing rate with chosen simulation parameters is ∼33 sp/s. G, There is a global reduction in FF during delay compared with foreperiod. H, Histogram of change in FF across neurons. FF during delay is significantly smaller than that during the foreperiod (p < 10–26). I, Intermittent bursting is obtained by adding spike-frequency adaptation and short-term facilitation on top of the model in A. J, Temporal evolution of the adaptation current and facilitation variables (0-3 s: foreperiod; 3-4 s: stimulus onset; 4-10 s: delay). K, Spike raster plot across a subset of neurons (40 from each pool) in a single trial with long-lasting stable intermittent bursting in Pool A during delay. L, With intermittent bursting, there is a global increase in FF from foreperiod to delay. M, Spike raster plot across trials during delay (only 3 s displayed) for a single neuron with intermittent bursting. The mean firing rate with chosen simulation parameters is ∼40 sp/s. *p < 0.05. **p < 0.01. ***p < 0.001.

Stable attractor networks are known to produce a reduction of spiking variability to sub-Poisson levels in neurons participating in the persistent WM activity, because excitatory and inhibitory synaptic inputs can become unbalanced in a high-activity WM state (Durstewitz and Gabriel, 2007; Renart et al., 2007; Roudi and Latham, 2007; Barbieri and Brunel, 2008). For instance, overall excitatory input can stabilize one specific attractor state, resulting in reduced variability. To address this issue, we simulated a network which realizes stable persistent delay activity through balanced slow excitation and fast inhibition (Lim and Goldman, 2013) (Model 2, Fig. 3E; see Materials and Methods). Such a network implements a negative derivative feedback mechanism which generates spikes in a persistent way, while preserving Poisson-like irregular firing which is driven by fluctuations of subthreshold mean synaptic input when excitation and inhibition cancel each other out.

To generate transient bursts, on the other hand, we adopted a synaptic mechanism for maintaining WM traces over relatively long intervals. Mongillo et al. (2008) proposed that WM signals can be maintained through the mechanism of short-term synaptic facilitation, mediated by increased residual calcium levels at presynaptic terminals. In the model proposed by Lundqvist et al. (2011), neuronal spike-frequency adaptation reduces the lifetime of high-activity WM states, which generates short bursts of spiking, whereas short-term synaptic facilitation allows synapses to maintain a selective WM trace which causes reactivation of the associated pyramidal-neuron pool. This combination of faster spike-frequency adaptation and slower synaptic facilitation generates a WM code consisting of intermittent bursting in a selective pyramidal-neuron pool. To implement this proposed mechanism in our models, here we incorporated a spike-frequency adaptation current and short-term synaptic facilitation in the NMDA synapses in the excitatory neurons to Model 1 (Model 3, Fig. 3I; see Materials and Methods). In this scenario, WM information is maintained by spike-induced changes in synaptic weights in between bursts of selective spiking activity.

We compared simulated spike trains displaying nonbursting persistent activity and intermittent bursts in the neural population that is responsive to the stimulus. With the parameters chosen, all three models are able to generate long-lasting WM delay activity (for illustrations, see Fig. 3B,F,K). Our Model 3 generates similar bursty spiking patterns to the ones shown in Lundqvist et al. (2011). Figure 3J shows the temporal evolution of the facilitation and adaptation variables. Sample single neuron raster plots were shown for Models 1 and 3, which are spatially uniform models (Fig. 3D,M). We calculated FF on the simulated data with the same time lines as in the WM tasks (see Three spike-train datasets recorded in monkey lateral PFC during WM tasks, foreperiod: 1 s, cue onset: 0.5 s, and delay: 3 s). With Model 1, we found that there was a global reduction in FF during the delay compared with foreperiod, and neurons with strongest persistent activity fired with low variability with FF much <1 (Fig. 3C). With Model 2, there was also a statistically significant shrink in FF during delay compared with foreperiod (Fig. 3G,H). In sharp contrast to Models 1 and 2, there was a dramatic increase in FF during the delay compared with foreperiod (Fig. 3L). The magnitudes of FF in the biophysical Model 3 were significantly smaller than those in the doubly stochastic Poisson model. This is because, in our mathematical model, the burst duration and the interburst intervals are completely random and follow exponential distribution, whereas synaptic biophysical models simulating real neuronal processes had much less flexibility and more constraints, including refractoriness. Nonetheless, both mathematical and biophysical models yielded the same qualitative predictions. These simulation results from a biophysically grounded circuit model further confirmed the conclusions in the previous section, thereby justifying that FF is a sensitive measure that can dissociate diverging WM mechanisms, and can be used to test on experimental data.

Three spike-train datasets recorded in monkey lateral PFC during WM tasks

We analyzed single neuron spike-train data recorded from lateral PFC of macaque monkeys in three parametric WM tasks: an ODR task (Constantinidis et al., 2001) (Fig. 4A); a VDD task (Romo et al., 1999) (Fig. 4B); and a spatial MNM task (Riley et al., 2018) (Fig. 4C). The three tasks differ in several features, allowing us to test the generality of our findings. They differ in stimulus modality (visual for ODR and MNM vs somatosensory for VDD), role of WM in guiding behavioral response (veridical report of location for ODR vs binary discrimination for VDD vs binary match/nonmatch report for MNM), and prototypical stimulus tuning curves of single PFC neurons (bell-shaped for ODR and MNM vs monotonic for VDD). Delay activity in the ODR task can be attributed to response preparation in some neurons (Funahashi et al., 1993; Markowitz et al., 2015), whereas the VDD and MNM task delay activity is free of this premotor signal as the action is not specified until after the WM delay. The tasks have a 0.5 s cue epoch followed by a 3 s (ODR and VDD) or a 1.5 s (MNM) delay epoch, which are relatively long and allow characterization of time-varying WM representations.

Figure 4.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 4.

WM tasks. A, In the ODR task, the subject fixates on a central point for a period of 1 s, and a visuospatial cue of variable spatial angle is presented for 0.5 s, followed by a 3 s mnemonic delay. After the delay, the subject makes a saccadic eye movement to the remembered location (Constantinidis et al., 2001). B, In the VDD task, after a foreperiod of 1 s, the subject receives a 0.5 s vibrotactile stimulus of variable mechanical frequency (cue 1, f1) to the finger, followed by a 3 s mnemonic delay. After the delay, a second stimulus (cue 2, f2) is presented and the subject reports, by level release, which stimulus had a higher frequency (Romo et al., 1999). C, In the spatial MNM task, the subject fixates for a period of 1 s, followed by a stimulus which is shown for 0.5 s. After a delay of 1.5 s, a second stimulus is shown for 0.5 s, either matching the first in location or of a different location. After a second delay of 1.5 s (not analyzed), the subject makes a saccade to one of the two choice targets representing match and nonmatch shown, respectively, on an additional display (not shown) (Riley et al., 2018).

FF observed in PFC inconsistent with burst-coding models

The key results we have established so far with statistical and biophysical modeling is that, if WM delay activity were encoded as sharp intermittent bursts, there should be a dramatic increase in FF during delay compared with the foreperiod, and FF should be positively correlated with rmean, with FF highest under the preferred stimulus. To test whether these predictions are true in empirical spiking data, we analyzed three single-neuron spike-train datasets recorded in lateral PFC during WM tasks described in the previous section, and characterized FF during different task epochs under each stimulus condition.

We selected neurons with at least consecutive 500 ms of stimulus-selective delay activity (for details, see Materials and Methods), and only analyzed time bins with stimulus-selective activity in correct trials. We divided task epochs into time bins of 250 ms, and computed rmean and FF in each time bin under each stimulus condition. We defined the preferred and least-preferred stimulus conditions according to rmean in the corresponding task epoch. A baseline level of FF was computed during the foreperiod by averaging across time for each neuron. Next, we measured FF averaged across trials during the delay under each stimulus condition, and compared with the baseline level. The correlation between rmean and FF was computed as the Pearson correlation by concatenating all stimulus-selective delay-period time bins.

We sought to identify neurons exhibiting task-related burstiness, and therefore classified all selected neurons with criteria based on the theoretical implications of the burst-coding models examined in the sections above. Figure 5 summarizes, for each task, the numbers of neurons with well-tuned delay activity that meet the following criteria: (1) higher FF during the delay than during the foreperiod under the preferred stimulus; (2) FF > 3 for the preferred stimulus (the minimum of FF shown in Fig. 2E-G); (3) positive correlations between rmean and FF across stimulus conditions during delay (p < 0.05) and the intersections between these categories. Burst-coding models predict that neurons should be at the intersections of these categories. It can be seen from the Venn diagram, however, that numbers of neurons at the intersections are very small compared with the number of well-tuned neurons in each task (ODR: 8 of 179; VDD: 2 of 139; MNM: 0 of 75). Even the neurons satisfying the conditions 2 or 3 alone are of small fraction in each one of the tasks. Conversely, most of the well-tuned neurons are in none of these categories (ODR: 113 of 179; VDD: 100 of 139; MNM: 41 of 75).

Figure 5.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 5.

Venn diagram showing for each task (A-C) the number of neurons categorized as follows: displaying significant positive correlation between rmean and FF, having FF > 3 under the preferred stimulus, having FF under the preferred stimulus higher than those in the foreperiod and under the least preferred stimulus (termed “Increased FF under pref stim” in the figure), and the combination of intersections between all three categories. The total number of well-tuned neurons and those that are in none of these categories are also shown. Numbers in each segment are exclusive.

Figure 6 displays, for each task, plots for a single neuron with nonbursting persistent activity (Fig. 6A,C,E), and for the single neuron with maximal bursting behaviors chosen from neurons at the intersections of the three categories defined above (Fig. 6B,D,F). In general, stationarity across time is not a universal feature of WM delay activity: the mean firing rates can have trends, such as ramping (Fig. 6A,C). We noticed that, even among the neurons chosen from the intersections of the Venn diagram, firing patterns were still not as bursty as hypothesized by burst-coding models. For instance, the neuron in Figure 6F has very low firing rates, and its sparsity was not because of burstiness but because of weak spiking activity. In Figure 6D, the neuron has high FF during the foreperiod and under the least preferred stimulus during the delay. Therefore, it is an intrinsically bursty neuron, and the burstiness is not related to the task (Compte et al., 2003; Shinomoto et al., 2009; González-Burgos et al., 2019). In addition, higher FF can be because of weak drift in overall levels of firing rate across trials.

Figure 6.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 6.

Sample neurons with their IDs in the datasets. The five panels from left to right are as follows: i, temporal evolutions of the mean firing rate under each stimulus condition, denoised by projecting to the first principal components capturing 80% of the variance of population activity of all well-tuned neurons in the same dataset (Murray et al., 2017); ii, spike raster plot under the preferred stimulus; iii, spike raster plot under the least preferred stimulus; iv, bar plot of average FF during foreperiod, during the stimulus-selective segments of the delay under the preferred stimulus, and during the stimulus-selective segments of the delay under the least preferred stimulus; and v, scatterplot of average FF versus mean firing rate across stimulus conditions during the stimulus-selective segments of the delay. A, A neuron with nonbursting stable activity in the ODR task. B, The neuron with maximal bursting behaviors in the ODR task. C, A neuron with nonbursting stable activity in the VDD task. D, The neuron with maximal bursting behaviors in the VDD task. E, A neuron with nonbursting stable activity in the MNM task. F, The neuron with maximal bursting behaviors in the MNM task.

We next characterized these FF features at the population level. We plotted the histograms of FF during foreperiod, and under the preferred and the least preferred stimulus conditions during delay for each task. We found that the majority of neurons had FF <3 in each one of the tasks (Fig. 7, top row), significantly lower than mathematical and circuit model predictions. We also created scatterplots of FF for preferred versus least-preferred stimuli during delay, and the histograms of their differences (Fig. 7, bottom row). We found that in ODR and MNM datasets, but not in the VDD dataset, FF values for the preferred stimulus were larger than those for least preferred stimulus, at relatively small effect sizes with statistical significance (two-sided Wilcoxon signed-rank test, p = 0.003, n = 179, Cohen's d = 0.24 for ODR; two-sided Wilcoxon signed-rank test, p = 0.15, n = 139, Cohen's d = 0.18 for VDD; two-sided Wilcoxon signed-rank test, p = 0.002, n = 75, Cohen's d = 0.35 for MNM).

Figure 7.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 7.

Population statistics of FF in each task. Top, Histograms of FF during the foreperiod, under the preferred stimulus during delay, and under the least preferred stimulus during delay, for all neurons with well-tuned delay activity. Bottom, Scatterplot of FF under the preferred stimulus versus under the least preferred stimulus, along with histogram of the difference in FF (preferred – the least preferred stimulus). A, Results of the ODR task. FF under the preferred stimulus was slightly greater than that under the least preferred stimulus. B, Results of the VDD task. FF under the preferred stimulus was not significantly different from that under the least preferred stimulus. C, Results of the MNM task. FF under the preferred stimulus was greater than that under the least preferred stimulus. *p < 0.05. **p < 0.01. ***p < 0.001.

We also characterized the distribution of the correlation between FF and mean firing rate (Corr(rmean,FF)) across stimulus of the entire population for each task. Similarly, in ODR and MNM datasets, the distributions are weakly biased to the positive side with large SDs, while in the VDD dataset, the distribution is centered around 0 (mean: 0.08, SD: 0.26, two-sided Wilcoxon signed-rank test: p = 0.001 for ODR; mean: −0.05, SD: 0.33, two-sided Wilcoxon signed-rank test: p = 0.23 for VDD; mean: 0.12, SD: 0.29, two-sided Wilcoxon signed-rank test: p = 0.002 for MNM). This shows that in the empirical data overall there is no strong positive correlation between rmean and FF as predicted by the burst-coding models.

As illustrated in the doubly stochastic model and the circuit models above, if WM delay activity were a result of transient bursts, there should be higher FF for the preferred stimulus during delay compared with the foreperiod baseline. We tested this through analysis of (neuron, bin) pairs during the delay for each task, to control for temporal variability of firing rates across the delay. Here the preferred and least-preferred stimuli were defined based on the firing rates in each time bin, and the FF of this bin was compared with the baseline foreperiod FF. We found that, in each task, FF values for both stimulus conditions during the delay were reduced compared with the foreperiod (Fig. 8A-C), at small effect sizes (ODR preferred stimulus: two-sided Wilcoxon signed-rank test, p = 1.4 × 10–16, n = 1151, Cohen's d = 0.15; ODR least preferred stimulus: two-sided Wilcoxon signed-rank test, p = 3.8 × 10–49, n = 1053, Cohen's d = 0.29; VDD preferred stimulus: two-sided Wilcoxon signed-rank test, p = 5.1 × 10–21, n = 726, Cohen's d = 0.21; VDD least preferred stimulus: two-sided Wilcoxon signed-rank test, p = 1.4 × 10–23, n = 726, Cohen's d = 0.32; MNM preferred stimulus: two-sided Wilcoxon signed-rank test, p = 0.004, n = 237, Cohen's d = 0.06; MNM least preferred stimulus: two-sided Wilcoxon signed-rank test, p = 1.7 × 10–7, n = 237, Cohen's d = 0.11). Combining with findings from Figure 7, we find that FF is quenched after stimulus onset in all the three tasks, while FF under the preferred stimulus is reduced by less amount compared with under the least preferred stimulus in ODR and MNM tasks. This suggests that WM delay activity is characterized by relatively high across-trial stability, which is in clear contradiction to burst-coding model predictions that are in the opposite direction.

Figure 8.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 8.

Statistics of FF of (neuron, bin) pairs in each task. Histograms of the change in FF (delay – foreperiod) under the preferred and least-preferred stimulus conditions, for all (neuron, bin) pairs with well-tuned delay activity. Empirical results contradict sharply the theoretical predictions of burst-coding models. A-C, In all three datasets, the change in FF was significantly <0 for both preferred and least-preferred stimulus conditions. *p < 0.05. **p < 0.01. ***p < 0.001.

To summarize, we have the following findings in the three empirical spike-train datasets from WM tasks. First, the magnitudes of FF were much smaller than model predictions. Second, only a small fraction of neurons displayed increased FF for the preferred stimulus, and have positive correlation between FF and rmean during delay. Third, FF values were overall reduced during WM delay compared with foreperiod baseline, although FF values for the preferred stimulus were reduced slightly less in ODR and MNM tasks. Fourth, WM delay activity is characterized by stability across trials, but not necessarily stationarity across time. Reasoning by contradiction, we conclude that coding by intermittent bursts in single-neuron spiking cannot well explain WM delay activity observed in PFC. On the other hand, FF can be sensitive to potential drift of overall firing rates across trials. However, across-trial drift, which is unrelated to bursting, would only inflate FF but not reduce it. In that sense, our estimation here provides an upper bound of FF contribution because of burstiness, which even further strengthens the constraints on bursting.

Empirical results are better explained by the nonbursting picture

In the previous section, we observed small magnitudes of FF (mostly < 3 ≈ the 95% quantile) at the population level in each task, and very few neurons had FF increasing with the mean firing rate rmean across stimulus conditions during delay. The doubly stochastic Poisson spiking model we built can provide insights on the interpretation of these results. The burst-coding model requires that the baseline firing should be weak; therefore, in terms of model parameters, we should have rL≪rH, with rL small and homogeneous across stimulus conditions. It also requires that the bursts are sparse, namely, τHτL≪1. If these conditions were not met, for instance with rL close to rH, or rL tuned with stimulus, or τHτL not low enough, FF would be lower at given observed mean firing rate level, according to Equation 27.

In Figure 9A, as an illustrative example, we depicted three sets of parameters all with the same rmean ≈ 20 sp/s. Scenario 1 corresponds to the bursting regime with high FF. Scenario 2 corresponds to the a regime with a smaller difference between rL and rH, and low FF. Scenario 3 corresponds to the a regime with burst duration longer than the low-state duration, and low FF. In Figure 9B, we reparametrized and plotted FF as a function of rmean and the scaled difference of firing rates between the two states rH−rLrmean. As shown there, the low FF (<3) region is characterized by a small difference between rL and rH, and the parameter set of Scenario 2 lies in this region. This means that the default baseline state may not be as quiescent as expected, but is at a relatively high firing rate. In Figure 9C, we reparametrized and plotted FF as a function of rmean and the mean dwell time ratio τHτL. As shown there, low FF (<3) region is characterized by very high τHτL, even far exceeding the sparse regime, and the parameter set of Scenario 3 lies in this region. This means that the actual burst duration is too long to be considered as a brief burst. Our mathematical formalism also provides an explanation for why most neurons had the same level of FF across stimulus conditions, although the firing rates were well tuned. It is possible to fix FF, and solve for rmean in terms of rL, τHτL, and τH. As shown in Figure 9D, we fixed FF (= 1.5 here, as most frequently observed in empirical data) and τH = 65 ms, and plotted rmean against rL and τHτL. Similarly, higher rmean, which corresponds to preferred stimuli, can be obtained through stronger baseline firing rL and/or longer burst durations (sacrifice of sparsity), without increasing FF.

Figure 9.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 9.

Interpretation of the empirical results in terms of the doubly stochastic model. A, Given the same mean firing rate rmean (≈ 20 sp/s in the schematics), less burstiness explains the observed low FF in empirical data (in correspondence to Fig. 1D). Scenario 1 with rL = 5 sp/s, rH = 100 sp/s, τL = 350 ms, and τH = 65 ms corresponds to the bursting regime with high FF. Low FF can be achieved either by smaller difference between rH and rL (Scenario 2 with with rL = 18 sp/s, rH = 30 sp/s, τL = 350 ms, and τH = 65 ms), and/or higher mean dwell time ratio τHτL (Scenario 3 with with rL = 5 sp/s, rH = 30 sp/s, τL = 350 ms, and τH = 525 ms). B, Heatmap of FF as a function of mean firing rate and the scaled difference of firing rates in the two states rH−rLrmean, with τL = 350 ms, τH = 65 ms, and time bin width Δ = 250 ms. Red line indicates the contour line of FF = 3. Region where FF < 3, which was mostly observed in empirical data, is characterized by smaller difference in rH and rL. Scenarios 1 and 2 in A are marked on the heatmap. C, Heatmap of FF as a function of mean firing rate and mean dwell time ratio τHτL, with τL = 350 ms, rL = 5 sp/s, and time bin width Δ = 250 ms. Red line indicates the contour line of FF = 3. Region where FF < 3, which was mostly observed in empirical data, is characterized by high τHτL. Scenarios 1 and 3 in A are marked on the heatmap. D, Heatmap of rmean as a function of rL and τHτL with FF = 1.5 fixed, which was typically observed in empirical data, and with τH = 65 ms, Δ = 250 ms. At fixed FF level, stronger mean firing rate under the preferred stimulus is achieved by elevating baseline firing and/or sacrificing spiking sparsity.

All these scenarios indicate that, with empirical data, the corresponding doubly stochastic model corresponds to the baseline L state being higher and/or the high-rate state being much longer than proposed for burst-coding. In other words, modulation of delay activity across stimulus conditions is not well accounted for by transient bursts, and is more consistent with nonbursting theories.

Discussion

In this study, we addressed an ongoing debate on the nature of WM delay activity between persistent activity and burst coding (Constantinidis et al., 2018; Lundqvist et al., 2018a), by examining implications of coding proposals on stability across trials of neuronal spiking activity during WM delays. Computational modeling of doubly stochastic Poisson processes and biophysically based spiking circuits demonstrated that FF is a sensitive measure that well dissociates these distinct coding regimes, and that burst-coding scenarios imply a substantial stimulus-tuned increase in FF during WM delays. In contrast, empirical spike-train data from PFC showed that FF is reduced during WM delays, with few neurons exhibiting positive correlation between FF and stimulus tuning of delay activity. We conclude that the contributions of burstiness to WM delay spiking activity in PFC are strongly constrained.

Relation to other empirical studies

We distinguish between across-trial stability and across-time stationarity as two important features to characterize WM delay activity, with across-trial stability as a key test of burst-coding models. Using FF to quantify trial-to-trial variability, time windows for analysis should be long enough relative to proposed characteristic burst durations. Of note, there are other measures of spiking variability which can better characterize spiking irregularity at faster timescales within trials (Compte et al., 2003; Shinomoto et al., 2009). Transient bursts in the context of proposals we are addressing are meant to be sharp and sparse, which could be distinct from bursting activity defined elsewhere, such as intrinsic burstiness of neurons (Compte et al., 2003; González-Burgos et al., 2019).

Observed reduction in FF at the population level during WM delays is consistent with prior studies of single-neuron spiking. FF is reduced following stimulus onset in multiple cortical regions (Nawrot et al., 2008; Churchland et al., 2010). In line with our findings, Chang et al. (2012) found that frontal eye field neurons had broad tuning of FF, which was dissociated from behavioral engagement and distinct from firing rates. Furthermore, Shafi et al. (2007) found reduction of trial-to-trial variability, quantified with coefficient of variation, during WM delays among delay-activated cells in lateral PFC. Our study contributes to this literature by further characterizing the relationship between FF and stimulus tuning of delay activity, and by relating empirical analysis to predictions for burst-coding proposals derived from statistical and circuit modeling.

Reevaluation of burst-coding proposals

Our analyses were based on burst-coding proposals from empirical studies (Lundqvist et al., 2016, 2018b) and on theoretical models for WM (Mongillo et al., 2008; Lundqvist et al., 2011; Mi et al., 2017). We showed that burst-coding proposals can require both the firing rate within burst states and FF to be very high in a biophysically constrained parameter space. Comparatively smaller FF can only be reached with stronger baseline firing and less sparseness (i.e., with less burstiness). Furthermore, regardless of FF magnitude, burst-coding models imply an increase in FF positively correlated with the mean firing rate during the delay, whereas opposite trends are displayed in empirical data. Our conclusion that intermittent transient bursts play a limited role in WM delay coding was established through reasoning by contradiction. We note that high FF does not necessarily imply bursting, as high FF can result from nonbursty overdispersion of interspike intervals or from slowly drifting excitability across trials.

Our findings are not inconsistent with previously observed intermittent bouts of elevated narrow-band LFP power (Lundqvist et al., 2016, 2018b). LFP signals do not reflect individual spike contributions or single-neuron dynamics (Kajikawa and Schroeder, 2011). Therefore, LFP bursting is not equivalent to bursting of neuronal spiking, and is not incompatible with persistent activity in neuronal spiking. We also note that, in the study of Lundqvist et al. (2016), LFP bursts were mostly observed during the stimulus presentation, and less so throughout the WM delay. It is thus unclear whether this LFP feature is directly related to WM maintenance, and whether those observations constitute evidence for burst-coding theoretical models of WM maintenance.

Implications for circuit models of WM

We simulated and compared three biophysically detailed circuit mechanisms for generating WM delay activity. A thorough characterization of each model in the high-dimensional parameter space is beyond the scope, but we showed that FF can well distinguish discharge patterns in distinct regimes. Attractor networks with recurrent connectivity have long been influential models of WM maintenance (Wang, 2001; Chaudhuri and Fiete, 2016). Wimmer et al. (2014) demonstrated that a diffusing bump attractor model captures key aspects of neuronal variability and correlations in spatial WM maintenance during the ODR task. There remain gaps between theoretical models and empirical recordings of neuronal spiking activity during WM. Empirical neuronal spiking data exhibit strong temporal variations in firing rates within trials, which fails to be captured by models with fixed-point attractors (Shafi et al., 2007; Druckmann and Chklovskii, 2012; Barak et al., 2013; Murray et al., 2017; Cueva et al., 2020). Furthermore, attractor network models do not necessarily capture levels of spike-train irregularity observed empirically (Renart et al., 2007; Roudi and Latham, 2007; Barbieri and Brunel, 2008).

Fixed-point attractor networks can reproduce reduced trial-to-trial variability during delay states, but with too low FF. Modeling studies found that this issue could be ameliorated in a regime of excitation-inhibition balance (Renart et al., 2007; Lim and Goldman, 2013, 2014; Hansel and Mato, 2013). We found that, in the model of Lim and Goldman (2013), the FF distribution resembled empirical data in magnitude, and reproduced the observed shrinkage of FF during delay compared with foreperiod better than the stable attractor model. Burst-coding models, by contrast, operate in a very different dynamic regime characterized by periodically occurring sharp bursts of spiking separated by longer periods of quiescence (Mongillo et al., 2008; Lundqvist et al., 2011; Fiebig and Lansner, 2017; Mi et al., 2017).

An important future direction is to directly incorporate further neurobiological data to inform and constrain WM circuit models. For instance, the chaotic attractor network model by Pereira and Brunel (2018) generated time-varying highly irregular activity within a WM attractor state. However, this model displayed strong firing-rate fluctuations and likely produces unrealistically high FF values. Our study provides a benchmark in terms of across-trial stability for future computational models of WM.

Limitations and future directions

One limitation of our datasets is that the number of trials per stimulus condition is on the order of 10, which can yield large variance in FF estimation. Future experiments can collect larger number of trials of single-neuron recordings in WM tasks with a long enough delay period. However, our main aim was to characterize the qualitative relationships between FF and bursting, mean firing rate, stimulus conditions, and WM delay activity, which was not dependent on high precision in estimation of FF. Furthermore, we simulated the empirically matched number of trials in the circuit models, which proved to be sufficient for this purpose.

Another limitation is that all three tasks require WM maintenance of a single item, whereas two- or three-item trials were analyzed by Lundqvist et al. (2016, 2018b). Our study therefore does not address the possibility that multiple-item WM utilizes a burst-dependent dynamic code, for example, in which short temporally separable bursts coding for each item enables multiplexing within WM circuits (Lundqvist et al., 2011; Mi et al., 2017). Future experimental and modeling studies should compare differences in PFC coding variability across different numbers of items in WM.

In addition to WM delays, foreperiod activity could exhibit characteristic variability and dynamics (Churchland et al., 2010; Murray et al., 2014). In particular, foreperiod activity could contain traces or reactivations of previous trials' stimulus information, which could elevate foreperiod FF and bias relative comparisons of delay versus foreperiod FF. Prior studies found that PFC neuronal coding of the previous trial's stimulus are typically weak and occur near the end of the foreperiod (Papadimitriou et al., 2017; Barbosa et al., 2020). Our results show that FF during delay is reduced relative to the foreperiod, which supports our inference about neuronal variability during maintenance of the more strongly encoded current trial's stimulus. Furthermore, model-derived predictions for the delay FF's magnitude and correlation with mean firing rates are unaffected by such foreperiod FF components.

We chose bottom-up generative modeling approaches, rather than directly fitting parameters of a statistical model to empirical data. Parametric statistical models rely on strong assumptions, and estimators can have large bias or variance. Developing statistical methods to capture and fit the complex discharge patterns on a single-trial basis is an important direction for future studies (Goris et al., 2014; Latimer et al., 2015; Charles et al., 2018), to more fully characterize the nature of spiking variability during WM.

Simultaneous recordings can inform the nature of shared neuronal variability (Crowe et al., 2010; Cohen and Kohn, 2011). The structure of neuronal correlations, and their modulation by task conditions, can constrain and inform models of cortical circuit dynamics and computation (Huang et al., 2019). Investigation of neuronal correlations (Leavitt et al., 2017b), in relation to coding subspaces (Murray et al., 2017), may provide crucial insights into the relationship between coordinated spiking dynamics and neural representations of cognitive states.

Footnotes

  • This work was supported by National Institutes of Health Grant R01 MH112746 to J.D.M.; National Science Foundation NeuroNex Grant 2015276 to J.D.M.; and National Institutes of Health, National Eye Institute Award R01 EY017077 to C.C. We thank Norman Lam for useful discussions.

  • The authors declare no competing financial interests.

  • Correspondence should be addressed to John D. Murray at john.murray{at}yale.edu

SfN exclusive license.

References

  1. ↵
    1. Amit DJ,
    2. Brunel N
    (1997) Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex. Cereb Cortex 7:237–252. doi:10.1093/cercor/7.3.237 pmid:9143444
    OpenUrlCrossRefPubMed
  2. ↵
    1. Barak O,
    2. Tsodyks M
    (2014) Working models of working memory. Curr Opin Neurobiol 25:20–24. doi:10.1016/j.conb.2013.10.008 pmid:24709596
    OpenUrlCrossRefPubMed
  3. ↵
    1. Barak O,
    2. Tsodyks M,
    3. Romo R
    (2010) Neuronal population coding of parametric working memory. J Neurosci 30:9424–9430. doi:10.1523/JNEUROSCI.1875-10.2010 pmid:20631171
    OpenUrlAbstract/FREE Full Text
  4. ↵
    1. Barak O,
    2. Sussillo D,
    3. Romo R,
    4. Tsodyks M,
    5. Abbott LF
    (2013) From fixed points to chaos: three models of delayed discrimination. Prog Neurobiol 103:214–222. doi:10.1016/j.pneurobio.2013.02.002 pmid:23438479
    OpenUrlCrossRefPubMed
  5. ↵
    1. Barbieri F,
    2. Brunel N
    (2008) Can attractor network models account for the statistics of firing during persistent activity in prefrontal cortex? Front Neurosci 2:114–122. doi:10.3389/neuro.01.003.2008 pmid:18982114
    OpenUrlCrossRefPubMed
  6. ↵
    1. Barbosa J,
    2. Stein H,
    3. Martinez RL,
    4. Galan-Gadea A,
    5. Li S,
    6. Dalmau J,
    7. Adam KC,
    8. Valls-Solé J,
    9. Constantinidis C,
    10. Compte A
    (2020) Interplay between persistent activity and activity-silent dynamics in the prefrontal cortex underlies serial biases in working memory. Nat Neurosci 23:1016–1024. doi:10.1038/s41593-020-0644-4 pmid:32572236
    OpenUrlCrossRefPubMed
  7. ↵
    1. Brody CD,
    2. Hernández A,
    3. Zainos A,
    4. Romo R
    (2003) Timing and neural encoding of somatosensory parametric working memory in macaque prefrontal cortex. Cereb Cortex 13:1196–1207. doi:10.1093/cercor/bhg100 pmid:14576211
    OpenUrlCrossRefPubMed
  8. ↵
    1. Brunel N,
    2. Wang XJ
    (2001) Effects of neuromodulation in a cortical network model of object working memory dominated by recurrent inhibition. J Comput Neurosci 11:63–85. doi:10.1023/A:1011204814320
    OpenUrlCrossRefPubMed
  9. ↵
    1. Cavanagh SE,
    2. Towers JP,
    3. Wallis JD,
    4. Hunt LT,
    5. Kennerley SW
    (2018) Reconciling persistent and dynamic hypotheses of working memory coding in prefrontal cortex. Nat Commun 9:3498. doi:10.1038/s41467-018-05873-3 pmid:30158519
    OpenUrlCrossRefPubMed
  10. ↵
    1. Chang MH,
    2. Armstrong KM,
    3. Moore T
    (2012) Dissociation of response variability from firing rate effects in frontal eye field neurons during visual stimulation, working memory, and attention. J Neurosci 32:2204–2216. doi:10.1523/JNEUROSCI.2967-11.2012 pmid:22323732
    OpenUrlAbstract/FREE Full Text
  11. ↵
    1. Charles AS,
    2. Park M,
    3. Weller JP,
    4. Horwitz GD,
    5. Pillow JW
    (2018) Dethroning the Fano factor: a flexible, model-based approach to partitioning neural variability. Neural Comput 30:1012–1045. doi:10.1162/neco_a_01062 pmid:29381442
    OpenUrlCrossRefPubMed
  12. ↵
    1. Chaudhuri R,
    2. Fiete I
    (2016) Computational principles of memory. Nat Neurosci 19:394–403. doi:10.1038/nn.4237 pmid:26906506
    OpenUrlCrossRefPubMed
  13. ↵
    1. Churchland AK,
    2. Kiani R,
    3. Chaudhuri R,
    4. Wang XJ,
    5. Pouget A,
    6. Shadlen MN
    (2011) Variance as a signature of neural computations during decision making. Neuron 69:818–831. doi:10.1016/j.neuron.2010.12.037 pmid:21338889
    OpenUrlCrossRefPubMed
  14. ↵
    1. Churchland MM,
    2. Abbott LF
    (2012) Two layers of neural variability. Nat Neurosci 15:1472–1474. doi:10.1038/nn.3247 pmid:23103992
    OpenUrlCrossRefPubMed
  15. ↵
    1. Churchland MM,
    2. Yu BM,
    3. Cunningham JP,
    4. Sugrue LP,
    5. Cohen MR,
    6. Corrado GS,
    7. Newsome WT,
    8. Clark AM,
    9. Hosseini P,
    10. Scott BB,
    11. Bradley DC,
    12. Smith MA,
    13. Kohn A,
    14. Movshon JA,
    15. Armstrong KM,
    16. Moore T,
    17. Chang SW,
    18. Snyder LH,
    19. Lisberger SG,
    20. Priebe NJ, et al
    . (2010) Stimulus onset quenches neural variability: a widespread cortical phenomenon. Nat Neurosci 13:369–378. doi:10.1038/nn.2501 pmid:20173745
    OpenUrlCrossRefPubMed
  16. ↵
    1. Cohen MR,
    2. Kohn A
    (2011) Measuring and interpreting neuronal correlations. Nat Neurosci 14:811–819. doi:10.1038/nn.2842 pmid:21709677
    OpenUrlCrossRefPubMed
  17. ↵
    1. Compte A,
    2. Brunel N,
    3. Goldman-Rakic PS,
    4. Wang XJ
    (2000) Synaptic mechanisms and network dynamics underlying spatial working memory in a cortical network model. Cereb Cortex 10:910–923. doi:10.1093/cercor/10.9.910 pmid:10982751
    OpenUrlCrossRefPubMed
  18. ↵
    1. Compte A,
    2. Constantinidis C,
    3. Tegner J,
    4. Raghavachari S,
    5. Chafee MV,
    6. Goldman-Rakic PS,
    7. Wang XJ
    (2003) Temporally irregular mnemonic persistent activity in prefrontal neurons of monkeys during a delayed response task. J Neurophysiol 90:3441–3454. doi:10.1152/jn.00949.2002
    OpenUrlCrossRefPubMed
  19. ↵
    1. Constantinidis C,
    2. Franowicz MN,
    3. Goldman-Rakic PS
    (2001) Coding specificity in cortical microcircuits: a multiple-electrode analysis of primate prefrontal cortex. J Neurosci 21:3646–3655. doi:10.1523/JNEUROSCI.21-10-03646.2001
    OpenUrlAbstract/FREE Full Text
  20. ↵
    1. Constantinidis C,
    2. Funahashi S,
    3. Lee D,
    4. Murray JD,
    5. Qi XL,
    6. Wang M,
    7. Arnsten AF
    (2018) Persistent spiking activity underlies working memory. J Neurosci 38:7020–7028. doi:10.1523/JNEUROSCI.2486-17.2018 pmid:30089641
    OpenUrlAbstract/FREE Full Text
  21. ↵
    1. Cox DR
    (1955) Some statistical methods connected with series of events. J R Stat Soc B 17:129–157. doi:10.1111/j.2517-6161.1955.tb00188.x
    OpenUrlCrossRef
  22. ↵
    1. Crowe DA,
    2. Averbeck BB,
    3. Chafee MV
    (2010) Rapid sequences of population activity patterns dynamically encode task-critical spatial information in parietal cortex. J Neurosci 30:11640–11653. doi:10.1523/JNEUROSCI.0954-10.2010 pmid:20810885
    OpenUrlAbstract/FREE Full Text
  23. ↵
    1. Cueva CJ,
    2. Saez A,
    3. Marcos E,
    4. Genovesio A,
    5. Jazayeri M,
    6. Romo R,
    7. Salzman CD,
    8. Shadlen MN,
    9. Fusi S
    (2020) Low-dimensional dynamics for working memory and time encoding. Proc Natl Acad Sci USA 117:23021–23032. doi:10.1073/pnas.1915984117 pmid:32859756
    OpenUrlAbstract/FREE Full Text
  24. ↵
    1. Deco G,
    2. Hugues E
    (2012) Neural network mechanisms underlying stimulus driven variability reduction. PLoS Comput Biol 8:e1002395. doi:10.1371/journal.pcbi.1002395 pmid:22479168
    OpenUrlCrossRefPubMed
  25. ↵
    1. Druckmann S,
    2. Chklovskii DB
    (2012) Neuronal circuits underlying persistent representations despite time varying activity. Curr Biol 22:2095–2103. doi:10.1016/j.cub.2012.08.058 pmid:23084992
    OpenUrlCrossRefPubMed
  26. ↵
    1. Durstewitz D,
    2. Gabriel T
    (2007) Dynamical basis of irregular spiking in NMDA-driven prefrontal cortex neurons. Cereb Cortex 17:894–908. doi:10.1093/cercor/bhk044 pmid:16740581
    OpenUrlCrossRefPubMed
  27. ↵
    1. Fano U
    (1947) Ionization yield of radiations: II. The fluctuations of the number of ions. Phys Rev 72:26–29. doi:10.1103/PhysRev.72.26
    OpenUrlCrossRef
  28. ↵
    1. Fiebig F,
    2. Lansner A
    (2017) A spiking working memory model based on Hebbian short-term potentiation. J Neurosci 37:83–96. doi:10.1523/JNEUROSCI.1989-16.2016 pmid:28053032
    OpenUrlAbstract/FREE Full Text
  29. ↵
    1. Funahashi S,
    2. Bruce CJ,
    3. Goldman-Rakic PS
    (1989) Mnemonic coding of visual space in the monkey's dorsolateral prefrontal cortex. J Neurophysiol 61:331–349. doi:10.1152/jn.1989.61.2.331 pmid:2918358
    OpenUrlCrossRefPubMed
  30. ↵
    1. Funahashi S,
    2. Chafee MV,
    3. Goldman-Rakic PS
    (1993) Prefrontal neuronal activity in rhesus monkeys performing a delayed anti-saccade task. Nature 365:753–756. doi:10.1038/365753a0 pmid:8413653
    OpenUrlCrossRefPubMed
  31. ↵
    1. Goldman MS,
    2. Compte A,
    3. Wang XJ
    (2008) Neural integrator models. In: Encyclopedia of neuroscience, pp 6165–6178. San Diego: Academic.
  32. ↵
    1. González-Burgos G,
    2. Miyamae T,
    3. Krimer Y,
    4. Gulchina Y,
    5. Pafundo DE,
    6. Krimer O,
    7. Bazmi H,
    8. Arion D,
    9. Enwright JF,
    10. Fish KN,
    11. Lewis DA
    (2019) Distinct properties of layer 3 pyramidal neurons from prefrontal and parietal areas of the monkey neocortex. J Neurosci 39:7277–7290. doi:10.1523/JNEUROSCI.1210-19.2019 pmid:31341029
    OpenUrlAbstract/FREE Full Text
  33. ↵
    1. Goris RL,
    2. Movshon JA,
    3. Simoncelli EP
    (2014) Partitioning neuronal variability. Nat Neurosci 17:858–865. doi:10.1038/nn.3711 pmid:24777419
    OpenUrlCrossRefPubMed
  34. ↵
    1. Hansel D,
    2. Mato G
    (2013) Short-term plasticity explains irregular persistent activity in working memory tasks. J Neurosci 33:133–149. doi:10.1523/JNEUROSCI.3455-12.2013 pmid:23283328
    OpenUrlAbstract/FREE Full Text
  35. ↵
    1. Huang C,
    2. Ruff DA,
    3. Pyle R,
    4. Rosenbaum R,
    5. Cohen MR,
    6. Doiron B
    (2019) Circuit models of low-dimensional shared variability in cortical networks. Neuron 101:337–348.e4. doi:10.1016/j.neuron.2018.11.034 pmid:30581012
    OpenUrlCrossRefPubMed
  36. ↵
    1. Hussar C,
    2. Pasternak T
    (2010) Trial-to-trial variability of the prefrontal neurons reveals the nature of their engagement in a motion discrimination task. Proc Natl Acad Sci USA 107:21842–21847. doi:10.1073/pnas.1009956107 pmid:21098286
    OpenUrlAbstract/FREE Full Text
  37. ↵
    1. Kajikawa Y,
    2. Schroeder CE
    (2011) How local is the local field potential? Neuron 72:847–858. doi:10.1016/j.neuron.2011.09.029
    OpenUrlCrossRefPubMed
  38. ↵
    1. Latimer KW,
    2. Yates JL,
    3. Meister ML,
    4. Huk AC,
    5. Pillow JW
    (2015) Single-trial spike trains in parietal cortex reveal discrete steps during decision-making. Science 349:184–187. doi:10.1126/science.aaa4056 pmid:26160947
    OpenUrlAbstract/FREE Full Text
  39. ↵
    1. Leavitt ML,
    2. Mendoza-Halliday D,
    3. Martinez-Trujillo JC
    (2017a) Sustained activity encoding working memories: not fully distributed. Trends Neurosci 40:328–346. doi:10.1016/j.tins.2017.04.004 pmid:28515011
    OpenUrlCrossRefPubMed
  40. ↵
    1. Leavitt ML,
    2. Pieper F,
    3. Sachs AJ,
    4. Martinez-Trujillo JC
    (2017b) Correlated variability modifies working memory fidelity in primate prefrontal neuronal ensembles. Proc Natl Acad Sci USA 114:E2494–E2503. doi:10.1073/pnas.1619949114 pmid:28275096
    OpenUrlAbstract/FREE Full Text
  41. ↵
    1. Lim S,
    2. Goldman MS
    (2013) Balanced cortical microcircuitry for maintaining information in working memory. Nat Neurosci 16:1306–1314. doi:10.1038/nn.3492 pmid:23955560
    OpenUrlCrossRefPubMed
  42. ↵
    1. Lim S,
    2. Goldman MS
    (2014) Balanced cortical microcircuitry for spatial working memory based on corrective feedback control. J Neurosci 34:6790–6806. doi:10.1523/JNEUROSCI.4602-13.2014 pmid:24828633
    OpenUrlAbstract/FREE Full Text
  43. ↵
    1. Litwin-Kumar A,
    2. Doiron B
    (2012) Slow dynamics and high variability in balanced cortical networks with clustered connections. Nat Neurosci 15:1498–1505. doi:10.1038/nn.3220 pmid:23001062
    OpenUrlCrossRefPubMed
  44. ↵
    1. Lundqvist M,
    2. Herman P,
    3. Lansner A
    (2011) Theta and gamma power increases and alpha/beta power decreases with memory load in an attractor network model. J Cogn Neurosci 23:3008–3020. doi:10.1162/jocn_a_00029 pmid:21452933
    OpenUrlCrossRefPubMed
  45. ↵
    1. Lundqvist M,
    2. Rose J,
    3. Herman P,
    4. Brincat SL,
    5. Buschman TJ,
    6. Miller EK
    (2016) Gamma and beta bursts underlie working memory. Neuron 90:152–164. doi:10.1016/j.neuron.2016.02.028 pmid:26996084
    OpenUrlCrossRefPubMed
  46. ↵
    1. Lundqvist M,
    2. Herman P,
    3. Miller EK
    (2018a) Working memory: delay activity, yes! Persistent activity? Maybe not. J Neurosci 38:7013–7019. doi:10.1523/JNEUROSCI.2485-17.2018 pmid:30089640
    OpenUrlAbstract/FREE Full Text
  47. ↵
    1. Lundqvist M,
    2. Herman P,
    3. Warden MR,
    4. Brincat SL,
    5. Miller EK
    (2018b) Gamma and beta bursts during working memory readout suggest roles in its volitional control. Nat Commun 9:394. doi:10.1038/s41467-017-02791-8 pmid:29374153
    OpenUrlCrossRefPubMed
  48. ↵
    1. Markowitz DA,
    2. Curtis CE,
    3. Pesaran B
    (2015) Multiple component networks support working memory in prefrontal cortex. Proc Natl Acad Sci USA 112:11084–11089. doi:10.1073/pnas.1504172112 pmid:26283366
    OpenUrlAbstract/FREE Full Text
  49. ↵
    1. Mi Y,
    2. Katkov M,
    3. Tsodyks M
    (2017) Synaptic correlates of working memory capacity. Neuron 93:323–330. doi:10.1016/j.neuron.2016.12.004 pmid:28041884
    OpenUrlCrossRefPubMed
  50. ↵
    1. Miller EK,
    2. Erickson CA,
    3. Desimone R
    (1996) Neural mechanisms of visual working memory in prefrontal cortex of the macaque. J Neurosci 16:5154–5167. pmid:8756444
    OpenUrlAbstract/FREE Full Text
  51. ↵
    1. Miller P
    (2006) Analysis of spike statistics in neuronal systems with continuous attractors or multiple, discrete attractor states. Neural Comput 18:1268–1317. doi:10.1162/neco.2006.18.6.1268 pmid:16764505
    OpenUrlCrossRefPubMed
  52. ↵
    1. Mongillo G,
    2. Barak O,
    3. Tsodyks M
    (2008) Synaptic theory of working memory. Science 319:1543–1546. doi:10.1126/science.1150769 pmid:18339943
    OpenUrlAbstract/FREE Full Text
  53. ↵
    1. Murray JD,
    2. Bernacchia A,
    3. Freedman DJ,
    4. Romo R,
    5. Wallis JD,
    6. Cai X,
    7. Padoa-Schioppa C,
    8. Pasternak T,
    9. Seo H,
    10. Lee D,
    11. Wang XJ
    (2014) A hierarchy of intrinsic timescales across primate cortex. Nat Neurosci 17:1661–1663. doi:10.1038/nn.3862 pmid:25383900
    OpenUrlCrossRefPubMed
  54. ↵
    1. Murray JD,
    2. Bernacchia A,
    3. Roy NA,
    4. Constantinidis C,
    5. Romo R,
    6. Wang XJ
    (2017) Stable population coding for working memory coexists with heterogeneous neural dynamics in prefrontal cortex. Proc Natl Acad Sci USA 114:394–399. doi:10.1073/pnas.1619449114 pmid:28028221
    OpenUrlAbstract/FREE Full Text
  55. ↵
    1. Nawrot MP,
    2. Boucsein C,
    3. Rodriguez Molina V,
    4. Riehle A,
    5. Aertsen A,
    6. Rotter S
    (2008) Measurement of variability dynamics in cortical spike trains. J Neurosci Methods 169:374–390. doi:10.1016/j.jneumeth.2007.10.013 pmid:18155774
    OpenUrlCrossRefPubMed
  56. ↵
    1. Oram MW,
    2. Wiener MC,
    3. Lestienne R,
    4. Richmond BJ
    (1999) Stochastic nature of precisely timed spike patterns in visual system neuronal responses. J Neurophysiol 81:3021–3033. doi:10.1152/jn.1999.81.6.3021 pmid:10368417
    OpenUrlCrossRefPubMed
  57. ↵
    1. Papadimitriou C,
    2. White RL,
    3. Snyder LH
    (2017) Ghosts in the machine: II. Neural correlates of memory interference from the previous trial. Cereb Cortex 27:2513–2527. pmid:27114176
    OpenUrlCrossRefPubMed
  58. ↵
    1. Pereira U,
    2. Brunel N
    (2018) Attractor dynamics in networks with learning rules inferred from in vivo data. Neuron 99:227–238.e4. doi:10.1016/j.neuron.2018.05.038 pmid:29909997
    OpenUrlCrossRefPubMed
  59. ↵
    1. Purcell BA,
    2. Heitz RP,
    3. Cohen JY,
    4. Schall JD
    (2012) Response variability of frontal eye field neurons modulates with sensory input and saccade preparation but not visual search salience. J Neurophysiol 108:2737–2750. doi:10.1152/jn.00613.2012 pmid:22956785
    OpenUrlCrossRefPubMed
  60. ↵
    1. Renart A,
    2. Moreno-Bote R,
    3. Wang XJ,
    4. Parga N
    (2007) Mean-driven and fluctuation-driven persistent activity in recurrent networks. Neural Comput 19:1–46. doi:10.1162/neco.2007.19.1.1 pmid:17134316
    OpenUrlCrossRefPubMed
  61. ↵
    1. Riley MR,
    2. Qi XL,
    3. Zhou X,
    4. Constantinidis C
    (2018) Anterior-posterior gradient of plasticity in primate prefrontal cortex. Nat Commun 9:3790. doi:10.1038/s41467-018-06226-w pmid:30224705
    OpenUrlCrossRefPubMed
  62. ↵
    1. Romo R,
    2. Brody CD,
    3. Hernández A,
    4. Lemus L
    (1999) Neuronal correlates of parametric working memory in the prefrontal cortex. Nature 399:470–473. doi:10.1038/20939 pmid:10365959
    OpenUrlCrossRefPubMed
  63. ↵
    1. Roudi Y,
    2. Latham PE
    (2007) A balanced memory network. PLoS Comput Biol 3:1679–1700. doi:10.1371/journal.pcbi.0030141 pmid:17845070
    OpenUrlCrossRefPubMed
  64. ↵
    1. Shadlen MN,
    2. Newsome WT
    (1998) The variable discharge of cortical neurons: implications for connectivity, computation, and information coding. J Neurosci 18:3870–3896. pmid:9570816
    OpenUrlAbstract/FREE Full Text
  65. ↵
    1. Shafi M,
    2. Zhou Y,
    3. Quintana J,
    4. Chow C,
    5. Fuster J,
    6. Bodner M
    (2007) Variability in neuronal activity in primate cortex during working memory tasks. Neuroscience 146:1082–1108. doi:10.1016/j.neuroscience.2006.12.072 pmid:17418956
    OpenUrlCrossRefPubMed
  66. ↵
    1. Shinomoto S,
    2. Kim H,
    3. Shimokawa T,
    4. Matsuno N,
    5. Funahashi S,
    6. Shima K,
    7. Fujita I,
    8. Tamura H,
    9. Doi T,
    10. Kawano K,
    11. Inaba N,
    12. Fukushima K,
    13. Kurkin S,
    14. Kurata K,
    15. Taira M,
    16. Tsutsui KI,
    17. Komatsu H,
    18. Ogawa T,
    19. Koida K,
    20. Tanji J, et al
    . (2009) Relating neuronal firing patterns to functional differentiation of cerebral cortex. PLoS Comput Biol 5:e1000433. doi:10.1371/journal.pcbi.1000433 pmid:19593378
    OpenUrlCrossRefPubMed
  67. ↵
    1. Softky WR,
    2. Koch C
    (1993) The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPs. J Neurosci 13:334–350. doi:10.1523/JNEUROSCI.13-01-00334.1993
    OpenUrlAbstract/FREE Full Text
  68. ↵
    1. Spaak E,
    2. Watanabe K,
    3. Funahashi S,
    4. Stokes MG
    (2017) Stable and dynamic coding for working memory in primate prefrontal cortex. J Neurosci 37:6503–6516. doi:10.1523/JNEUROSCI.3364-16.2017 pmid:28559375
    OpenUrlAbstract/FREE Full Text
  69. ↵
    1. Tiganj Z,
    2. Cromer JA,
    3. Roy JE,
    4. Miller EK,
    5. Howard MW
    (2018) Compressed timeline of recent experience in monkey lateral prefrontal cortex. J Cogn Neurosci 30:935–950. doi:10.1162/jocn_a_01273 pmid:29698121
    OpenUrlCrossRefPubMed
  70. ↵
    1. Tolhurst D,
    2. Movshon JA,
    3. Thompson I
    (1981) The dependence of response amplitude and variance of cat visual cortical neurones on stimulus contrast. Exp Brain Res 41:414–419. doi:10.1007/BF00238900 pmid:7215502
    OpenUrlCrossRefPubMed
  71. ↵
    1. Tolhurst DJ,
    2. Movshon JA,
    3. Dean AF
    (1983) The statistical reliability of signals in single neurons in cat and monkey visual cortex. Vision Res 23:775–785. doi:10.1016/0042-6989(83)90200-6
    OpenUrlCrossRefPubMed
  72. ↵
    1. Wang XJ
    (1999) Synaptic basis of cortical persistent activity: the importance of NMDARs to working memory. J Neurosci 19:9587–9603. doi:10.1523/JNEUROSCI.19-21-09587.1999
    OpenUrlAbstract/FREE Full Text
  73. ↵
    1. Wang XJ
    (2001) Synaptic reverberation underlying mnemonic persistent activity. Trends Neurosci 24:455–463. doi:10.1016/s0166-2236(00)01868-3 pmid:11476885
    OpenUrlCrossRefPubMed
  74. ↵
    1. Wang XJ
    (2002) Probabilistic decision making by slow reverberation in cortical circuits. Neuron 36:955–968. doi:10.1016/s0896-6273(02)01092-9 pmid:12467598
    OpenUrlCrossRefPubMed
  75. ↵
    1. Wimmer K,
    2. Nykamp DQ,
    3. Constantinidis C,
    4. Compte A
    (2014) Bump attractor dynamics in prefrontal cortex explains behavioral precision in spatial working memory. Nat Neurosci 17:431–439. doi:10.1038/nn.3645 pmid:24487232
    OpenUrlCrossRefPubMed
Back to top

In this issue

The Journal of Neuroscience: 41 (43)
Journal of Neuroscience
Vol. 41, Issue 43
27 Oct 2021
  • Table of Contents
  • Table of Contents (PDF)
  • About the Cover
  • Index by author
  • Ed Board (PDF)
Email

Thank you for sharing this Journal of Neuroscience article.

NOTE: We request your email address only to inform the recipient that it was you who recommended this article, and that it is not junk mail. We do not retain these email addresses.

Enter multiple addresses on separate lines or separate them with commas.
Trial-to-Trial Variability of Spiking Delay Activity in Prefrontal Cortex Constrains Burst-Coding Models of Working Memory
(Your Name) has forwarded a page to you from Journal of Neuroscience
(Your Name) thought you would be interested in this article in Journal of Neuroscience.
CAPTCHA
This question is for testing whether or not you are a human visitor and to prevent automated spam submissions.
Print
View Full Page PDF
Citation Tools
Trial-to-Trial Variability of Spiking Delay Activity in Prefrontal Cortex Constrains Burst-Coding Models of Working Memory
Daming Li, Christos Constantinidis, John D. Murray
Journal of Neuroscience 27 October 2021, 41 (43) 8928-8945; DOI: 10.1523/JNEUROSCI.0167-21.2021

Citation Manager Formats

  • BibTeX
  • Bookends
  • EasyBib
  • EndNote (tagged)
  • EndNote 8 (xml)
  • Medlars
  • Mendeley
  • Papers
  • RefWorks Tagged
  • Ref Manager
  • RIS
  • Zotero
Respond to this article
Request Permissions
Share
Trial-to-Trial Variability of Spiking Delay Activity in Prefrontal Cortex Constrains Burst-Coding Models of Working Memory
Daming Li, Christos Constantinidis, John D. Murray
Journal of Neuroscience 27 October 2021, 41 (43) 8928-8945; DOI: 10.1523/JNEUROSCI.0167-21.2021
del.icio.us logo Digg logo Reddit logo Twitter logo CiteULike logo Facebook logo Google logo Mendeley logo
  • Tweet Widget
  • Facebook Like
  • Google Plus One

Jump to section

  • Article
    • Abstract
    • Introduction
    • Materials and Methods
    • Results
    • Discussion
    • Footnotes
    • References
  • Figures & Data
  • Info & Metrics
  • eLetters
  • PDF

Keywords

  • bursting
  • computational model
  • Fano factor
  • prefrontal cortex
  • spike-train statistics
  • working memory

Responses to this article

Respond to this article

Jump to comment:

No eLetters have been published for this article.

Related Articles

Cited By...

More in this TOC Section

Research Articles

  • Identification of the Acid-Sensitive Site Critical for Chloral Hydrate (CH) Activation of the Proton-Activated Chloride Channel
  • Disentangling Object Category Representations Driven by Dynamic and Static Visual Input
  • Irrelevant Threats Linger and Affect Behavior in High Anxiety
Show more Research Articles

Behavioral/Cognitive

  • Disentangling Object Category Representations Driven by Dynamic and Static Visual Input
  • Irrelevant Threats Linger and Affect Behavior in High Anxiety
  • Multisession Anodal Transcranial Direct Current Stimulation Enhances Adult Hippocampal Neurogenesis and Context Discrimination in Mice
Show more Behavioral/Cognitive
  • Home
  • Alerts
  • Visit Society for Neuroscience on Facebook
  • Follow Society for Neuroscience on Twitter
  • Follow Society for Neuroscience on LinkedIn
  • Visit Society for Neuroscience on Youtube
  • Follow our RSS feeds

Content

  • Early Release
  • Current Issue
  • Issue Archive
  • Collections

Information

  • For Authors
  • For Advertisers
  • For the Media
  • For Subscribers

About

  • About the Journal
  • Editorial Board
  • Privacy Policy
  • Contact
(JNeurosci logo)
(SfN logo)

Copyright © 2023 by the Society for Neuroscience.
JNeurosci Online ISSN: 1529-2401

The ideas and opinions expressed in JNeurosci do not necessarily reflect those of SfN or the JNeurosci Editorial Board. Publication of an advertisement or other product mention in JNeurosci should not be construed as an endorsement of the manufacturer’s claims. SfN does not assume any responsibility for any injury and/or damage to persons or property arising from or related to any use of any material contained in JNeurosci.