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.
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:
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
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:
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:
Hence the variance of the spike count is given by the following:
The autocovariance of the underlying telegraph firing rate process has the following form:
Therefore, the characteristic time of the exponential decay is
Therefore, FF is given by the following:
FF therefore depends on the four parameters
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
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:
We used the following values for the recurrent synaptic conductances (in nS) in the N = 2000 neurons network: for pyramidal cells,
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,
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:
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.
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.
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
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
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
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
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
If WM were subserved by sharp intermittent bursts, we would have rH ≫ rL and
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
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:
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
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.
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.
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).
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.
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).
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(
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.
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
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
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
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
