## Abstract

Delay-period activity of prefrontal cortical cells, the neural hallmark of working memory, is generally assumed to be sustained by reverberating synaptic excitation in the prefrontal cortical circuit. Previous model studies of working memory emphasized the high efficacy of recurrent synapses, but did not investigate the role of temporal synaptic dynamics. In this theoretical work, I show that biophysical properties of cortical synaptic transmission are important to the generation and stabilization of a network persistent state. This is especially the case when negative feedback mechanisms (such as spike-frequency adaptation, feedback shunting inhibition, and short-term depression of recurrent excitatory synapses) are included so that the neural firing rates are controlled within a physiological range (10–50 Hz), in spite of the exuberant recurrent excitation. Moreover, it is found that, to achieve a stable persistent state, recurrent excitatory synapses must be dominated by a slow component. If neuronal firings are asynchronous, the synaptic decay time constant needs to be comparable to that of the negative feedback; whereas in the case of partially synchronous dynamics, it needs to be comparable to a typical interspike interval (or oscillation period). Slow synaptic current kinetics also leads to the saturation of synaptic drive at high firing frequencies that contributes to rate control in a persistent state. For these reasons the slow NMDA receptor-mediated synaptic transmission is likely required for sustaining persistent network activity at low firing rates. This result suggests a critical role of the NMDA receptor channels in normal working memory function of the prefrontal cortex.

- working memory
- prefrontal cortex
- persistent activity
- NMDA receptor
- synaptic dynamics
- short-term plasticity
- rate control
- synchronization
- spiking neuron model

Working memory is a fundamental cognitive function, by virtue of which information can be actively retained for seconds and used in the brain (Baddeley, 1986; Fuster, 1988; Goldman-Rakic, 1995). Its neuronal correlate, delay-period activity, has been widely documented by unit recording studies of behaving monkeys (Fuster and Alexander, 1971; Kubota and Niki, 1971;Miyashita and Chang, 1988; Gnadt and Andersen, 1988; Funahashi et al., 1989; Miller et al., 1996; Chafee and Goldman-Rakic, 1998; Rainer et al., 1998; Romo et al., 1999). For example, in a visuospatial delayed-response experiment (Funahashi et al., 1989), the animal's delayed saccadic eye movement is guided by the short-term memory of a visual cue. Neurons in the dorsolateral prefrontal (PFC) cortex were found to display elevated firing activity during the entire delay period. This persistent activity is tuned to the spatial location of the cue in some cells, but not in other cells. Therefore, there are two distinct aspects of the mnemonic coding by the PFC cells: the persistent nature of the delay-period activity and the formation of the tuned “memory field”.

It is generally assumed that persistent activity is sustained by some kind of reverberating discharge within a recurrent neural network (Hebb, 1949; Amit, 1995). The characteristic horizontal connections found in the superficial layers II–III of the dorsolateral PFC may provide the anatomical substrate for such a recurrent circuit (Levitt et al., 1993; Kritzer and Goldman-Rakic, 1995). However, it remains unknown what are the realistic synaptic properties and circuit dynamics that are required for a robust network-induced persistent activity. Indeed, most previous model studies used simple firing-rate models (Wilson and Cowan, 1973; Amari, 1977; Zipser et al., 1993; Amit et al., 1994; Camperi and Wang, 1998; Moody et al., 1998). Amit and collaborators (Amit et al., 1990; Amit and Tsodyks 1991; Amit and Brunel, 1997) used leaky integrate-and-fire (LIF) spiking neurons but did not take into account realistic postsynaptic current time courses.

In this paper I present a network model of spiking neurons in which synapses are endowed with realistic gating kinetics, based on experimentally measured dynamical properties of cortical synapses. I will focus on how delay-period activity could be generated by neuronally plausible mechanisms; the issue of memory field formation will be addressed in a separate study. A main problem to be investigated is that of “rate control” for a persistent state: if a robust persistent activity necessitates strong recurrent excitatory connections, how can the network be prevented from runaway excitation in spite of the powerful positive feedback, so that neuronal firing rates are low and comparable to those of PFC cells (10–50 Hz)? Moreover, a persistent state may be destabilized because of network dynamics. For example, fast recurrent excitation followed by a slower negative feedback may lead to network instability and a collapse of the persistent state. It is shown that persistent states at low firing rates are usually stable only in the presence of sufficiently slow excitatory synapses of the NMDA type. Functional implications of these results for the role of NMDA receptors in the PFC working memory function are discussed.

## MATERIALS AND METHODS

*The leaky integrate-and-fire model*. To simulate a local recurrent cortical network, I used a network model of leaky integrate-and-fire neurons (Tuckwell, 1988), with either all-to-all or sparse connectivity. Such a network can be viewed as a cortical cell assembly that stores a particular memory item. As a result of Hebbian learning, the internal excitatory recurrent connections are strong and homogeneous, whereas the interactions between this cell assembly and the rest of the circuit are relatively weak and are neglected.

The network model consists of two populations of neurons (N_{e} pyramidal cells and N_{i} inhibitory interneurons). Each pyramidal cell obeys the following equation:
Equation 1
Equation 2where C_{m} is the capacitance, I_{app} represents the afferent input, and the leak current I_{L} = g_{L}(V_{m}− V_{L}). I_{AHP} = g_{AHP}[Ca^{2+}](V_{m} − V_{K}) describes a calcium-activated potassium current for spike-frequency adaptation. [Ca^{2+}] is incremented by an amount α_{Ca} with each spike discharge, and decays with a time constant τ_{Ca} afterwards (cf. Treves, 1993; Y. H. Liu and X.-J. Wang, unpublished observations). I_{syn,ee} and I_{syn,ie} are the recurrent synaptic inputs from pyramidal cells and interneurons, respectively.

A spike is discharged each time V_{m} is driven to reach a firing voltage threshold V_{th}. Then V_{m} is reset to V_{reset} and stays there for an absolute refractory period τ_{ref}. The intrinsic parameters were calibrated based on the intracellular data of cortical pyramidal neurons (McCormick et al., 1985; Mason and Larkman, 1990; Troyer and Miller, 1997): C_{m} = 0.5 nF, g_{L} = 0.025 μS (so that the time constant τ_{m} = C_{m}/g_{L} = 20 msec); V_{L} = −70, V_{th} = −52, V_{reset} = −59 (in mV); τ_{ref} = 2 msec. The frequency–current curve of an isolated cell has a current threshold I_{c} = g_{L}(V_{th} − V_{L}) = 0.45 nA. For the adaptation current V_{K} = −85 mV, α_{Ca}= 0.2 μM, τ_{Ca} = 80 msec (Helmchen et al., 1996), and g_{AHP} will be specified in the text whenever it is not zero.

The interneuron model represents fast-spiking GABAergic cells that do not display spike-frequency adaptation (McCormick et al. 1985). Each interneuron obeys the equation:
Equation 3which is similar to Equation 1, except that I_{AHP} is absent, and I_{syn,ei} is the recurrent synaptic input from pyramidal cells. Mutually inhibitory interactions among interneurons were not included. The parameter values for the interneurons are (cf.McCormick et al., 1985) C_{m} = 0.2 nF, g_{L} = 0.02 μS (τ_{m}= C_{m}/g_{L} = 10 msec); V_{L} = −65, V_{th} = −52, V_{reset} = −60 (in mV); τ_{ref} = 1 msec. The frequency–current curve of an isolated interneuron has a current threshold I_{c} = 0.26 nA.

*Synaptic kinetics and short-term depression.* The EPSC originating from a presynaptic pyramidal cell consists of two components, I_{AMPA} and I_{NMDA}. The AMPA receptor (AMPAR)-mediated current I_{AMPA} = g_{AMPA}s(V_{m} − V_{E}), with V_{E} = 0 mV. The gating variable s (the fraction of open channels) is described by two first-order kinetics:
Equation 4
Equation 5where the sum is over presynaptic spike times. The scaling factor φ controls the speed of synaptic kinetics without affecting the steady state, φ = 1 unless specified otherwise. For the AMPAR channels, I used τ_{x} = 0.05 msec and τ_{s} = 2 msec (the time-to-peak is ∼0.2 msec); α_{x} = 1 (dimensionless), and α_{s} = 1 (in msec^{−1}). The NMDA receptor (NMDAR)-mediated current I_{NMDA}= g_{NMDA}s(V_{m} − V_{E})/(1 + [Mg^{2+}] exp(−0.062 V_{m})/3.57) (Jahr and Stevens, 1990), with a voltage dependence controlled by the extracellular magnesium concentration [Mg^{2+}] = 1.0 mm. The gating variable s obeys the same types of equations (Eqs. 4, 5), but with τ_{x} = 2 msec and τ_{s} = 80 msec (the time-to-peak is ≃8 msec).

This model of excitatory synapses was chosen for three reasons. First, it is based on a plausible kinetic scheme (Wang and Rinzel, 1992;Destexhe et al., 1994). In response to a presynaptic spike, the time course of s has a smooth rising phase and an exponential decay with time constant τ_{E} = τ_{s}/φ, that can be matched to the experimental data (Hestrin et al., 1990a; Lester et al., 1990). Second, there is temporal summation and, if the presynaptic firing frequency is high compared to 1/τ_{E}, s will saturate in the steady state (s ≤ 1) (Fig. 1). The saturation effect is much more significant for the slow NMDAR-mediated EPSC than for the fast AMPAR-mediated EPSC and has important implications for the network dynamical behavior. Finally, the model is sufficiently simple to allow detailed analysis of the network activity.

The IPSC originating from an interneuron is assumed to be mediated by GABA_{A} receptors (GABA_{A}Rs), I_{GABA} = g_{GABA}s(V_{m}− V_{I}), with V_{I} = V_{L} = −70 mV (“shunting inhibition”). The gating variable s obeys a simple first-order kinetics with saturation (Wang and Rinzel, 1992):
Equation 6with α_{I} = 0.9 and τ_{I} = 10 msec. The superscript in t_{j}^{−} indicates that the increment of s by a spike should be calculated using the value of s immediately before the spike on the right hand side of the equation, Δs = s(t_{j}^{+}) − s(t_{j}^{−}) = α_{I}(1 − s(t_{j}^{−})).

Most simulations were done with all-to-all connectivity. In that case a neuron receives synaptic inputs from all neurons in the network, and the summation of synaptic currents is normalized by the number of neurons N. Sparse connectivity was also considered (see Fig.11). There, the coupling between neurons is randomly assigned, with an average number of synapses per neuron M_{syn}(which is much smaller than N ), and the summation of synaptic currents is normalized by M_{syn}. The probability that a pair of neurons are connected in either direction is p = M_{syn}/N.

In some of the model simulations, short-term depression was incorporated for the pyramid-to-pyramid recurrent excitatory synapses (Markram and Tsodyks, 1996; Abbott et al., 1997; C. M. Hempel, K. H. Hartman, X.-J. Wang, G. G. Turrigiano, and S. B. Nelson, unpublished observations). Short-term depression is assumed to be caused by transmitter vesicle depletion at the presynaptic terminals (Stevens and Wang, 1995). It is introduced into the synapse model as follows. The parameter α_{x}, which mimicks the amount of transmitter release per spike, is multiplied by a quantity D (the fraction of available vesicles). D obeys the dynamical equation (Abbott et al., 1997):
Equation 7That is, D is reduced by a factor (1 − p_{υ}) for each spike, ΔD = D(t_{j}^{+}) − D(t_{j}^{−}) = −p_{υ}D(t_{j}^{−}), or D(t_{j}^{+}) = (1 − p_{υ})D(t_{j}^{−}). It recovers toward 1 with time constant τ_{D} in the absence of stimuli. In a simple biophysical model of vesicle depletion in which the release probability is proportional to the number of available vesicles, p_{υ} is identified with the release probability per vesicle (Wang, 1999). I used τ_{D} = 500 msec and p_{υ} = 0–0.35.

*Asynchronous States.* In this work, persistent activity is assumed to be achieved by a bistability between a rest state and an active state of the network. We shall see that the persistent activity often occurs as an asynchronous network state, in which the discharges of neurons are spread out in time uniformly so that at any time there is a same fraction of neurons firing (Amit and Tsodyks, 1991; Abbott and van Vreeswijk, 1993; Gerstner, 1999). In the presence of the voltage dependence of the NMDAR channels, the nonlinear LIF model cannot be solved explicitly, and the analysis of the asynchronous states is intractable. However, as we shall see, none of our conclusions in this work depends on the voltage sensitivity of the NMDAR-activated conductance. Therefore, the calculations of the asynchronous state were done with [Mg^{2+}] = 0.

The firing rates R_{E} and R_{I} of pyramidal cells and interneurons in an asynchronous state were calculated as follows. Let us denote the average synaptic drives by s_{E} and s_{I}. Each of the two is an average over neural population, and is constant in time for an asynchronous state. It is the same as the time average of each individual s(t) over a period 1/R. For s_{E} (Eqs. 4,5), an approximation is obtained by substituting ∑_{j}δ(t − t_{j}) with R. The steady state is:
Equation 8where ν = α_{x}α_{s}τ_{x}τ_{s}msec^{−1}. This approximation is accurate when the synaptic current kinetics are sufficiently slow (Ermentrout, 1994), hence reasonable for the NMDAR channels (Fig.1*D*). On the other hand, it is also correct as long as saturation is negligible, which is the case for the fast AMPAR channels. The average s_{E}only depends on the product ν and is independent of the scaling factor φ. It becomes nonlinear in R_{E} at R_{E} ≥ 1/ν and saturates at R_{E} ≫ 1/ν. For AMPAR channels ν = 0.1 msec^{−1}, 1/ν = 10 kHz; so s_{AMPA} does not saturate at realistic firing rates, s_{AMPA} ≃ νR_{E}. For NMDAR channels ν = 160 msec^{−1}, 1/ν = 6.25 Hz, and s_{NMDA} is a highly nonlinear function of R_{E}.

With short-term depression (p_{ν} ≠ 0), the parameter α_{x} is multiplied by the steady-state value of D. It is worth noting that the amount of synaptic transmission is given by the value of D immediately preceding a spike (denoted by D_), and not the time average over a period. The steady state value of D_ is (Abbott et al., 1997; Wang, 1999):
Equation 9where the approximation is obtained by ∑_{j}δ(t − t_{j}^{−}) = R_{E} in Equation 7.

For the GABA_{A}R-activated synaptic drive, the average was calculated over a periodic firing pattern of rate R_{I} (compare Eq. 6):
Equation 10At realistic firing rates, the steady-state approximation obtained by ∑_{j}(t − t_{j}^{−}) = R_{I} in Equation 6, s_{I} = α_{I}τ_{I}R_{I}/(α_{I}τ_{I}R_{I}+ 1), is not accurate for the moderately slow IPSCs.

Given s_{E}(R_{E}), the equation for an interneuron is the same as that of a single LIF neuron,
Equation 11with g̃_{L} = g_{L} + g_{syn,ei}s_{E}, andĨ = I_{app} − g_{L}(V_{reset} − V_{L}) − g_{syn,ei}s_{E}V_{reset}.

For a constant input current I_{app}, the firing rate is given by:
Equation 12which is itself a function of s_{E}(R_{E}) (i.e. the interneurons are driven by recurrent excitation). Similarly, the voltage equation for a pyramidal cell can be solved for a constant input current. The adaptation current has a steady-state average I_{AHP} = g_{AHP}[Ca^{+2}]_{aυ}(V_{m}− V_{K}), where [Ca^{+2}]_{aυ} = α_{Ca}τ_{Ca}R_{E} according to Equation 2 with ∑_{j} δ(t − t_{j}) = R_{E}. The same formula in Equation 12 applies to R_{E}, except thatg̃_{L} = g_{L} + g_{AHP}[Ca^{+2}]_{aυ} + g_{syn,ee}s_{E} + g_{syn,ie}s_{I}, andĨ = I − g_{L}(V_{reset} − V_{L}) − g_{AHP}[Ca^{+2}]_{aυ}(V_{reset}− V_{K}) − g_{syn,ee}s_{E}V_{reset} − g_{syn,ie}s_{I}(V_{reset} − V_{I}).

In simulations, noise was added by including a random component I_{λ} = i_{λ}s_{λ} in the external current, I_{app} = I_{0}+ I_{λ}; I_{0} is a constant current, and I_{λ} is a stochastic synaptic current of the AMPA type. With a Poisson input train of rate λ, s_{λ} is incremented by 1 with each input and decays with a time constant τ_{λ} = 2 msec. At a high rate λ, this Poisson current is approximated by a Gaussian white noise with a mean μ = i_{λ}λτ_{λ} and a variance ς^{2} = i_{λ}^{2}τ_{λ}λ. Unless noted otherwise, i_{λ} = 0.06 nA and λ = 2500 Hz for pyramidal cells; and i_{λ} = 0.04 nA and λ = 2000 Hz for interneurons. Given a fixed μ = 0.06 × 2.5 × 2 = 0.3 nA, the mean input current to pyramidal cells I = I_{0} + μ can be varied by changing I_{0}, whereas the noise amplitude remains the same.

In the presence of noise, the neural discharges are described by the first-passage times across the firing threshold (Ricciardi, 1977), instead of Equation 12. The expression for the firing rate is:
Equation 13where er f(x) = (2/
) ∫_{0}^{x} exp (−x′^{2})dx′ is the error function, V_{ss} = I_{eff}/g̃_{L} andτ̃_{m} = C_{m}/g̃_{L} (g̃_{L}as given above). The effective current I_{eff} = I_{0} + μ + g_{L}V_{L} for interneurons, and I_{eff} = I_{0} + μ + g_{L}V_{L} + g_{AHP}[Ca^{+2}]_{aυ}V_{K}+ g_{syn,ie}s_{I}V_{I} for pyramidal cells.

The neural firing rates of the asynchronous state were approximately computed in two steps. First, Equation 13 applied to R_{I} is a function of R_{E}, R_{I} = g(R_{E}). Then, Equation 13for R_{E} becomes self-consistent,
Equation 14which is solved to yield R_{E}. Note that f is the input–output function of a LIF neuron, another way of writing Equation 14 is:
Equation 15where I_{tot} depends on I_{app}, R_{E} and R_{I}.

In numerical simulations, the initial condition can be prescribed to be near the asynchronous state, by assuming that the neural output patterns are periodic with the phases uniformly distributed in a time period [0, T = 1/R] (Abbott and van Vreeswijk, 1993). Note that in the presence of noise, the time course of the neural membrane potential is not exactly periodic. However, this initial condition should be close to the actual asynchronous state. If the latter is a stable attractor, with this initial condition the network should quickly converge to it.

*Numerical integrations.* The model was numerically integrated using a second order Runge–Kutta method, with an interpolation procedure to determine the spike times (Hansel et al., 1998). The time step dt = 0.02 − 0.05 msec. Typically I used N_{e} = 1000 and N_{i} = 200, some conclusions were checked with N_{e}= 5000 and N_{i} = 1000 (N_{i}/N_{e} = 20%).

In simulations, the network activity was measured by the instantaneous firing rate R_{E}(t) of the pyramidal cell population as follows. The time was divided into small bins (Δt = 1–10 msec were used). Then,

For example, in an asynchronous state R_{E}would be constant in time and equals the firing rate of each individual cell. A coherent network oscillation would be reflected by a rhythmic time course of R_{E}(t).

## RESULTS

### NMDA receptor channels and persistent activity at low rates

Persistent activity is produced by an excitatory neural network, when the recurrent synapses are sufficiently strong. In Figure2*A*, the network is initially in a rest state. In response to a transient input pulse, neurons start to discharge spikes that activate recurrent synapses, which in turn elicit more spikes. This positive feedback loop between the spike firing and the recurrent synaptic drive leads to a self-sustained network activity, outlasting the input. In the persistent state, neurons fire spikes asynchronously in time: at any given moment there is always a fixed fraction of cells firing. Therefore, the synaptic drive to each cell is tonic (constant in time). Moreover, the average firing rate of neurons is ∼40 Hz, within the physiological range of the persistent activity of PFC cells during the delay period (Funahashi et al., 1989; Rainer et al., 1998). The network is turned off by a brief hyperpolarizing input, from the persistent state back to the rest state. In this simulation, the leak conductance g_{L} differs from cell to cell according to a Gaussian distribution. Cells with the smallest g_{L} values are the most excitable and display spontaneous firing in the rest state; whereas cells with the largest g_{L} values are the least excitable and only show transient responses to the input pulse but no persistent activity (Fig.2*A*).

The bistability between the rest state and active state is a network phenomenon. As illustrated in Figure 2*B*, during persistent activity, a neuron can be temporally hyperpolarized by an applied current pulse, but its activity resumes itself immediately after the perturbation, because the firing of any single neuron is sustained by synaptic inputs from the circuit. Such a manipulation would be feasible experimentally only with intracellular recording from a behaving animal during a working memory task. The prediction is that if bistability is not a single cell property but is instead induced by the network circuit, a hyperpolarizing current pulse should be incapable of switching a neuron off from its persistent activity.

In model simulations, the NMDAR-mediated synaptic transmission was necessary to generate network persistent activity, at low firing rates such as in Figure 2. For the purpose of illustration, consider first the simplified situation of a perfectly synchronous network state in which all neurons behave exactly the same in time. Therefore, the population of identical excitatory neurons can be reduced to a single neuron endowed with an autapse (Fig.3*A*). Suppose that the synaptic transmission is of the NMDA type (decay time τ_{E} = 80 msec). The cell is switched onto a firing state by a transient input. At the end of input pulse, the NMDAR-mediated current decays slowly, and after the time span of an interspike interval (ISI), it remains large enough to trigger another spike, which in turn generates more EPSC. This process between the spiking and synaptic activation can continue indefinitely, provided that the decay of the NMDAR-mediated current is not too fast compared to a typical ISI, i.e. the τ_{E}/ISI ratio is sufficiently large. Otherwise, if the synaptic current generated by one spike decays back to zero before the next spike is triggered, the cell will return to the rest state instead. This is shown in Figure3*B*, where the synapse is now assumed to be of the AMPA type (τ_{E} = 2 msec). The peak AMPAR-mediated EPSC here is ∼10× that of the NMDAR-mediated EPSC in Figure3*A*, but it decays rapidly between spikes during the input pulse and does not give rise to persistent activity. Using a considerably stronger synaptic conductance, the AMPAR-mediated current can be large enough to generate a persistent activity, but at a very high firing rate, so that the τ_{E}/ISI ratio is again large (see below).

The above argument applies to the network, if the neural firing patterns are partially synchronous. For example, this can happen because of the interplay between rapid synaptic excitation and slower inhibition in a two-population network of pyramidal cells and interneurons (Fig. 4). Powerful AMPAR-activated synapses between pyramidal cells amplify the network activity, which is damped afterwards by recurrent inhibition, leading to synchronous network oscillations at ∼8 Hz (the oscillation frequency ranges from 8 to 65 Hz, when the pyramid-to-interneuron coupling strength is varied gradually). Note that the AMPAR-activated synaptic drive s_{AMPA} fluctuates between zero and a peak level during the oscillation. Without NMDAR channels, clearly this synchronous network state would not be self-sustained, because when s_{AMPA} is almost zero the network would have to collapse back to the rest state. The slow NMDAR-mediated current does not decay back to zero during the waning phases of the network oscillation. As a result, the tonic component of the NMDAR-activated synaptic drive s_{NMDA} can sustain a synchronous persistent state. Here, the requirement is that the oscillation period T must not be too long compared to the NMDAR channel decay time constant (τ_{E}/T must be large).

Therefore, by virtue of its temporal summation, NMDAR channels (but not AMPAR channels) can provide sufficient tonic drive to maintain a synchronous persistent state at low rates. On the other hand, if the persistent state is asynchronous, a tonic synaptic drive can be realized by a spatial summation over neurons. In the latter case, because the synaptic drive is constant in time regardless of the τ_{E}/ISI ratio, it would seem that the fast AMPAR channels alone might be sufficient to maintain a persistent network state at any firing rate. As we will see below, this is not the case because of the problem of rate control with the AMPAR channels.

### Frequency–current relation of a bistable network

Persistent activity in our network model is realized as a bistability between a rest state and an active state, where the network can be switched on from the rest state by a transient stimulus and remains in the persistently active state afterwards. Consider for example the case where synaptic connections are mediated by the fast AMPARs. For a fixed synaptic coupling (g_{AMPA} = 1.05) and a given external drive (I = 0.3 nA), the neuronal firing rate of an asynchronous network is given by the nonlinear equation R = f(I_{tot}(R)) (Eq. 15in Materials and Methods). The function f is the neuronal input–output relation, and the total input I_{tot}is a function of R caused by the recurrent synaptic interactions. When the left and right hand sides of the equation are plotted on a same graph, the solutions for R correspond to the intersection points of the two curves. As shown in Figure5*A* (top panel, solid curve), there are three states of different firing rates: a rest state (in which synapses are not activated), an active (persistent) state, and a middle state (which is always unstable, thus not observable in network simulations). The instability of a steady state can be intuitively understood as follows. When f(I_{tot}) < R, the total current acts to decrease firing, whereas when f(I_{tot}) > R the total current acts to increase firing. Therefore, if the rate R happens to be slightly higher than the middle steady state, f(I_{tot}) > R and R will increase further; whereas if R is lower than the middle state, f(I_{tot}) < R and R will decrease further. In either case the system will drift away, and the middle steady state is not stable against small perturbations.

The bistability occurs within a certain range of the I values (Fig. 5*A, top panel*). If the external drive is too small (I = 0.1), the combined external and recurrent drive is not sufficient to maintain a persistent state. On the other hand, if it is too large (I = 0.5), the rest state no longer exists. By plotting the steady states as function of I, an S-shaped frequency–current curve is obtained for a bistable asynchronous network (Fig. 5*A, bottom panel*). Let us denote by I_{a} and I_{b} the two I values delimiting the bistable range. I_{a} is the smallest I value for an active state, and I_{b} is the largest I value for the rest state. I_{b} ≃ 0.4 nA is close to the threshold current for an isolated neuron, because recurrent synapses are not activated in the rest state. The firing rate of the active state increases with I; the lowest possible rate corresponds to I_{a}, at the left-knee of the curve. In our example neuronal firing rates of persistent activity are above 110 Hz, much higher than those observed in the PFC neurons (10–50 Hz).

Can the firing rate of persistent activity be reduced by weaker recurrent synaptic connections? In Figure 5*B* are shown the frequency–current curves of the network at various coupling strengths (g_{AMPA}). We see that bistability becomes possible only with sufficiently strong g_{AMPA}. With larger g_{AMPA}, persistent state can be realized at smaller I (I_{a} shifts to the left), so the bistability range (I_{b} − I_{a}) is wider (the persistent state is more robust). On the other hand, the lowest firing rate of a persistent state (at I_{a}) dramatically increases with g_{AMPA} (Fig. 5*B, filled square*). Therefore, there is a tradeoff between the lowest firing rate possible and the robustness of the phenomenon: if we require that the bistable range be reasonably large (at least 0.1–0.3 nA, for example), the firing rate of a persistent state is always 100–200 Hz or higher. Furthermore, the stability of the active state is not guaranteed. Indeed, the persistent state close to I_{a} is usually not observed in direct simulations of the network model, presumably because it is not stable in the presence of noise. The stability issue will be discussed in more detail below, when negative feedback processes are included.

In contrast to the case with AMPAR-activated synaptic transmission, with only NMDAR-activated synaptic transmission, robust persistent states at low firing rates are possible (Fig. 5*C,D*). The bistable range increases nearly linearly with the NMDAR-activated conductance g_{NMDA}, whereas the lowest firing rate of the persistent state remains <40 Hz (Fig. 5*D, filled square*). The dramatically different input–output relations obtained with the AMPA- or NMDA-type synapses can be explained in terms of their respective gating kinetics. As shown in Figure 5*C*(*top panel*), for a given synaptic coupling (g_{NMDA} = 0.006) and external drive (I = 0.3 nA), the input–output relation f(R) saturates at low firing rates with NMDAR channels, in contrast to the case with AMPAR channels (Fig. 5*A, top panel*). This is because the dependence of f on R is via the synaptic drive s_{E}(R) (Eq. 8). The fast-decaying AMPAR channels do not accumulate over time, hence do not saturate except at very high firing rates (∼500 Hz). By contrast, the slowly decaying NMDAR-mediated current saturates at firing rates within the physiological range (Fig.1*D*). At >50 Hz or so s_{NMDA}becomes independent of the input rate, so it can no longer be increased further to sustain higher firing rates. (The actual firing rate, which also depends on g_{NMDA} and the input I, can of course be >50 Hz.) For this reason, NMDAR (not AMPAR) channels are well suited to realize persistent states at low firing rates in a robust manner.

### Negative feedback mechanisms for rate control

Can some negative feedback mechanisms be used to resolve the problem of rate control with the AMPAR channels alone? This question is addressed next, by considering consecutively spike-frequency adaptation, recurrent shunting inhibition, and short-term synaptic depression.

#### Spike-frequency adaptation

Spike-frequency adaptation, a common property of (“regular spiking”) cortical pyramidal neurons (McCormick et al., 1985; Mason and Larkman, 1990; Wang, 1998), is added to the model neuron by including an I_{AHP}. To assess the effects of I_{AHP} on a persistent state sustained by the AMPAR channels, the frequency–current curve is calculated for different g_{AHP} values (Fig.6*A*). For a fixed I the firing rate of the active state is reduced by g_{AHP}, (Fig. 6*A, vertical dotted line*). At the same time, however, the bistable range shrinks dramatically and eventually disappears with large g_{AHP} values (for g_{AHP} ≥ 0.005).

This effect of I_{AHP} is readily explained in term of a negative current that counterbalances the excitatory synaptic current. Suppose that the firing rate is given by the input–output relation R = f(I_{tot}) (Eq. 15), where I_{tot} = I_{app} − I_{syn} − I_{AHP}. The average membrane potential of a firing neuron is approximately half-way between V_{reset} and V_{th}, V_{aυ} ≃ (V_{reset} + V_{th})/2 = −55.5 mV. Then, one has I_{syn} ≃ g_{AMPA}s_{E}V_{aυ} ≃ g_{AMPAυ}V_{aυ}R = −g̃_{AMPA}R, with s_{E} ≃ υR and g̃_{AMPA} = −g_{AMPAυ}V_{aυ}. On the other hand, I_{AHP} ≈ g_{AHP}[Ca^{2+}]_{aυ} (V_{aυ} − V_{K}) =g̃_{AHP}R, with [Ca^{2+}]_{aυ} = α_{Ca}τ_{Ca}R andg̃_{AHP} = g_{AHP}α_{Ca}τ_{Ca}(V_{aυ}− V_{K}). Taken together, we have
Equation 16Therefore, the addition of I_{AHP} amounts to a subtractive reduction of the effective recurrent synaptic excitation. For example, if g_{AMPA} = 1.2 and g_{AHP} = 0.0025, g̃_{AMPA} = −g_{AMPA}νV_{aυ} = 6.66 andg̃_{AHP} = g_{AHP}α_{Ca}τ_{Ca}(V_{aυ}− V_{K}) = 1.18. Thus,g̃_{AMPA} − g̃_{AHP} = 6.66 − 1.18 = 5.48. This is equivalent to a reduced g_{AMPA} value (g̃_{AMPA}− g̃_{AHP})/(−νV_{aυ}) = 0.99 in the absence of I_{AHP}. Indeed, the frequency–current curve with g_{AMPA} = 1.2 and g_{AHP} = 0.0025 and that with g_{AMPA} = 0.99 and g_{AHP} = 0 are essentially superimposable (Fig. 6*A*).

Note that the specific form of this subtraction depends on the model details. For example, if I_{AHP} has the following functional form I_{AHP} = g_{AHP}[Ca^{2+}]^{n}/([Ca^{2+}]^{n}+ D_{K}^{n})(V_{m} − V_{K}), n > 1, then the subtractive term in Equation 16 will be nonlinear.

It is important to emphasize that the stability of an active state is not guaranteed. In the it is shown that the stability of an asynchronous state depends critically on the synaptic time constant. In fact, with the fast AMPAR-mediated synapses, any active state in the presence of an I_{AHP} is expected to be unstable if its firing rate is below the lowest possible firing rate of an active state with I_{AHP} = 0. This is true regardless whether the active state belongs to a bistable range or not. For example, at I = 0.45 nA and g_{AHP} = 0.01 there is a single state with R = 30 Hz (Fig. 6*A*, cross). As shown in Figure 6*B*, this asynchronous state is not stable. Instead, neurons fire synchronously repetitive bursts of spikes that alternate with quiescent phases in time, the network oscillation has a frequency of 3 Hz. Such rhythmic bursting has also been reported in other studies that are not related to persistent activity (van Vreeswijk and Hansel, 1997; G. B. Ermentrout, personal communication). Synchronous burst oscillation is a common phenomenon in neurons and networks, usually when a strong and rapid autocatalytic process is combined with a slower negative feedback (here, the recurrent AMPAR-activated synaptic excitation and the I_{AHP}). Clearly, because the fast AMPAR-activated synaptic drive goes to zero between the bursts, the network would have to collapse onto the rest state, if the latter existed. In other words, when the active state in a bistable range is unstable (Fig. 6*A*, open circle), it is not observable, and the only stable behavior is the rest state. From these results it is concluded that I_{AHP} cannot subserve as a rate control mechanism unless additional slow synaptic transmission is present, such as that mediated by the NMDARs.

#### Recurrent shunting inhibition

Synaptic shunting inhibition has been suggested as a rate control mechanism in the neocortex (Douglas et al., 1995). When a neuron is at rest, shunting inhibition does not produce a net hyperpolarizing current because its reversal potential V_{I} is close to the resting potential. Instead, it causes an increase in membrane conductance, which divides the excitatory synaptic current (Carandini and Heeger, 1994). However, as it was recently pointed out by Holt and Koch (1997), the situation is different when the cell is in a repetitively firing state. In that case, the spiking mechanism essentially clamps the average membrane potential roughly half way between V_{reset} and V_{th}, well above V_{I}(for example, V_{I} = −70 mV, whereas V_{aυ} = (V_{reset} + V_{th})/2 = −55.5 mV), and the effect of inhibitory synapses is hyperpolarizing. For example, suppose that the model network is in a persistently active state, and each neuron receives a feedforward synaptic inhibition with a given input rate R_{I}. Then, this input is equivalent to a negative current I_{GABA} = g_{GABA}s_{I}(V_{aυ} − V_{I}), where the synaptic drive s_{I} as a function of R_{I} is given by Equation 10. Therefore, the addition of feedforward inhibition simply shifts a frequency–current curve to the right by the fixed amount I_{GABA}, without changing the range of network bistability or the lowest firing rate of a persistent state. This conclusion was confirmed by simulations (data not shown).

In the case of feedback synaptic inhibition, the firing of inhibitory interneurons is driven by pyramidal cells, and R_{I} is a function of R_{E}, R_{I} = g(R_{E}) (Fig.7*C*). In this case I_{b} remains the same, because g_{GABA} has no effect on the rest state. On the other hand, a larger I is needed to counterbalance I_{GABA} for the persistent activity (I_{GABA} shifts I_{a} to the right). Therefore the range of network bistability (I_{b} − I_{a}) is reduced. Note that, with increasing g_{GABA}, although the firing rate at a given I is reduced, the lowest possible rate of a persistent state (Fig. 7*A, filled square*) remains almost the same. Therefore, recurrent inhibition acts in a subtractive manner, in the sense that is produces a negative current that counterbalances the recurrent excitatory synaptic current. In terms of the firing rate equation R_{E} = f(I_{tot}), we have:
Equation 17with g̃_{GABA} = g_{GABA}(V_{aυ} − V_{I}) and R_{I} = g(R_{E}) (Eq. 14). The subtractive term is nonlinear in R_{E}.

If inhibitory neurons are not near the firing threshold, they will fire spikes only when their excitatory drive is sufficiently strong, e.g. R_{I} = 0 unless R_{E} is above a critical value ∼25 Hz (Fig. 7*C*). As a result, the portion of the frequency–current curve of the pyramidal cell with R_{E} < 25 Hz (on the middle branch) cannot be altered by feedback inhibition. With sufficiently strong g_{GABA}, network bistability is always preserved, and the lowest firing rate of a persistent state remains 25 Hz (Fig. 7*B*). In this way, persistent activity with reasonably low firing rates becomes possible.

I also considered an additional effect that may be caused by shunting inhibition. Suppose that shunting inhibition produces an increase in membrane conductance along a dendritic cable of length L, between the excitatory synapses and the spike triggering zone. The effective characteristic cable length λ is then expected to decrease like λ ∼ (g_{L} + g_{GABA}s_{I})^{−1/2}. To take into account the exponential attenuation of excitatory synaptic inputs along a passive cable, the excitatory conductance g_{E} should be multiplied by a factor ∼ exp(− L/λ) ∼ exp(− β(g_{GABA}s_{I})^{1/2}), where β is given in terms of the cable properties (Abbott, 1991). This highly nonlinear effect was suggested to provide a solution to the high firing rate problem in neural networks (Abbott, 1991). When this effect is included in the model, persistent states with low firing rates can be obtained, the frequency–current curve of the pyramidal cell is similar to Figure 7*B* (data not shown).

In any case, when a persistent state with low firing rate is realized with synaptic inhibition, its stability still remains to be determined. In fact, such a state was never observed in the network simulations, if the recurrent excitation was mediated exclusively by the fast AMPARs. Again, intuitively, such an active state is expected to be unstable because of the interplay between a fast recurrent excitation and a slower negative feedback. This is shown mathematically in the . To illustrate this point by computer simulations, I used the scaling parameter φ for the synaptic kinetics (Eqs. 4, 5) to change systematically the EPSC gating rates, whereas the average synaptic drive s_{E} and the firing rate R_{E} remained the same. φ was varied so that τ_{E} = τ_{s}/φ was between 2 and 80 msec. Let us choose g_{GABA} = 0.03 and I = 0.34 nA, the persistent state has a firing rate of 33 Hz (Fig. 7*B*). As shown in Figure8, when the excitatory synapses are slow (τ_{E} = 80 msec; comparable to that of the NMDAR channels), a persistent state can be sustained in the network (Fig. 8*A*). Because of the slow synaptic build-up, the network firing activity gradually ramps up during the input pulse. Moreover, in contrast to partially synchronous activity of Figure 4, with slow synaptic excitation (in the absence of a fast component) the persistent state is asynchronous. When τ_{E} is sufficiently reduced, the network activity in the persistent state displays increasingly large temporal fluctuations (Fig.8*B*). If τ_{E} is decreased below a critical value (τ_{E} ≃ 18 msec), the persistent state becomes unstable, because synchronous fluctuations eventually bring the network too close to the rest state, and the activity terminates (Fig. 8*C*).

To conclude, the effect of GABA_{A }synaptic inhibition is largely subtractive rather than divisive in repetitively firing neurons. Therefore, the phenomenon of persistent activity becomes less robust and can be abolished completely by strong recurrent inhibition. Moreover, when a persistent state with low rate does exist, it cannot be stably maintained unless the excitatory synapses are sufficiently slow (the ratio τ_{E}/τ_{I} must not be too small).

#### Short-term synaptic depression

I now turn to short-term depression of the excitatory synapses as a rate control mechanism. A typical simulation result is shown in Figure 9. In the absence of short-term depression (the parameter p_{υ} = 0; see Materials and Methods), a persistent activity state has a firing rate close to 200 Hz (Fig. 9*A*). The addition of short-term depression (p_{υ} = 0.35) reduces the firing rate to ∼40 Hz, back to the physiological range of PFC neurons (Fig.9*B*). Note that, because of short-term depression, the neuronal firing shows an exponential decrease during the depolarizing input pulse (Fig. 9*B*, top and middle panels); and immediately after the pulse there is a trough in the neural activity during which time the synapses recover from depression (Fig.9*B*, bottom panel). In this simulation both fast AMPAR and slow NMDAR channels are included, and the dynamics is asynchronous in the persistent state.

The frequency–current curve is calculated for different degrees of short-term depression (Fig.10*A*). In this case the undepressed AMPAR-mediated currents are so strong that with p_{υ} = 0 the firing rates of the persistent states are ∼500 Hz, near the neuronal saturation (data not shown). As we see in Figure 10*A*, short-term depression dramatically decreases the lowest firing rate of the active states (Fig. 10*A, filled square*). The range of bistability also shrinks (I_{a} shifts to the right) with increasingly strong short-term depression; but for some p_{υ} values this range remains reasonably large while the physiological firing rates are achieved. Short-term depression gives rise to synaptic saturation, which occurs at lower firing rates with larger p_{υ} (Fig.10*B*). Indeed, for AMPAR channels s_{E} ≃ νR_{E}. With short-term depression s_{E} = νR_{E}/(1 + p_{υ}τ_{D}R_{E}) (Eq. 9). In terms of the firing rate equation R = f(I_{tot}), we have:
Equation 18Therefore, the effect of short-term depression divides the amplitude of the excitatory synaptic drive. Unlike a subtractive mechanism (spike-frequency adaptation or recurrent inhibition), which is equally strong at all rates, a divisive mechanism affects high rates disproportionally. This leads to the flattening of the f(I_{tot}(R)) curve (Fig.10*B*). At high frequencies [R_{E} ≫ 1/(p_{υ}τ_{D})], the synaptic current becomes independent of the firing rate (Abbott et al., 1997). As a result, the positive feedback between firing and synaptic excitation has to stop at some firing rate, well below the neuronal saturation level (∼500 Hz).

The dynamical stability of these asynchronous persistent states with short-term depression was checked in direct network simulations. Simulations were performed with both all-to-all and sparse couplings. In a sparse network, unlike an all-to-all network, the number of synaptic connections varies widely from cell to cell, with an average M_{syn}. One might expect that such heterogeneity would favor an asynchronous persistent state against instability and synchrony. In fact, I found that as long as M_{syn}is not too small (≥100), the network behaves similarly with sparse or all-to-all coupling. This is true independent of the network size N_{e}. In other words, what matters is the absolute number of connections per neuron M_{syn}, not the connection probability p = M_{syn}/N_{e}. Similar to the case of spike-frequency adaptation or recurrent inhibition, it was found that fast AMPAR channels could not sustain such a low rate state, and that slower synapses were required (See for stability analysis). To be quantitative, for a given persistent state I varied the synaptic time constants systematically in network simulations by changing the scaling parameter φ (Eqs. 4, 5). This way, the smallest value of τ_{E} = τ_{s}/φ that was needed for the persistent state to be observable was determined. For example, consider the persistent states at I = 0.3 nA of Figure 10*A*, which have the firing rate ranged from 100 to 35 Hz as p_{υ} is varied from 0.15 to 0.35. The minimal τ_{E} required for the stability of each of these states is plotted as function of the firing rate R in Figure11*A*. The critical τ_{E} is larger with lower R, it also depends on the time constant of the depression process τ_{D} (see ).

Figure 11*A* was obtained with M_{syn} = 100 for a sparse network. The dependence on M_{syn} is shown in Figure11*B*, for R = 35 Hz. One observes that the minimal τ_{E} is not sensitive to M_{syn}, as long as M_{syn} ≥ 100. At very small M_{syn}, there is an abrupt increase of the required minimal τ_{E}, i.e. even slower synapses are needed to stabilize the active state. This is because, if a neuron receives a very small number of synaptic inputs, each at a low rate, the synaptic current must be long-lasting in order to produce a sustained tonic drive to the postsynaptic cell. Figure 11, *C*and *D*, illustrates the network dynamical behavior for τ_{E} around the critical minimum for R = 35 Hz (M_{syn} = 100). In these simulations, the network was initially set to be very close to the asynchronous state (see Materials and Methods for the asynchronous initial condition). Below the critical value (Fig. 11*C*; τ_{E} = 49 msec), the synchronous state is unstable. The network activity fluctuates in time, and R oscillates with growing amplitude. When R gets close to zero, the synaptic excitation becomes too weak to bring the network back up again, and the network activity dies out (Fig. 11*C*). On the other hand, above the critical value (Fig.11*D*; τ_{E} = 50 msec), fluctuations of the network activity decay with time, and the asynchronous persistent state is stable. In this random and sparse network, the firing rate of a neuron is a linear function of its number of synaptic connections (Fig. 11*D*, bottom panel), and is widely distributed across the neural population (20–60 Hz).

To conclude, unlike spike-frequency adaptation or synaptic inhibition, short-term depression acts as a divisive mechanism for rate control. The resulting persistent states at low firing rates are not stable, unless τ_{E} is larger than a critical value, which depends on both the short-term depression time constant and the firing rate. For the firing rates in the physiological range of PFC cells, the required synaptic kinetics is much slower than that of the AMPAR channels.

## DISCUSSION

The general finding of this work is that memory processes performed in strongly recurrent cortical circuits, such as delay-period activity, depend on the temporal dynamics as much as on the efficacy of recurrent synapses. Three main conclusions are: (1) the asynchronous dynamics is generally not stable in a fast recurrent excitation/slow negative feedback system; (2) slow NMDAR-activated synapses are powerful for maintaining a stable persistent activity at low firing rates; (3) short-term depression of excitatory synapses provides an efficient mechanism for rate control.

### NMDA receptors and persistent activity

NMDAR channels were found to be crucial to persistent activity in the network model for two reasons. First, their slow gating kinetics naturally leads to synaptic saturation at low firing rates, as observed experimentally (Fig. 1*A*), thereby contributing to the rate control of network activity. This saturation of the steady-state response to repetitive stimulation should be distinguished from receptor saturation by a single vesicle of transmitter; the latter is not supported by recent data (Mainen et al., 1999) and is not assumed in the present model. Second, slow synapses usually suppress network instability and oscillations, but are also able to sustain a partially synchronized network dynamics realized by other (fast) mechanisms. The voltage dependence of gating kinetics represents another interesting feature of I_{NMDA}to be explored in the context of working memory processes (Lisman et al., 1998).

Is there experimental evidence for a critical role of NMDARs in delay-period activity of the prefrontal cortex? Scherzer et al. (1998)reported a much higher expression of the NMDAR subunit mRNAs in the prefrontal cortex than in other cortical areas (such as primary visual cortex) of the human brain; which raises the interesting question of whether this regional difference could be correlated with the conspicuous occurrence of persistent activity in the association cortices in contrast to sensory cortices. NMDARs have been demonstrated to contribute to synaptic transmission at intracortical connections of sensory cortices (Thomson et al., 1985; Larson-Prior et al., 1991;Armstrong-James et al., 1993; Thomson and Deuchars, 1995; Markram et al., 1997) and frontal cortex (Sutor and Hablitz, 1989; Hirsch and Crépel, 1990; Kang, 1995). In certain cortical area, this contribution may overwhelm that of AMPARs and dominate recurrent horizontal excitations (Fleidervish et al., 1998). More quantitative analysis of the NMDAR- and AMPAR-mediated synaptic currents in the PFC has been lacking, for both the monkey and the rodent. On the other hand, in behavioral experiments with rats performing a spatial delayed alternation task, systematical administration (Verma and Moghaddam, 1996) or microinjection into the prefrontal cortex (Romanides et al., 1999) of NMDAR antagonists impaired working memory. These observations are consistent with our hypothesized importance of NMDARs to working memory. A direct experimental test, however, will need to be done on behaving animals, by combining pharmacological manipulation of NMDARs with neuronal recordings from the prefrontal cortex.

Note that in a model of persistent activity in the gaze control system,Seung (1996) also suggested that slow synaptic transmission is of crucial importance, but for quite different reasons. That network model is only weakly nonlinear, and slow synapses are useful to prolong the lifetime of transient memory storage.

### Rate control and robustness of network bistability

I have tested three candidate rate control mechanisms: spike-frequency adaptation, feedback inhibition, and synaptic short-term depression. I argue that a rate control mechanism should be assessed based on its effect on the entire frequency–current curve of the network. A rate control mechanism is judged effective if it reduces the lowest firing rate of persistent activity down to a physiologically plausible range; and at the same time the network bistability should remain robust within a reasonable parameter range. By these criteria, it was found that both spike-frequency adaptation and feedback inhibition are not adequate. Both act in a subtractive way, in the sense that each produces a negative current that counterbalances the recurrent excitatory synaptic current (Eqs. 16,17), and they readily abolish the persistent activity phenomenon. A note of caution is warranted there, because this study used the simple LIF neuron model that does not take into account more complex features of cortical neurons, such as dendritic morphology or other ionic currents that may contribute to single neuron dynamics. In particular, it would be worth re-examining the issue of feedback inhibition in a more realistic situation where, for example, shunting inhibition is located near the soma of a neuron, spatially separated from the excitatory inputs at dendritic sites. Moreover, our conclusion on shunting inhibition follows from the required preservation of network bistability, hence it does not deny the importance of recurrent inhibition as a rate control mechanism in situations without persistent activity, such as sensory processes in the primary visual cortex (Douglas et al., 1995;Borg-Graham et al., 1998). Finally, synaptic inhibition is likely indispensible for the formation of memory fields of the PFC neurons (Goldman-Rakic, 1995; Camperi and Wang, 1998; Rao et al., 1999).

In contrast to spike-frequency adaptation or synaptic inhibition, short-term synaptic depression acts as a divisive mechanism, in the sense that it divides the recurrent synaptic conductance (Eq. 18). Short-term depression reduces the firing rate not by preventing the neuronal saturation, but by saturating the synaptic drive at low firing rates (Fig. 10*B*). In the divisive but not subtractive case, firing rate of persistent activity is reduced effectively, whereas bistability is preserved in a robust way. Recent *in vitro* experiments have indicated that short-term depression is a general property of the rat PFC synapses (Hempel, Hartman, Wang, Turrigiano, and Nelson, unpublished observations). It would be interesting to see whether there is evidence for short-term depression in firing patterns of PFC cells of the behaving animal. Similar to our model simulation (Fig. 9*B*), in a delayed-response task, PFC neurons often display an exponential decrease of the firing rate during the cue presentation, followed by a trough of activity (Chafee and Goldman-Rakic, 1998; Romo et al., 1999; G. Rainer and E. K. Miller, personal communication). Such an effect needs to be measured quantitatively, and its underlying cellular mechanism remains to be elucidated.

### Stability and synchronization

To sustain persistent activity, a tonic synaptic drive is required to remain significantly above zero at any moment. This can be achieved by the fast AMPA-type synapses alone, if neuronal firings are asynchronous. However, previous work has shown that the asynchronous state is dynamically unstable if the excitatory synapses are too fast (Abbott and van Vreeswijk, 1993). I found that this problem is much more serious in the presence of a strong negative feedback mechanism for rate control. A pertinent question is to what extent this conclusion holds true in the presence of additional factors that increase the disorder of the network. Previous work has shown that noise has a stabilizing effect on the asynchronous dynamics of a network of excitatory neurons (Abbott and van Vreeswijk, 1993;Gerstner, 1999). In another study, a random network of excitatory and inhibitory neurons, coupled with instantaneous synapses, was found to be less synchronous with sparser connectivity (Brunel, 1999). None of these models contains a slow negative feedback mechanism. Here, in a network where recurrent excitation interacts with slow short-term depression, I found that asynchronous dynamics is not stable if the excitatory synapses are fast, even in the presence of synaptic noise and when the network connectivity is very sparse and the neuronal firing properties are widely heterogeneous (Fig. 11). Further analysis is needed to see if asynchronous dynamics are generally unstable in such fast recurrent excitation/slow negative feedback systems, even in the presence of heterogeneity and noise. The problem of stability of the asynchronous dynamics is of interest in the larger context of balanced excitatory–inhibitory neural networks (Shadlen and Newsome, 1994; van Vreeswijk and Sompolinski, 1996).

Therefore, a general finding here is that when an asynchronous persistent state has a low firing rate, its stability requires that the excitatory synaptic time constant be comparable to the effective time constant of the negative feedback mechanism. For a recurrent network of pyramidal cells and interneurons, the stability of a persistent state critically depends on whether the GABA_{A}R-mediated inhibition is as fast as the AMPAR-activated excitation. For both AMPARs (Geiger et al., 1995) and GABA_{A}Rs (Macdonald and Olsen, 1994), the deactivation kinetics is regulated by the subunit composition and thus may be specific for each cell type. In hippocampal pyramidal neurons of the rat, the decay time constant of the AMPAR-mediated EPSCs is ≃2 msec (at 35°C) (Hestrin et al., 1990a), whereas that of the fast component of the GABA_{A}R-mediated IPSCs is ≃6–10 msec (Banks et al., 1998). Hence IPSCs are approximately three to five times slower than EPSCs. The present study showed that such a mismatch of synaptic time constants does not favor the stability of an asynchronous dynamics at low firing rates, and for this reason the slow NMDAR channels could be required for the maintenance of a persistent state.

It is an open question whether completely asynchronous dynamics is indeed the *modus operandi* of delay-period activity in the PFC circuit. Funahashi (1998) recently reported that simultaneously recorded PFC cells displayed significant temporal correlations in a spatial working memory task. In my model simulations, when both the fast AMPA and slow NMDAR-mediated synaptic components are present, the fast AMPAR-activated recurrent excitation in interplay with slower negative feedback processes often leads to synchronous neural firings and network oscillations. In such a synchronous persistent state, the decay time constant of the slow synaptic component must not be too small compared to the average interspike interval (or oscillation period) of neurons. For typical firing rates of PFC cells of 10–50 Hz, ISI ≃20–100 msec, the NMDAR channels are needed.

### From cellular physiology to behavior

The present study raised and highlighted a number of experimental questions, their answers will contribute to bridge the gap between behavior-related neural activity and its underlying biological mechanisms.

#### Synaptic physiology of the prefrontal cortex

(1) What are the precise time courses of the AMPAR-mediated EPSCs and GABA_{A}R-mediated IPSCs? Is there a mismatch between the two? (2) In response to a repetitive train of stimuli, is the NMDAR-mediated EPSC a linear function of the stimulus frequency in the steady state? If not, what is the frequency above which the current saturates? (3) What are the relative amplitudes of the AMPAR- and NMDAR-mediated EPSCs? Can they be differentially modulated by neuromodulators such as dopamine (Cepeda et al., 1992)? (4) Across the somatodendritic membrane of a pyramidal neuron, is there a spatial segregation of excitatory and inhibitory synapses? (5) What are the short-term plasticity properties of PFC synapses?

#### Neural delay-period activity of the behaving animal

(6) Is there evidence for adaptation/depression of neuronal discharges? (7) How variable/random is the neuronal persistent activity? Do spike trains display some regular temporal structure? (8) Do neurons fire asynchronously, or is there synchronization within neural assemblies? (9) Can persistent activity of a neuron be switched off by an intracellularly injected current pulse? (no, if delay-period activity is network-induced; yes, if there is bistability at the single cell level) (10) Would local blockade of NMDA receptors in the PFC impair an animal's working memory performance? What are the correlated changes in the delay-period activity of PFC neurons?

### Implications for schizophrenia

In recent years, there is growing evidence that working memory impairments are prominent symptoms in schizophrenia (Goldman-Rakic, 1994; Weinberger and Berman, 1996), and that dysfunction of the NMDAR-mediated neurotransmission in the cortex may be at the origin of these cognitive deficits (Javitt and Zukin, 1991; Coyle, 1996). For example, a noncompetitive NMDA antagonist such as phencyclidine or ketamine produces working memory deficits in healthy human subjects that closely resemble schizophrenia (Javitt and Zukin, 1991; Krystal et al., 1994). Moreover, significant alternations in gene expression of the NMDA receptor subunits were found in PFC of schizophrenics (Akbarian et al., 1996). However, the cellular mechanisms through which working memory relies on the NMDAR channels are largely unknown. The present theoretical work suggests a candidate scenario for the working memory malfunction in PFC, namely, an imbalance between the fast AMPAR- and the slow NMDAR-mediated components of the recurrent synaptic transmission within the PFC circuit can give rise to network dynamical instability and disruption of delay-period persistent activity.

## Appendix

### Stability of an asynchronous state

In this Appendix, I show that, in general, an asynchronous persistent state is not stable in a fast recurrent excitation/slow negative feedback system. Using a heuristic approach, I will write a dynamical equation for the population activity, in each of the three cases: spike-frequency adaptation, synaptic shunting inhibition, and short-term synaptic depression. Then I will discuss in detail the stability analysis of such a dynamical system.

### General remark

Because the excitatory network has a large number of dynamical variables (at least as many as the number of pyramidal cells), a rigorous stability analysis of the network involves as many degrees of freedom (Abbott and van Vreeswijk, 1993; Treves, 1993; Gerstner, 1999). However, our approach is to focus on the fastest and most stable of all dynamical modes for the system (when decoupled from the negative feedback). This would yield a single dynamical equation for the population firing rate which, combined with another equation describing the negative feedback, forms a two-variable system. The idea is that if a steady state is not stable for the two-variable system, it must be unstable for the original network. On the other hand, if it is stable by this description, it still is not necessarily stable for the full network system.

### Spike-frequency adaptation

The starting point is the firing rate equation R = f(R,[Ca^{2+}]_{aυ}). f is the neuronal input–output relation, and the input current includes contributions from the recurrent excitatory synaptic current, which itself depends on R, the adaptation current that is proportional to [Ca^{2+}], as well as the external current I (not explicitly shown). Suppose that one can write a dynamical equation for R, like:
Equation 19where τ_{E} is a characteristic time constant for the most stable dynamical mode of the network (when g_{AHP} = 0). It is assumed to be dominated by the time constant of the excitatory synapses rather than the membrane time constant (Abbott and van Vreeswijk, 1993; Treves, 1993; Gerstner, 1999). The steady state is given by dR/dt = 0, from which R = f(R,[Ca^{2+}]_{aυ}) is recovered.

The equation for the calcium concentration averaged over a typical interspike interval 1/R is:
Equation 20Equations 19 and 20 constitute a dynamical system of two variables; it can be analyzed by the phase-plane technique (Strogatz, 1994; Rinzel and Ermentrout, 1998). It is convenient to choose R and g_{AHP}[Ca^{2+}] as independent variables, so that the function f(R, g_{AHP}[Ca^{2+}]_{aυ}) is the same for different adaptation strengths g_{AHP}. The nullcline dR/dt = F(R,[Ca^{2+}]_{aυ}) = 0 plotted on the phase plane has three branches (Fig.12). The first one is R = 0. The second one is a decreasing function of [Ca^{2+}]_{aυ}, signifying a reduction of R by adaptation. The third one is the middle branch connecting the other two branches. On the other hand, the nullcline d[Ca^{2+}]_{aυ}/dt = G(R,[Ca^{2+}]_{aυ}) = 0 is a straight line (Eq. 20). A steady state is given by an intersection of the two nullclines. For small g_{AHP}, there are three intersection points, one on each of the three branches. The third intersection point with the highest firing rate corresponds to the persistent state. With increasingly larger g_{AHP}, the R -nullcline remains the same (because it only depends on the product g_{AHP}[Ca^{2+}]_{aυ}), whereas the [Ca^{2+}]_{aυ}nullcline has a decreasing slope. With sufficiently large g_{AHP}, the persistent steady state is moved from the top branch to the middle branch of the R -nullcline (Fig. 12).

### Recurrent shunting inhibition

In the case of shunting inhibition, we can heuristically write two coupled equations for the firing rates R_{E} and R_{I} of the excitatory and inhibitory neurons,
Equation 21
Equation 22where τ_{E} and τ_{I}are the excitatory and inhibitory synaptic time constants, respectively. The steady states are given by dR_{E}/dt = 0 and dR_{I}/dt = 0, from which the neuronal input–output relations (Eq. 14) are recovered. Two examples with different inhibition strengths (g_{GABA}) are shown in Fig. 13, using the same parameters as in Figure 7*B*. As we see in Figure13*A*, the R_{E}-nullcline has three branches, and the asynchronous active state (R_{E}^{*},R_{I}^{*}) with a high rate is located on the top one. When the firing rate R_{E} is reduced sufficiently by strong recurrent inhibition, (R_{E}^{*},R_{I}^{*}) is moved to the middle branch of the R_{E}-nullcline (Fig. 13*B*).

### Short-term synaptic depression

In this case, the equations are:
Equation 23
Equation 24the steady states are given by dR/dt = 0 and dD/dt = 0, which yield R = f(R, D) and D = 1/(1 + p_{υ}Rτ_{D}). These nullclines are plotted in Fig. 14, for two different p_{υ} values. Again, the R -nullcline has three branches. The asynchronous persistent state with a high firing rate is on the top branch, and is moved to the middle branch when its rate is reduced by strong synaptic depression.

### Stability analysis

We have seen that in each of the three cases, the R_{E}-nullcline has three branches (the usual situation when there is a bistability between the rest state and an active state). The asynchronous active state, whose rate is reduced into the physiological range by a negative feedback mechanism, is typically located on the middle branch of the R -nullcline. I would like to show that this persistent state is not stable if the excitatory synapses are much faster than the effective time constant of the negative feedback. I will consider only the case of recurrent synaptic inhibition. The other two cases can be treated in a similar manner.

To study the local stability of a persistent state (R_{E}^{*}, R_{I}^{*}), in the presence of shunting inhibition, we shall linearize Equations 21 and 22. The local stability is determined by the matrix:
Equation 25evaluated at (R_{E}^{*}, R_{I}^{*}). This steady state is locally stable, if the two eigenvalues λ_{1} and λ_{2} of M have a negative real part. A qualitative analysis can be performed as follows. Let us define Tr = ∂F/∂R_{E} + ∂G/∂R_{I}, and Det = (∂F/∂R_{E})(∂G/∂R_{I}) − (∂F/∂R_{I}) (∂G/∂R_{E}), then λ_{1} and λ_{2} are solutions of the algebraic equation λ^{2} − Trλ + Det = 0; and λ_{1} + λ_{2} = Tr, λ_{1}λ_{2} = Det.

Because F decreases with R_{I} and G increases with R_{E}, we have:
Equation 26and the sign of ∂F/∂R_{E} depends on which branch of the R_{E}-nullcline is the steady state located. Moreover, we shall make use of the information about the slope of each nullcline at (R_{E}^{*}, R_{I}^{*}) (Fig. 13). For the R_{E}-nullcline,
Equation 27and for the R_{I}-nullcline,
Equation 28Suppose first that (R_{E}^{*}, R_{I}^{*}) is on the top branch (Fig. 13*A*), where the slope of the R_{E}-nullcline is negative,
R_{E}< 0. Combining Equation 27 with Equation 26, we have:
Equation 29From Equations 26 and 29, we deduce that λ_{1} + λ_{2} = Tr < 0 and λ_{1}λ_{2} = Det > 0. If λ_{1} and λ_{2} are real, clearly they must both be negative, because their sum is negative, and their product is positive. If λ_{1} and λ_{2} are complex, their real part is Tr/2 < 0. In both cases, we conclude that the steady state is stable.

Let us now assume that the steady state (R_{E}^{*}, R_{I}^{*}) is on the middle branch (Fig.13*B*), where the slope of the R_{E}-nullcline is positive,
R_{E}> 0. Thus,
Equation 30Because the slope of the R_{E}-nullcline is larger than that of the R_{I}-nullcline at (R_{E}^{*}, R_{I}^{*}), we have:
Equation 31This, combined with Equations 26 and 30, leads to Det = λ_{1}λ_{2} > 0.

Note that ∂F/∂R_{E}(> 0) and ∂G/∂R_{I}(< 0) are proportional to 1/τ_{E} and 1/τ_{I}; respectively. Therefore, the sign of Tr = ∂F/∂R_{E} + ∂G/∂R_{I} depends on the relative speeds of recurrent excitation and feedback inhibition. Suppose that τ_{E} is much smaller than τ_{I}, the positive term dominates and Tr = λ_{1} + λ_{2} > 0. Because λ_{1}λ_{2} > 0, if λ_{1} and λ_{2} are real, then they must be positive; if they are complex, then their real part is positive. In both cases, the asynchronous active state is unstable. If τ_{E} is much larger than τ_{I}, the negative term dominates, thus Tr = λ_{1} + λ_{2} < 0. Together with λ_{1}λ_{2} > 0, we conclude that the asynchronous active state is stable within this framework of the population-activity description.

## Footnotes

This work was supported by the National Science Foundation (IBN-9733006), the Alfred P. Sloan Foundation, and the W. M. Keck Foundation. I thank P. S. Goldman-Rakic, E. Marder, J. Lisman, and N. Brunel for discussions and helpful comments on this manuscript.

Correspondence should be addressed to Xiao-Jing Wang, Center for Complex Systems, Brandeis University, Waltham, MA 02454. E-mail:xjwang{at}volen.brandeis.edu.