| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
Volen Center for Complex Systems, Brandeis University, Waltham, Massachusetts
Submitted 11 May 2005; accepted in final form 14 October 2005
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIXGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIXGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
; Compte et al. 2003; Shadlen and Newsome 1994 Softky and Koch 1993). Fluctuations in spike discharges are often correlated between simultaneously recorded cells on a trial-by-trial basis (Bair et al. 2001; Gawne and Richmond 1993; Lee et al. 1998; Zohary et al. 1994); this covariance greatly impacts the extent to which noise can be averaged out over a neuronal pool and so determines the accuracy with which a stimulus feature can be extracted from a large population of neurons (Abbott and Dayan 1999; Averbeck and Lee 2004; Shadlen and Newsome 1996; Shamir and Sompolinsky 2004; Sompolinsky et al. 2002). Correlations arise from common inputs caused by both overlapping afferents and local synaptic interconnections. For example, cells in the primary visual cortex become correlated because of shared thalamic inputs as well as lateral connections within the cortex. Exactly how correlated noise enhances or decreases the coding efficiency for an oriented visual stimulus critically depends on the (feedforward vs. recurrent) network architecture (Seriès et al. 2004).
Correlated fluctuations provide a valuable probe into the dynamic nature of a neural network. For responses of sensory cells to external stimuli, the effects of direct sensory inputs and local connectivity are confounded. In contrast, neurons in working memory systems exhibit persistent activity in the absence of external stimulation, when the subject must hold information transiently in the mind. Hence fluctuations of persistent activity (Compte et al. 2003; Constantinidis and Goldman-Rakic 2002) are believed to result from the intrinsic dynamics of a working memory microcircuit. The cellular mechanisms of mnemonic persistent activity represent a topic of intense current interest (Brody et al. 2003; Major and Tank 2004; Wang 2001). Of special interest are working memory circuits that encode an analog quantity, such as those integrating a velocity signal into a persistent change in position or direction (Aksay et al. 2000; Taube and Bassett 2003) or storing in active short-term memory the spatial location or temporal frequency of a sensory stimulus (Funahashi et al. 1989; Romo et al. 1999). After a transient input, neurons in such a network exhibit sustained changes of activity that are tuned with a stimulus feature in a graded manner. This suggests the system could possess a continuous family of stable self-sustained activity states. Neuronal tuning curves that arise from such activity states are often bell-shaped functions or monotonic functions of the encoded quantity.
For a parametric working memory system, if the delay firing rates of two monotonically tuned neurons are plotted against each other, a curved line results, as shown in Fig. 1A. Such a curved line is referred to as a line attractor (Seung 1996). The realization of a continuum of monotonically tuned persistent firing states requires fine-tuning of network parameters (Miller et al. 2003; Seung et al. 2000b). In an alternative scenario, the network could display multiple persistent states, not truly a continuum, by virtue of a set of discrete stable activity patterns that can be robustly realized without fine-tuning of parameters (Goldman et al. 2003; Koulakov et al. 2002). Different stimuli cause the network to change discontinuously from one pattern of activity to another, whereas background noise does not cause a switch in the state of persistent activity. Hence the discrete network model has a limited sensitivity to small differences in stimuli, but is robust to noise in the circuit. Whether the continuous or discrete model better describes neural integrators remains an open question (Major and Tank 2004).
|
), we reasoned that unusual fluctuation properties should be a conspicuous feature of continuous "parametric" memory networks, in which mnemonic neural firing is monotonically tuned to an analog quantity (Aksay et al. 2000; Miller et al. 2003; Romo et al. 1999; Seung 1996; Seung et al. 2000b). A more robust discrete integrator network (Goldman et al. 2003; Koulakov et al. 2002) would not show the same random drifts in response to noise, unless the noise were strong enough to cause transitions between the multiple discrete states (Miller 2006). If the scale of noise fluctuations sets the coarseness of the system, a system is considered as quasi-continuous whenever the difference in discrete rates is less than the noise-induced variation in rate during a trial. The purpose of this study is to analyze temporal fluctuations in firing patterns of working memory models with monotonically tuned persistent states. We show that power-law fluctuations are a salient feature of continuous attractor networks.
|
METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIXGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
We developed a cortical microcircuit model for the task of parametric working memory, as published previously (Miller et al. 2003). The task requires macaque monkeys to encode and maintain in memory a vibrational frequency, before comparison with a second vibrational frequency. In the secondary somatosensory cortex of a monkey, the spiking response of a neuron to the vibrotactile stimuli varies either positively monotonically or negatively monotonically with the stimulus frequency. Neurons in the prefrontal and premotor cortices, which receive inputs from the secondary somatosensory cortex, show tuned mnemonic activity that can be sustained for up to 6 s during a delay after the end of the initial stimulus (before a 2nd comparison stimulus). Our model network is intended to reflect the activity of prefrontal cortical neurons. We included two sets of populations, corresponding to sets of neurons that receive inputs separately from the secondary somatosensory cortex. One set of inputs is positively monotonically tuned while the other set is negatively monotonically tuned to the vibrational stimulus frequency.
We used interconnected integrate-and-fire neurons, grouped into populations that had strong intrapopulation connections and weaker connections between populations. The connection strengths decrease exponentially with the difference between population numbers (Fig. 1B). The decay in connection strength was more gradual from high-to-low population numbers, so that the neurons in populations with low index received more excitatory current from other cells than those with higher index. This led to a range of thresholds for responses to external input, with the low-index populations being the most excitable and having the highest spontaneous rates.
We couple the two sets of oppositely tuned populations by cross-inhibition. The positively monotonic populations are labeled 1 for the most excitable to 12 for the least excitable, and the negatively monotonic neurons are labeled 1* for the most excitable to 12* for the least excitable. Cross-inhibition is not uniform, but strongest between populations labeled by i and (14 – i)*, where i runs from 2 to 12. Hence the increase in activity of a highly excitable positively monotonic population, after a weak stimulus, helps to reduce the activity of a less excitable negatively monotonic population. A larger stimulus results in a strong increase in activity for the less excitable positively monotonic populations and a concurrent decrease in activity of the more excitable negatively monotonic populations.
We realized a continuous attractor, with moderate recurrent excitation, where the total synaptic weight from all other excitatory populations is of the same order as the self-excitation within a population. For each population, the background external input amplitude and intrapopulation recurrent excitation strength were adjusted to a particular combination, corresponding to a point called a "cusp" in the two-dimensional parameter space, indicated by C in Fig. 2A. For a strength of recurrent excitation greater than at the cusp, bistability can exist between distinct up (persistently active) and down (resting) states, with a gap between them (A in Fig. 2A). For weaker recurrent excitation than at the cusp, the population firing rate varies smoothly with the excitatory drive (B in Fig. 2A). At the cusp, there is a possibility of an approximately vertical line segment, indicating a range of stable firing rates for the population with a fixed background excitatory input (position C, Fig. 2A, right).
|
) by creating many bistable populations with a range of thresholds. We achieved this by increasing the recurrent connection strength within each population while reducing the connection strengths between populations (see point A in Fig. 2A) and increasing the average leak conductances. The bistable populations each converge onto a readout population of 400 neurons, whose firing rates encode in a step-like manner the number of active bistable populations. We use the statistical properties of the readout cells for a fair comparison between the discrete model and the continuous model. The set of bistable populations should just be considered as one particular mechanism for generating the more realistic properties of the readout cells, which are suitable for experimental comparison. We chose stimulus strengths for the discrete attractor, so that each stimulus would cause the discrete system to be in the same state for all trials with that stimulus. In this way, we can use the discrete system as a control to show that typical networks with a single stable state after a stimulus do not exhibit the same type of fluctuations seen in the continuous attractor network.
Single neuron parameters
Individual neurons are simulated using the single-compartment leaky-integrate-and-fire model (Tuckwell 1988), such that the membrane potential, Vi, of cell, i, follows the current-balance equation
 | (1) |
ref, before continuing to follow Eq. 1.
The total synaptic drive for excitation or inhibition (SE or SI) is the sum of synaptic inputs from all presynaptic neurons j
 | (2) |
i is the relative synaptic weight from cell j to cell i, and sj is the synaptic current gating variable activated by the presynaptic neuron j firing spikes at times tspike,j. Specifically, for excitatory synapses, we have
 | (3) |
 | (4) |
s. The probability of vesicular release, 
PR(t)
, is described in the next subsection.
Background noisy input to all neurons is simulated using uncorrelated Poisson spike trains at a rate, rext, through nonsaturating synapses, of conductance gext, which are gated according to
 | (5) |
ext following spikes at times, tspike,ext.
Similarly, during the stimulus, Poisson spike trains of rate, 
, generate additional excitation through 
-amino-3-hydroxy-5-methylisoxazole-4-propionic acid receptor (AMPAR)-mediated synapses of conductance, gcue, multiplied by a gating variable, scue, which follows
 | (6) |
, represents the combined spike rate of multiple afferents from secondary somatosensory cortex, and increases or decreases linearly with stimulus frequency for the positively monotonic or negatively monotonic populations, respectively.
In the network models presented here, background and stimulus inputs are mediated by -amino-3-hydroxy-5-methylisoxazole-4-propionic acid (AMPA) receptors, with ext = 2ms, recurrent excitation through N-methyl-D-aspartate (NMDA) receptors (Wang 1999, 2001) with s = 100 ms and VE = 0 mV, and inhibition through GABAA receptors with s = 10 ms and VI = –70 mV. In the continuous attractor network, the cellular parameters are as follows for excitatory cells: Cm = 0.5 nF, gL = 38.4 nS, VL = –70 mV,Vreset = –60 mV, Vthr = –45 mV, ref = 2 ms, gext = 6 nS, rext = 1.2 kHz, gE = 36 nS, gI = 12 nS. For inhibitory cells, they are as follows: Cm = 0.2 nF, gL = 17.6 nS, VL = –70 mV, Vreset = –60 mV, Vthr = –50 mV, ref = 1 ms, gext = 1.6 nS, rext = 1.8 kHz,gE = gcue = 36 nS, gI = 12 nS. The discrete integrator network has a range of leak conductances for the excitatory cells, equally spaced from gL = 30.4 nS for cells in population-1 to gL = 40 nS for cells in population-12. The inhibitory neurons have gL = 20 nS, gext = 3 nS, and rext = 1.0 kHz. Otherwise single-cell properties are the same as in the continuous attractor network.
Short-term plasticity of excitatory synapses
All excitatory synapses exhibit short-term presynaptic facilitation and depression (Hempel et al. 2000; Varela et al. 1997). We implement the scheme described by Matveev and Wang (2000), which assumes a docked pool of vesicles containing neurotransmitter, where each released vesicle is replaced with a time constant, d. The finite pool of vesicles leads to synaptic saturation, because when the presynaptic neuron fires more rapidly than vesicles are replaced, no extra excitatory transmission is possible. Such synaptic depression contributes to stabilizing persistent activity at relatively low rates.
We assume that there is at most one vesicle release per spike; hence the release probability at any individual synapse, PR(t), is
 | (7) |
PR(t), simply scales the amplitude of synaptic transmission, as shown in Eq. 3. Similarly, we do not keep track of a discrete N(t) for every individual synapse, but assume that the several synapses between two neurons have a binomial distribution with an average number of docked vesicles, N(t) (where the brackets represent the average over this binomial distribution, with mean N(t)) and maximum N0. Hence N(t) is a continuous variable obeying
 | (8) |
PR(t) after a spike at time tspike. By averaging over the binomial distribution, we have
 | (9) |
 | (10) |
PR(t) is used in Eq. 3.
The vesicular release probability is given by the product of three gating variables, pv(t) = O1(t)O2(t)O3(t). A gating variable Oi(t) (i = 1,2,3) increases because of calcium influx triggered by an action potential, followed by a decay with time constant fi between spikes. Specifically, the following simple update rule is used: a gating variable Oi(t) (i = 1,2,3) follows
 | (11) |
Our simulations use the following values for the parameters in the continuous network: N0 = 16, d = 0.5 s, Cf1 = 0.45, f1 = 50 ms,Cf2 = 0.75, f2 = 200 ms, Cf3 = 0.9, f3 = 2 s. Excitatory neurons in the discrete attractor network are identical except with d = 0.1 s.
Network connectivity
Connection strengths between neurons depend only on their group numbers and are all-to-all. All weights are normalized by (i.e., divided by) the number of neurons in the presynaptic group, so that average network properties should be independent of the system size.
The set of weights WEI, WIE, WII all follow the same form
 | (12) |
EI determines the breadth of connections to other groups. WmaxIE, IE, WmaxII, and II have similar definitions.
The recurrent excitation has a slightly different form. First, the connections within the same group are significantly stronger than those between groups, so we define a separate set of parameters for the WiiEE =Wi. Second, the connection strengths between different groups i and j are asymmetric
 | (13) |
 | (14) |
The cross-inhibition, WcrossEI, is the strength of connection from each inhibitory population, labeled by i to the excitatory population of opposite tuning labeled by 14 – i for 2 
i 12.
The full set of parameters are as follows for the continuous network:
W0EE = 0.16, W1 = 0.244, W2 = 0.239, W3 = 0.237, W4 = 0.238, W5 = 0.239, W6 = 0.24, W7 = 0.241, W8 = 0.242, W9 = 0.243,W10 = 0.244, W11 = 0.245, W12 = 0.246; 1 = 0.5, 2 = 0.4, 3 = 0.39, 4 = 0.385, 5 = 0.385, 6 = 0.388, 7 = 0.392, 8 = 0.397, 9 = 0.402, 10 = 0.408, 11 = 0.414, 12 = 0.42; AEE = 1.5;WmaxEI = 1.65, EI = 0.25, WmaxIE = 0.5, IE = 0.2, WmaxII = 2.0, II = 0.5. WcrossEI = 0.25.
For the discrete network: W0EE = 0.14, W1 = 0.35, W2 = 0.365,W3 = 0.378, W4 = 0.39, W5 = 0.401, W6 = 0.412, W7 = 0.423,W8 = 0.434, W9 = 0.445, W10 = 0.455, W11 = 0.465, W12 = 0.475; 1 = 10, 2 = 10, 3 = 10, 4 = 10, 5 = 10, 6 = 10, 7 = 10, 8 = 10, 9 = 10, 10 = 10, 11 = 10, 12 = 10; AEE = 1.0;WmaxEI = 0.3, EI = 0.4, WmaxIE = 0.3, IE = 0.4, WmaxII = 0.5, II = 0.5; WcrossEI = 0.25.
The readout cells of the discrete network are excited by cells from the populations, i, with weights WiER as follows: W1ER = 0.45, W2ER = 0.4, W3ER = 0.35, W4ER = 0.4, W5ER = 0.25, W6ER = 0.2, W7ER = 0.2, W8ER = 0.2, W9ER = 0.2, W10ER = 0.2, W11ER = 0.2, W12ER = 0.2.
Covariance functions
In the following sections, a key quantity that enters the calculations and affects all the statistics is the covariance in spike rates, covij(r) (t1, t2) between neurons i and j (which can refer to the same neuron with i = j) at times t1 and t2. The special case covii(r) (t1, t2) is the variance in spike rates across trials for a neuron i at time t1. When analyzing data, we count the number of spikes by neuron, i, in trial, n, in specific bins of width 
t (where t = 25 ms unless otherwise stated) between t = (k – 1/2)t and t = (k + 1/2)t as Nin(k). Hence the relevant covariance in spike count becomes
 | (15) |
, explicitly in the final line. We achieve the average across different trials, n, in our computer simulations by using different initializations of the random number generator. Correlation functions and correlograms
We calculate the unnormalized cross-correlogram by averaging the covariance function over the measurement interval according to the formula (Bair et al. 2001; Brody 1998, 1999)
 | (16) |
and the total temporal window of data measurement is T. Both T and t are integral multiples of . The integration window for detection of spikes offset by is limited to T – . The above definition is for positive t. For negative , the summation limits are from k = –/ to T/ and prefactor is 1/(T + ) such that Cij(–) = Cji().
We normalize the covariance functions by dividing by the geometric mean of the variance at the two time-points (Brody 1999). This leads to the correlation coefficient, rij(t,t + ), between spike times t and t + of neurons i and j
 | (17) |
 | (18) |
|
replaced by T' to remove any artificial dependence on because of its inclusion in the integration limit. Different values of T' produce the four different curves in Fig. 4B.
|
The appropriate power spectrum for a nonstationary process is the Wigner-Ville spectrum (Flandrin 1989), which is well established in signal processing and physics (Mallat 1999). It is defined as a function of frequency, 
, and time, t, as
 | (19) |
, the time-difference ( = t2 – t1) leads to a result independent of t, the mean time [t = (t1 = t2)/2] for a stationary process.
For a nonstationary process, we can calculate a power spectrum in frequency alone, P(), independent of time, t, by averaging the Wigner-Ville spectrum across the measurement period (Flandrin 1989) as
 | (20) |
Calculation of power spectrum from spike trains
In Fig. 4C we plot the temporally averaged Wigner-Ville spectrum, P(), for neurons in our simulations, using a measurement interval of T = 10 s. Specifically, we calculated the average power spectrum by combining spikes from different neurons in a population and using Eqs. 15, 19, and 20 to give
 | (21) |
n = (n)/T, and the sums are over all binned spikes [where k is the time index for the bin that contains spikes at times t such that k – (t)/2 < t k + (t)/2)] from Npop = 25 different neurons from the same population (with identical network inputs, but different background noise) in trial- out of Ntrials. In the log-log plot of Fig. 4C, we show values for low-n up until statistical noise causes the power to be negative for some data points (in which case we cannot take the logarithm) and estimate the power coefficient, , from a straight line fit of log-power versus log-frequency through the points shown. Fano factor
The Fano factor, F(T), is a measurement of trial-to-trial variability of the spike count of an individual neuron. It is obtained from the variance in spike count during a measurement window of temporal length, T, and is expressed as a function of time
 | (22) |
Noise correlation
We evaluate the noise correlation by counting the total number of spikes during a delay of 6 s after the cue offset, for a given neuron, i, in the nth trial, as Ni,n. The noise correlation, Xij, between two neurons, i and j, is defined as (Lee et al. 1998)
 | (23) |
Ni = (
nNi,n)/Ntrial is the trial-averaged number of spikes for the ith neuron. Clearly Xii = 1, the maximum possible correlation when i = j, and in general, –1 Xij 1. The noise correlation is calculated for each cue and averaged across cues in Fig. 5, A and B.
|
 |
|
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIXGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
). The model was designed to simulate prefrontal cortical neurons in monkeys during a somatosensory delayed discrimination task (Romo et al. 1999). In this task, the monkeys were trained to discriminate the frequencies of two vibrotactile stimuli, presented before and after a delay period of 3–6 s, so that the behavioral response (a lever press to indicate whether the 1st or 2nd stimulus was at a higher frequency) depended on remembering the first stimulus frequency across the delay period. It was discovered that neurons in prefrontal and premotor cortices exhibited persistent delay activity that was tuned monotonically to the initial frequency. Such neurons maintain the necessary information for the animal to perform the delayed comparison task.
The model network contains 12,000 neurons that fire spontaneously and stochastically (with a CV averaging >1 at firing rates <10 Hz, dropping to 0.5 as firing rates increase to 40 Hz) because of noisy excitatory inputs. The neurons are arranged in two sets of 12 neuronal populations. In line with the experimental data (Romo et al. 1999), the activity of one set of populations increases with the stimulus frequency, whereas the activity of the other set decreases with increasing stimulus frequency. Hence the tuning of populations is positively monotonic or negatively monotonic, respectively. Recurrent excitation dominates within each set of populations, whereas the two oppositely tuned sets are connected by cross-inhibition (see Fig. 1B). The results comparing the mnemonic activity of our model network with the delay activity of monkey neurons are published elsewhere (Miller et al. 2003).
In our model, after a transient stimulus, neurons are able to fire persistently, with slow drift, over a range of firing rates. The feedback to each population is tuned (schematically to position C in Fig. 2A) so that each population is on the cusp of bistability, enabling the system to be close to possessing a continuous attractor (Aksay et al. 2000; Durstewitz 2003; Loewenstein and Sompolinsky 2003; Miller et al. 2003; Seung 1996; Seung et al. 2000b). The delay activity after different strengths of stimulus is shown in Fig. 2, B and D. Neurons are able to fire persistently over a wide range of rates. Initial stimuli of different intensities cause a transient increase in firing rates, before the network settles into a state of approximately constant firing rate, but with fluctuations caused by noise. We label the system in this case as quasi-continuous.
Alternatively, we can adjust parameters so that each population provides strong excitatory feedback to itself and receives less excitation from other sources (position A in Fig. 2A). In such a case, the population can be strongly bistable, often with a large gap in firing rates between a down state and an up state. If groups of neurons are bistable, but have differing excitabilities or thresholds, the network can possess a number of robust, discrete stable states (Goldman et al. 2003; Koulakov et al. 2002), with the number of states approximately equal to the number of bistable groups. Figure 2, C and E, shows the persistent activity for such a multistable network. The possible firing rates of an individual neuron have a large gap, between states where the rate is almost zero, and states where the firing rate is high (Fig. 2E, open symbols). Such a large gap in firing rates results from the strong recurrent feedback within a population, so that once the neurons in a population are able to fire a little, they excite each other with strong synaptic input, causing increasing activity in an escalating manner. Synaptic depression and receptor saturation limit the rise in firing rate to a persistent level of 40–50 Hz in our model. Any transition between the multiple discrete states in the network involves a discrete jump in the rate of at least one population of neurons, from their spontaneous value to their lowest persistent rate (across the gap in Fig. 2A). Such a large jump in firing rates of neurons in a population between stable states makes the discrete system more stable to noise fluctuations (Koulakov et al. 2002).
To produce more realistic cells in the discrete attractor network, to be compared with those in the continuous attractor network, we added a readout population to the discrete network. The readout population receives excitatory input that is a weighted sum from all populations of the same sign of tuning in the network. The tuning curve of the readout population (Fig. 2E, solid symbols) more closely resembles those of neurons in the continuous network and experimental data (Romo et al. 1999). Cells from the readout population also have a similar, rate-dependent CV (although about 10% lower) to neurons in the continuous attractor network. In all statistical analyses, we use such cells from the readout population of the discrete network.
Hence using the same framework of interconnected populations of neurons, we compared the noise and correlations of these two types of memory system. Our main results are independent of details of the model. However, they rely on noise moving the system along a continuous range of states but being unable to cause large jumps in activity between discrete states. A set of discrete attractors with small gaps in rate between states, or with large noise fluctuations, so that noise could cause spontaneous transitions between the discrete states, would show the same behavior we describe here for a continuous attractor (Miller 2006).
Random walk of rate
If a continuous attractor underlies our parametric working memory model, fluctuations in the state of the network are described by a random walk. To show this conceptual point, consider a simplified, linear model of a typical cell in a neural group whose firing rate r(t) obeys the following dynamic equation
 | (24) |
(t) is an uncorrelated noise, 
is the noise amplitude, and s is a synaptic or membrane time constant (a few to tens of milliseconds). Noise leads to an autocorrelation which decays exponentially
 |
; Egorov et al. 2002; Goldman et al. 2003; Loewenstein and Sompolinsky 2003) or through a reverberatory network (Koulakov et al. 2002; Miller et al. 2003; Seung 1996; Seung et al. 2000b). In either case, a positive feedback implies a component of the input current that increases with its own firing rate, so making the simplification of linear feedback, I = I0 + Wr, the above equation becomes
 | (25) |
eff(W) = s/(1 – GW) > s. If there is no feedback, W = 0, and in time s, the firing rate reaches its steady-state value, GI0, as in Eq. 24. With strong feedback (GW 1), one has eff >> s. In the limit eff 
, the first term in the right hand side of the equation vanishes. The firing rate integrates the input I as well as noise over time, in the sense of calculus: 
. Moreover, after the input is withdrawn, the firing rate is maintained stationary at any level (within an operating range) except for noise-induced drifts; hence the system behaves like a line attractor. The average rate is given by < r > = GI0 tI/s, where tI is the time when the input is withdrawn, whereas the trial-to-trial variance in rate increases linearly with time: (r – r)2 = At. The analogy between the random walk model and our parametric working memory model lies in the fact that the continuum of persistent activity states form a one-dimensional curved line in the space of population firing rates, as schematically shown in Fig. 1A. If the firing rates are temporarily perturbed by noise to move off this line, the network's dynamics cause the firing rates to return back to the line, but do not prevent any drifts along the continuous attractor. Such noise-induced drifts along the attractor are similar to a random walk. To assess whether the simple random walk model can quantitatively capture the behavior of our large-scale and highly nonlinear network model, we analyze the fluctuation properties of population firing rates and spike counts.
Random walk analysis of a line attractor
A closer analysis of the effects of noise on a system with a continuous attractor reveals that the correlation functions are determined by the products of gradients of the tuning curves of two neurons (Ben-Yishai et al. 1995; Pouget et al. 1998). The argument is shown in Fig. 1A, where the top three figures depict the tuning curves (average firing rate as a function of stimulus) for three neurons. The first two are positively monotonic; the third is negatively monotonic. When the system has a continuous attractor, the only fluctuations that are not rapidly damped out are those along the attractor. The attractor can be seen as a curve in a plot of the firing rate of one neuron versus another, with each point in the curve corresponding to the two firing rates of persistent activity after a particular stimulus. Noise has the same effect as the stimulus in shifting the firing rates along the attractor. Two positively monotonic neurons [with firing rates r1,r2 and tuning curves f1(s),f2(s)] either both increase or decrease their firing rates together, after a larger or smaller stimulus, and after noise fluctuations. Hence their noise correlations are expected to be positive.
If we consider the firing rate of mnemonic persistent activity as a function of stimulus, s, and noise, (t), expanded to first order in the noise, we have
 | (26) |
 | (27) |
, which is given by the product of the noise terms in Eqs. 26 and 27. Hence the cross-correlation is proportional to the product of gradients of the two tuning curves [(df1/ds). (df2/ds)] (Ben-Yishai et al. 1995; Pouget et al. 1998). If the gradients of tuning curves of two neurons have opposite sign [as in Fig. 1A; (df1/ds)(df3/ds) < 0], the noise correlation is negative between the corresponding neurons. That is, a fluctuation along the attractor simultaneously causes the positively monotonic neurons to increase their firing rates and negatively monotonic neurons to decrease their rates or vice versa. Similarly, for an individual neuron, the covariance function is proportional to (dfi/ds)2, so all effects of noise are greater after cues that fall on a steep part of a neuron's tuning curve.
Noise during the stimulus presentation leads to variation in the initial encoded values of s and has a correlated effect on the firing rate of neurons in the same way as noise during the delay. The main effect of both stimulus noise and drift noise in a line attractor is encoded in the one-dimensional variable, s. The firing rate of each neuron responds to the network noise by an amount proportional to the square of the gradient of its tuning curve. In the preceding analysis, we assumed that fluctuations are small compared with changes in the gradient of the tuning curve, so dfi/ds remains constant for each neuron for a complete trial, allowing Eqs. 26 and 27 to be expanded only to first order. If the first-order approximation is valid, the proportionality constant, A, is the same for both stimulus noise and drift noise. That allows us to write the total noise at time t as
 | (28) |
(dfi/ds)2.
Other results for a random walk are described in detail in the APPENDIX. We summarize the key features that can be compared with our network results here. The covariance function, cov(t1,t2), for a random walk starting at t = 0 is independent of the time difference, |t2 – t1|, and instead equals the variance of the process given above in Eq. 28 where t is the earlier of t1 or t2 (Gillespie 1992) (cf. Eqs. A4 and A5).
The other statistical measures presented in this paper are derived from the covariance or variance functions. The time-average of the covariance function is proportional to the integration limits used in the time-averaging (see Eqs. A6–A9). The average (Eq. 16), leads to a linear decay in ||
 | (29) |
, but linearly increasing with T'
 | (30) |
A hallmark of random walk behavior is an inverse-square power law for the Wigner-Ville power spectrum (Eq. 19) (Flandrin 1989). When averaged over a finite time interval, T, specific frequencies, n = n
/T should be used, so that the inverse-square power law is apparent in the averaged spectrum. We find that initial variance in the random walk leads to oscillations in the power spectrum between even and odd frequencies (Eqs. A11 and A12), such that
 | (31) |
 | (32) |
The Fano factor, F(T), for a random walk, also shows power-law behavior, increasing quadratically with time. Initial variance in the starting points contributes a linear term, such that for a Poisson-like renewal process
 | (33) |
Correlation function
Fluctuations are commonly quantified by autocorrelation functions of single neurons and cross-correlation functions between them (Bair et al. 2001; Brody 1999; Perkel et al. 1967a,b; Singer and Gray 1995). The autocorrelation function is a function of two times (1 for each spike) and should only be written as a function of the time-difference, or time-lag, , if the underlying process is stationary (Gardiner 1985). Many processes measured in neuroscience are not stationary, so the shuffle-correction is used to remove the principle effects of systematic variations in the average rate. We use the term correlogram to denote the shuffle-corrected covariance function, normalized by the product of SD in spike count (producing a correlation coefficient), averaged across the measurement interval (Brody 1999). We use the term unnormalized correlogram when referring to a temporal average of the covariance function.
Trial-to-trial differences in the underlying rate have been shown to give rise to artifacts in cross-correlograms at short time scales (Brody 1998, 1999). In contrast, we report that the correlation functions are affected on long time scales by the linear increase in variance of firing rate as a function of time that arises from the temporal integration of noise in a quasi-continuous attractor.
The autocorrelograms for three cells in our graded memory network are shown in black, green, and blue in Fig. 3A. Strikingly, temporal correlations between spikes persist for many seconds and seem to show a gradual, linear decrease with the time lag. In contrast, a network with robust, discrete states does not show any long-term correlations (Fig. 3A, red curve). In fact, the autocorrelograms for all neurons fall to zero during a time scale on the order of the synaptic time constant (100 ms for NMDA receptors in this work; Fig. 3A, inset).
The autocorrelogram in Fig. 3A was calculated by averaging over memory activity after different stimuli. If we compute the correlations separately after each stimulus, we find that only a small subset of stimuli result in such large, long-term persistence in the correlations. This is because, if the stimulus lies on a relatively flat part of a neuron's tuning curve, small network fluctuations do not cause significant changes in the firing rate of that neuron. Tuning curves of neurons in our network are not linear in general, but resemble sigmoidal functions, as seen for the majority of neurons in the experimental data (Miller et al. 2003; Romo et al. 1999). It is the steepest part of the tuning curve that is most sensitive to fluctuations. If the tuning curve of a neuron is linear (unlike those seen in our model), the effect of fluctuations during the delay should be independent of stimulus.
In Fig. 3B, we calculated the unnormalized autocorrelogram separately for delay activity after three stimuli, corresponding to the positions X, Y, and Z in Fig. 2D. To obtain better statistics, we averaged cross-correlations between 25 neurons in a single population. Because the population fires asynchronously and all neurons within one population receive the same network input, their autocorrelation functions are essentially the same. The average cross-correlation between all pairs of neurons in the same population is identical to the average autocorrelation of those neurons, except at very short time scales (less than 
25 ms). It is only the -function peak at = 0 and the dip in the autocorrelation at short (because of the refractory period and recovery from reset; Fig. 3A) that distinguishes autocorrelations from cross-correlations between cells within a population of our network.
Figure 3C shows that the cross-correlogram is greater between neurons with similar tuning curves (such as 2 neurons in the same population, 2–2, or neighboring population, 2–3) than between neurons with very different thresholds (2–12). The cross-correlogram is negative between neurons of opposite tuning (2–12*) and is temporally asymmetric. The asymmetry is a sign that firing rates of one population can drift for several seconds before the change in rate is strong enough to affect the populations with opposite sign of tuning. This is because cross-inhibition is not the dominant feature in our network and suggests that, on a time scale of a few seconds, our system resembles more closely two weakly coupled continuous attractors (where drift in 1 attractor provides input to the other) than a single set of stable points.
In all cases, autocorrelograms and cross-correlograms do not decay exponentially as exp(–t/) with a characteristic time constant, (see the text after Eq. 24), but decay more linearly, as is typical for a random walk.
Increasing variance of rate
We found that the variance of single-trial population firing rates (calculated with 100 trials) increases approximately linearly from the stimulus offset for up to 10 s of the mnemonic period (Fig. 4A) . The nonzero intercept at t = 0 reflects trial-to-trial noise in the system's response to the stimulus.
The linear time course with nonzero intercept for the firing rate variance is consistent with a random walk (see Eqs. A4 and A5 and Noise during the stimulus), with an initial distribution of starting points caused by variability of neuronal response during the stimulus presentation (Eq. A2). Note that the variance of neurons in the discrete system remains constant and low (Fig. 4A, red curve).
Covariance independent of time lag
Such a random walk also provides an explanation for the observed near-linear behavior of the correlograms reported in Fig. 3. Because random walks are nonstationary processes (the variance increases with time), correlation functions depend separately on the two times of measurement (t1,t1 + ) and not just on their difference, (see METHODS). Hence when evaluating correlations as a function of the time lag, , by integrating over t1, the precise measurement interval affects the result. It can be shown mathematically that a linear increase in the variance of rates as a function of time leads to a linear increase in the unnormalized correlogram as a function of the integration interval (Miller 2006; Saleh 1978). In the standard formulation, the integration interval is (T – ||), where is the time-lag and T is the total time of measurement. With fixed T, this leads to a linear decrease with , as seen in Fig. 3B.
In Fig. 4B, we tested whether the unnormalized autocorrelogram is truly a function of the measurement interval using an alternative integration window to calculate the temporal average (the sum over k in Eq. 16). We calculate the unnormalized autocorrelogram by including spikes over a range of T' + for each value of so that the integral over t1 is over a fixed value, T' = T – max, rather than T – . Hence the measured correlation interval is independent of and proportional to T'. We plotted curves for the average unnormalized autocorrelogram of 100 neurons (25 neurons in each of 4 populations) using increasing values of T' (4, 8, 12, and 16 s for the ascending curves). For a perfect random walk, the curves in Fig. 4B would be constant with and equally spaced along the y-axis (see METHODS). The curves do have a y-offset that increases monotonically with T', whereas they vary little with , so the network's activity is behaving qualitatively like a random walk during the delay.
Power law of power spectra
Furthermore, for a random walk, the power spectrum, P(n), of the spike train should include a delta-function at the origin and at low frequencies an inverse-square law decay, such that P(n) n– where a equals 2 (Flandrin 1989) and n = n/T. A pure power law decay should appear as a straight line on a log-log plot, with the slope of the line giving the exponent (Fig. 4C). We estimated the exponent, , for several cells by fitting the linear region of log-log power spectra. We fitted curves separately through odd and even data points, because these were offset due to noise during the stimulus (see Eq. A12 and Noise during the stimulus). The fits yielded values ranging from = 1.6 to = 2.1 ( = 1.77 and 1.66 for the odd and even frequencies, respectively, of population-2, shown in Fig. 4C). The fact
, by stimulus is the lowest for neural population 1 and highest for neural population 12. The 2 networks interact through pyramid-to-interneuron connections, resulting in cross-inhibition. Stimulus inputs are given equally to all excitatory cells within a network, but they differ across the 2 oppositely tuned networks.
