Journal of Neurophysiology

Temporal Integration by Stochastic Recurrent Network Dynamics With Bimodal Neurons

Hiroshi Okamoto, Yoshikazu Isomura, Masahiko Takada, Tomoki Fukai


Temporal integration of externally or internally driven information is required for a variety of cognitive processes. This computation is generally linked with graded rate changes in cortical neurons, which typically appear during a delay period of cognitive task in the prefrontal and other cortical areas. Here, we present a neural network model to produce graded (climbing or descending) neuronal activity. Model neurons are interconnected randomly by AMPA-receptor–mediated fast excitatory synapses and are subject to noisy background excitatory and inhibitory synaptic inputs. In each neuron, a prolonged afterdepolarizing potential follows every spike generation. Then, driven by an external input, the individual neurons display bimodal rate changes between a baseline state and an elevated firing state, with the latter being sustained by regenerated afterdepolarizing potentials. When the variance of background input and the uniform weight of recurrent synapses are adequately tuned, we show that stochastic noise and reverberating synaptic input organize these bimodal changes into a sequence that exhibits graded population activity with a nearly constant slope. To test the validity of the proposed mechanism, we analyzed the graded activity of anterior cingulate cortex neurons in monkeys performing delayed conditional Go/No-go discrimination tasks. The delay-period activities of cingulate neurons exhibited bimodal activity patterns and trial-to-trial variability that are similar to those predicted by the proposed model.


Graded neuronal activity that steadily increases or decreases during a delay period is typically seen in the prefrontal cortex when animals are anticipating sensory cues (Rainer et al. 1999), motor responses (Hoshi et al. 2000), or rewards (Watanabe 1996). In addition, such a graded activity may provide neural substrates for decision making (Barraclough et al. 2004; Brody et al. 2003; Constantinidis et al. 2002; Hanes and Schall 1996; Hernández et al. 2002; Lo and Wang 2006; Mazurek et al. 2003; Shadlen and Newsome 2001; Takeda and Funahashi 2002), invariant recognition of stimulus sequence (Hopfiled and Brody 2001), and perception of interval timing (Durstewitz 2003; Kitano et al. 2003; Leon and Shadlen 2003; Reutimann et al. 2004). It is widely accepted that delay-period activity of prefrontal neurons may be generated by reverberating synaptic input (Goldman-Rakic 1995) and the graded activity was modeled by recurrent networks with slow synapses (Wang 2002), special wiring patterns, or neuron-dependent activation thresholds (Miller et al. 2003; Seung et al. 2000), as well as by single-neuron mechanisms (Durstewitz 2003). The mechanism underlying graded neuronal activity or, more generally, temporal integration of an external input has been the focus of extensive investigation.

Here, we propose a recurrent neural network of stochastic neurons as a model for neural mechanisms for temporal integration. In the network, spike-triggered afterdepolarizing potentials enable the model neurons to exhibit distinct firing states, a baseline firing state, and an elevated firing state in the presence of an external input. Then, driven by this input, the neurons show transitions between these two states in a temporally organized manner. If we appropriately tune the weight of uniform recurrent synapses and the variance of background input, the number of neurons in the elevated firing state changes slowly at an external input-dependent speed to yield a graded population activity. The present temporal integration mechanism differs from the previously proposed deterministic mechanisms using bistable neurons (Koulakov et al. 2002) or bistable dendritic compartments (Goldman et al. 2003; Loewenstein and Sompolinsky 2003; Teramae and Fukai 2005). The present stochastic model does not assume any mechanism that relies on differences between individual processing units, such as neuron-dependent activation threshold. All neurons, and even recurrent synapses, can be identical. From trial to trial, stochastic noise recruits neurons sequentially for graded activity in a randomized order.

We examined the validity of the proposed mechanism by analyzing the graded activity recorded previously from the anterior cingulate cortex (ACC) of monkeys performing delayed Go/No-go discrimination tasks (Isomura et al. 2003). The ACC, which engages in motor (Isomura and Takada 2004; Shima and Tanji 1998), cognitive (Botvinick et al. 2004; Brown and Braver 2005; Greene et al. 2004; Shidara and Richmond 2002), and emotional (Bishop et al. 2004) functions, receives major inputs from the prefrontal and limbic cortices (Dum and Strick 1993). These inputs might be integrated during the Go/No-go discrimination tasks. The graded neuronal activity in ACC displayed bimodal frequency distributions and large trial-to-trial variability analogous to those predicted by the present model, suggesting that it has a biological reality.


Recurrent network model

Typically, we used a network of 500 excitatory and 100 inhibitory neurons. Each excitatory neuron projects to 10% of randomly chosen other excitatory neurons and to all inhibitory neurons, whereas each inhibitory neuron projects to all excitatory neurons, but not to other inhibitory neurons. We subsequently show the mathematical details of the present network model. The values of parameters are listed in a separate subsection.

The membrane potential dynamics of an excitatory neuron is described by leaky integrate-and-fire neuron models as C(dV/dt) = −GL(VVL) − ICATIAMPAIGABA + Iext + IBG, where C is the membrane capacitance, and GL and VL represent the conductance and reversal potential of leak current, respectively. When V reaches threshold Vθ, the neuron generates a spike. Then, V is reset to Vreset and is clamped at this voltage during a refractory period. Each postsynaptic spike triggers Ca2+ entry into the cell body and raises the concentration of intracellular Ca2+. The time course of the Ca2+ concentration is given as d[Ca2+]/dt = −[Ca2+]/τCa + ΔCaδ(tts), where ts denotes the time at which the postsynaptic neuron fires. Then, Ca2+-dependent cationic current generates a depolarizing potential that follows the spike generation. The current is modeled as ICAT = GCAT(VVCAT)[Ca2+]4/([Ca2+]4 + KCAT4) (Aoyagi et al. 2002). We note that the Ca2+-dependent cation current is known to be abundant in cortical neurons (Haj-Dahmane and Andrade 1997; Okada et al. 1999).

α-Amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA)–receptor-mediated recurrent synaptic input is described as IAMPA = GAMPA(VVAMPA) ∑n pn, where the sum is taken over presynaptic neurons. When the nth presynaptic neuron fires, the gating variable pn is instantaneously increased by rp(1 − pn) with rp = 0.5; otherwise, it obeys dpn/dt = −pnAMPA. We describe γ-aminobutyric acid type A (GABAA)–receptor-mediated synaptic input from inhibitory neurons in a similar way: IGABA = GGABA(VVGABA) ∑n qn, dqn/dt = −qnGABA, and qnqn + rq(1 − qn) at each presynaptic spike.

Iext represents an external input, which may represent internally or externally driven information that should be integrated during a delay period. Iext is depolarizing during climbing activity, whereas Iext = 0 during descending or spontaneous activity. To terminate a climbing activity, we apply a transient hyperpolarizing current at the end of the delay period (the downward arrow in Fig. 3A, top). To induce a descending activity, we apply a transient depolarizing current that lifts the neuronal activity to a prerequisite peak level at the beginning of the delay period (the upward arrow in Fig. 4A, top).

FIG. 1.

Bimodal firing states of model excitatory neuron. A: responses of a single excitatory neuron to a brief stimulus are shown in the absence of recurrent synaptic inputs and the fluctuating components of background synaptic inputs (the “frozen” condition). External input was set as Iext = 0 nA (top), for which the response was not bistable, and Iext = 0.025 nA (bottom), for which the response was bistable. In the latter case, neuronal firing was terminated by a hyperpolarizing input. Horizontal bars show the duration of the stimuli. B: model neuron with bistability repeats noise-driven transitions between the baseline and elevated firing states under the influences of continuous synaptic bombardments (top). Monitoring the intracellular calcium density enables us to distinguish the epochs of the elevated firing state (bottom, gray shades). C: presence of the 2 distinct firing states results in a bimodal consecutive firing-rate distribution. D: bimodal firing-rate distribution is shown at Iext = 0.035 nA.

FIG. 2.

Comparison between different temporal integration mechanisms. A: graded activity may be modeled as a trial- or an ensemble-average of gradually increasing firing rates of individual neurons. B: climbing activity (bottom) was constructed from nonstationary Poisson spike trains with a gradually increasing mean firing rate (top). C: graded activity in our model consists of temporally organized bimodal transitions between the baseline and elevated firing states. In the individual neurons, the transitions should occur at arbitrary temporal positions with equal probabilities. Both trial average and ensemble average give equally good representations of graded activity in the present model. D: climbing activity (bottom) was constructed from artificial bimodal Poisson spike trains showing stepwise increases in the mean firing rate (top). E: consecutive firing-rate distribution (see methods) exhibits a single peak in the climbing activity shown in B. F: by contrast, the firing-rate distribution is bimodal in the climbing activity shown in D.

FIG. 3.

Climbing activity in a recurrent network model. A: driven by a constant depolarizing input (Iext = 0.02 nA; top), a model neuron showed bimodal state transitions at various temporal positions in repeated trials of a simulated delay period (middle). Each transition to the elevated firing state could be detected by the Ca2+ transients (gray shades). Arrow indicates a brief hyperpolarizing input to terminate climbing activity. Averaging the single-neuron activity over 50 trials gives climbing activity (bottom). B: similarly, averaging spike trains over different neurons gives climbing activity in a single trial. C: consecutive firing-rate distributions of the single-cell activity shown in A were bimodal in 3 different periods. D: a similar bimodal consecutive rate distribution was obtained from the population neural activity shown in B. E: serial orders of activation in 50 trials are shown for randomly chosen 5 neurons. F: duration of the elevated firing state in each neuron was averaged over 50 trials. Individual neurons exhibit almost the same average duration, implying that the temporal order of their activations was determined randomly in different trials.

FIG. 4.

Descending activity in a recurrent network model. A and B: descending activity can be induced by a brief depolarizing input given at the beginning of a delay period (arrow). Here, Iext = 0. Model neurons showed bimodal state transitions at various temporal positions (gray shades in spike raster) during a simulated delay period. Both trial-averaged single-cell activity and single-trial population neural activity give similar descending profiles. C and D: both trial- and ensemble-averaged activities produce bimodal consecutive firing-rate distributions. E: serial orders of inactivation in 50 trials are shown for randomly chosen 5 neurons. F: average duration of the elevated firing state was almost the same in the individual neurons.

IBG represents the net background synaptic input from external neuron pools: IBG = μ + IE + II, where μ is a constant and IE and II represent excitatory and inhibitory inputs mediated by AMPA and GABAA receptors, respectively. They are described as IE = −GEsE(VVAMPA) and II = −GIsI(VVGABA). The gating variables obey dsE,I/dt = −sE,IAMPA,GABA + XE,I(t), with XE and XI representing independent Poisson processes of the rates RE and RI, respectively. The constant bias μ is introduced for the convenience in adjusting the firing rate of spontaneous activity in the baseline state (typically ∼1–3 Hz), although, in principle, we could do this without μ by adjusting the values of GE, GI, RE, and RI. The bistability of the neuron model was assessed under the “frozen” condition, in which IAMPA = IGABA = 0 and IBG is replaced with its average, that is, the fluctuating components are eliminated: IE + II → −GEREτAMPA(VVAMPA) − GIRIτGABA(VVGABA).

The membrane potential dynamics of an inhibitory neuron are described as C(dV/dt) = −GL(VVL) − IAMPA + IBG, where IAMPA stands for the synaptic input from excitatory neurons and IBG for the background synaptic input modeled similarly in the case of excitatory neurons. We include the inhibitory neurons to make the present network model realistic. However, these neurons do not play a crucial role in the proposed mechanism of temporal integration.

Parameter values

In the excitatory neuron model, C = 0.5 nF, GL = 0.025 μS, VL = −70 mV, Vθ = −52 mV, and Vreset = −62 mV. The refractory period was 4 ms. The AMPA-receptor–mediated synaptic current was modeled with τAMPA = 5 ms, VAMPA = 0 mV, and rp = 0.5. A suitable range of the maximum conductance of recurrent synapses depends on the network size. Unless otherwise stated, we set as GAMPA = 0.02 nS at excitatory-to-excitatory as well as excitatory-to-inhibitory AMPA synapses. We modeled the Ca2+-dependent cationic current with GCAT = 6 nS, VCAT = −35 mV, and KCAT = 1 μM (Aoyagi et al. 2002; Kang et al. 1998), and the intracellular calcium dynamics with τCa = 70 ms and ΔCa = 1.3 μM (Wang 1999). The background synaptic input to each excitatory neuron was modeled with GE = GI = 0.7 nS and RE = RI = 1.8 kHz. The constant bias was adjusted as μ = 0.2 nA.

In the inhibitory neuron model, C = 0.2 nF, GL = 0.02 μS, VL = −65 mV, Vθ = −52 mV, and Vreset = −60 mV. The refractory period was 2 ms. In GABAA-receptor–mediated inhibitory-to-excitatory synapses, τGABA = 10 ms, VGABA = −80 mV, rq = 0.5, and GGABA = 0.02 nS. The background synaptic input to an inhibitory neuron was modeled with GE = GI = 0.35 nS and RE = RI = 1.8 kHz. The constant bias was adjusted as μ = 0.5 nA.

Consecutive firing-rate distribution

We examined whether graded neuronal activity might consist of bimodal rate changes by using the consecutive firing-rate distribution. We first chose a target period to analyze the spike trains of a single neuron. The distribution visualized how the firing rate of a neuron was distributed during the target period, which did not necessarily involve sufficiently many spikes. In this study, the typical interval of the target period was 500 to 600 ms. We divided the target period into consecutive 200- or 300-ms-long intervals, allowing a 100- or 150-ms overlap between the neighboring intervals, respectively. Then, in each trial we calculated the firing rates averaged over each of these intervals. For instance, if the target period is 500 ms long and the interval is 200 ms long, the period can be divided into four consecutive intervals consisting of 0–200, 100–300, 200–400, and 300–500 ms. If spike trains of a given neuron are available for 20 trials, we may obtain (4 × 20 =) 80 samples of the firing rate. Constructing a normalized distribution of these samples of firing rate gives the consecutive firing-rate distribution.

The temporal length of each interval sets the resolution of the rate changes detectable by the consecutive firing-rate distributions. If the resolution should be <5 Hz, the interval should be >200 ms. We must, however, make the interval sufficiently small compared with a delay period, to see how the rate distributions evolve during that period. Within these conflicting demands, an interval of 200–300 ms usually yielded excellent results.

The consecutive rate distribution played a major role in the statistical tests of the present study. The distribution of the instantaneous firing rate (i.e., the inverse interspike intervals) is given as ρλ(ν) = (λ/ν2) exp(−λ/ν) for a Poisson spike train of rate λ. Then, we can derive the following rate distribution for a nonstationary Poisson spike train with the mean rate increasing steadily from λmin to λmax Math The corresponding consecutive firing-rate distribution should approximately coincide with this. This formula demonstrates that the rate distribution would have a single peak near at a frequency of (λmax − λmin)/2 (ln λmax − ln λmin)2 (see Fig. 2D), if spike trains exhibit graded rate changes in each trial. This fact was used in the statistical tests to discriminate the two hypotheses proposed for graded activity of ACC neurons.

Statistical test

In experiments, the consecutive firing-rate distributions of climbing or descending activity exhibited prominent troughs, which possibly imply that the neurons underwent discontinuous rate changes (see Fig. 2, C and F). Indeed, these troughs cannot be consistent with the hypothesis that the graded activity constitutes truly graded activities in the individual trials (see Fig. 2, A and E). Rather, the troughs represent statistically significant deviations from the predictions of this hypothesis. To prove this, we calculated the distribution expected for the firing rate ν, hypothesizing that graded activity represented the average of single-trial Poisson spike trains with such a graded rate change as indicated by the peristimulus time histogram (PSTH). We divided the time axis of graded activity into B consecutive bins of width Δ, which was typically 200 or 300 ms, and derived how many such bins might contain k = Δν spikes in each spike train. Let λb be the average firing rate in the bth bin of the PSTH. On the preceding hypothesis, the probability of finding k spikes in the bth bin of a single-trial spike train is πk,b = e−λbΔbΔ)k/k!. Then, the probability of finding nk bins in the spike train that contain k spikes is Pk(nk) = ∑s1, …, sB=0,1 pk,1(s1)pk,2(s2) … pk,B(sB), where pk,b(1) = πk,b and pk,b(0) = 1 − πk,b, and the sum is taken under the constraint s1 + s2 +… + sB = nk. By using the generating function gk(u) = ∑m=0B umPk(m) = ∏b=1Bk,bu + (1 − πk,b)], we can obtain the average and variance of nk as μk = ∂gk(u)/∂u|u=1 = ∑b=1B πk,b and σk2 = ∂2gk(u)/∂u2|u=1μk2 = ∑b=1B πk,b(1 − πk,b), respectively. We then plotted μk as a function of ν, which gives the firing-rate distribution for the observed graded activity on the hypothesis that the neuron exhibits a truly graded rate change in single trials (Fig. 6, GI, thick gray curves).

FIG. 5.

Parameter dependency of modeled graded activity. A1: temporal integration can be performed within suitable ranges of the intensity of background synaptic input and the synaptic weight: GAMPA = 0.02 nS and GE = GI = 0.7 nS (center). Inset: directions of positive and negative parameter changes with respect to the central panel. Value of GAMPA was increased (top) or reduced (bottom) by ±10%, or instead, the values of GE and GI were increased (right) or decreased (left) by ±5%. In each panel, we show the time evolution of the firing rate averaged over the excitatory neuron population for 3 different values of Iext: the larger the input, the steeper the slope of the rate increase. If the bias current is adjusted at 0.19, Iext may typically take its value within 0.01–0.04 nA in these examples. A2A9: value of GAMPA was increased (top row) or reduced (bottom row) by ±50%, and the value of GE (=GI) was increased (right column) or reduced (left column) by ±30% relative to the values used in A1. In these outmost panels, the sum of the bias and external currents typically takes a value within the range 0.14–0.29. B: slope of descending activity can be changed by the intensity of external input. At 1 s, the value of Iext was abruptly lowered from 0.012 nA to 0.008, 0.007, 0.006, 0.005, or 0.004 (the steepest curve).

FIG. 6.

Graded activity of monkey anterior cingulate cortex (ACC) neurons. A: ACC neuron showed climbing activity during the second delay period (2.8–4.0s) in the color Go trials (red-red). B: consecutive firing-rate distribution of the climbing activity manifested a bimodal shape with a trough at approximately 7 Hz. C: another ACC neuron displayed descending activity during the second delay period in the color Go trials (red-red). D: corresponding consecutive firing-rate distribution manifested a bimodal shape with a trough at approximately 5 Hz. E: this neuron exhibited descending activity during the second delay-period, which was preceded by a relatively high-frequency tonic firing. F: bimodal firing-rate distributions were obtained for the descending activity. GI: consecutive firing-rate distributions of the above 3 neurons were calculated in the entire periods of graded activity. Solid and dashed curves display the average and SE, respectively, of the bin counts predicted from the hypothesis that the graded activity consists of truly graded activity of single neurons. At the locations of the peaks (asterisks), the deviations from the predictions are statistically significant (P < 0.05). J: frequency distribution was calculated in the late epoch of the climbing activity shown in A. Peristimulus time histogram of this neuron exhibited an accelerated increase at times >3.4 s, which resulted in frequency components >20 Hz. K: rate distribution was shown in the early epoch of the descending activity shown in C. This neuron exhibited slightly steeper rate decreases at times <3 s, presumably being influenced by a cue stimulus. Rate distribution showed a rich spectrum of frequency components >40 Hz.

Thus obtained rate distribution in general has a single peak at some firing rate ν̄ (or some spike count = Δν̄). Therefore we performed statistical tests of the actual bin count at this firing rate. Assuming that the count obeys a Gaussian distribution with mean μ and variance σ2, we calculated the z-score to examine whether the difference between the actual and expected counts is statistically significant. Furthermore, we performed a t-test, in which the variance in the number of bins containing spikes was calculated from the trial-to-trail fluctuations around μ. This test enables us to examine the hypothesis of truly graded single-neuron activity allowing the trial-to-trial variations in the slope of graded rate changes. We note that the coefficient of variation (CV) values calculated from background or task-unrelated activities of several ACC neurons ranged from 0.73 to 1.02 with an average value of 0.91. This justifies the use of Poisson spike trains in the preceding statistical analysis.

Behavioral tasks and electrophysiological recordings

Two female Japanese monkeys (Macaca fuscata) were trained to perform self-paced, delayed conditional Go/No-go discrimination (Konorski-type) tasks using spatial (location) and nonspatial (color) visual cues, as previously described (Isomura et al. 2003). Briefly, the task started once the monkeys pressed a key and fixated on a fixation square on the CRT monitor. In the spatial-discrimination task, location-related visual cues were displayed (300 ms) on either the left or the right side of the fixation square twice at 1 and 2.5 s after the start. Subsequently, a go signal was displayed at the fixation position 4 s after the start. Thus each trial had two delay periods, 1.3–2.5 and 2.8–4.0 s, during which the visual cues were not present. If the two visual cues appeared in the same position, the monkeys had to release the key (Go trials) and if they appeared in different positions, the monkeys had to keep pressing the key until the fixation square disappeared 5 s after the start (No-go trials), to get a drop of juice as reward. In the nonspatial-discrimination task, color-related visual cues, blue and red squares, were displayed at the position of the fixation square instead of the location-related visual cues and the two cues in the same or different colors indicated the Go or No-go trials, respectively. The spatial- and nonspatial-discrimination tasks were alternately changed after every correct trial and therefore the monkeys were able to presume a task type (location or color) in advance. Go and No-go trials were systematically randomized within each of the spatial and nonspatial conditions. During the Go/No-go discrimination tasks, single-unit activity was recorded from the CMAr (area 24c), CMAd (area 6c), and CMAv (area 23c) on the side contralateral to the hand used for key pressing. All experiments were carried out in accordance with the Guide for the Care and Use of Laboratory Animals (National Institutes of Health 1996) and the Guideline for Care and Use of Animals (Tokyo Metropolitan Institute for Neuroscience 2000).


Bimodal response of single excitatory neurons

We first clarify the bimodal response properties of the single excitatory neurons. Figure 1A (top) shows the response of a single excitatory neuron to a brief depolarizing input in the absence of external and recurrent input. In the simulation, we eliminated the fluctuating components of the background synaptic input, keeping its mean (“frozen” condition: see methods). An action potential triggers Ca2+ entry through the voltage-dependent Ca2+ channels, which in turn activates Ca2+-dependent cation current to induce an afterdepolarizing potential. However, the afterdepolarizing potential is not large enough to produce regenerative spike discharge and thus the membrane potential returns the resting level (Fig. 1A, top). If we apply a weak constant external input Iext, the neuronal response turns to be bistable (Fig. 1A, bottom). That is, the depolarizing pulse induces a regenerative spike generation until it is terminated by a hyperpolarizing input. Because of this bistability, if we introduce the fluctuating components of background synaptic input, the neuron repeats distinct periods of baseline and elevated firing states (Fig. 1B, top). As illustrated by the shaded epochs ([Ca2+] > 2 μM), we can detect the transitions between the bimodal states by monitoring the concentration of intracellular Ca2+ (Fig. 1B, bottom). The consecutive firing-rate distribution (see methods) showed a peak at low frequencies (<10 Hz) and another broad peak at higher frequencies (20–30 Hz) (Fig. 1C). As Iext is increased, the peak at higher frequencies grows whereas that at lower frequencies decays (Fig. 1D), implying that the probability for the transition from the baseline to elevated firing state is increased with Iext. As we will see later, this property allows us to control the slop of graded activity by changing Iext.

We now propose a mechanism to obtain graded activity with such a binary neuronal response. Without this neuronal response property, graded activity may be obtained by averaging the gradual changes in firing rate of single neurons (Fig. 2A). The average may be taken across either different trials or different neurons, depending on how a specific model achieves the graded activity. Figure 2B displays the spike raster that was artificially generated by nonstationary Poisson spike trains with a mean firing rate increasing gradually from 3 to 30 Hz. Averaging these spike trains gives rise to a climbing activity, as shown by a PSTH. By contrast, graded activity can be obtained by averaging abrupt transitions between a baseline (3 Hz) and an elevated firing state (30 Hz) in a network of binary units (Fig. 2C). This mechanism works well when the occurrence of each transition is equally probable at an arbitrary time point during a delay period. The consecutive rate distributions of the different graded-activity types exhibited quite different profiles (Fig. 2, E and F). Most characteristically, the distribution obtained from the stepwise rate changes in single neurons exhibits a trough near the peak of the distribution obtained from the truly graded rate changes. Thus the rate distribution enables us to examine which type of graded activity given spike trains are more likely to represent.

Graded activity in recurrent neural networks

We then constructed a recurrent network consisting of 500 excitatory neurons and 100 inhibitory neurons (see methods). In the network model, excitatory neurons receive excitatory and inhibitory recurrent synaptic inputs, excitatory and inhibitory background synaptic inputs, and an external input to induce graded activity. Inhibitory neurons receive synaptic input from excitatory neurons as well as excitatory and inhibitory background synaptic inputs. Each excitatory neuron projects to 10% of randomly chosen other excitatory neurons and to all inhibitory neurons, whereas each inhibitory neuron projects to all excitatory neurons, but not to other inhibitory neurons. We note that the temporal integration performance was relatively independent of the connectivity of synapses. In the presence of the background synaptic inputs, the excitatory neurons in the network showed a low- (∼5 Hz) and a higher-frequency (∼10–30 Hz) firing state, respectively. Below, we demonstrate that our network model actually exhibits the type of graded activity proposed in Fig. 2C.

Depending on stimulus protocols, excitatory model neurons show climbing or descending activity in the simulations of the network. A constant depolarizing input applied to these neurons induced climbing activity, which was obtained either by averaging spike trains of a single neuron over repeated trials (Fig. 3A) or by averaging spike trains of a neuron ensemble in a single trial (Fig. 3B). In each trial, climbing activity was terminated by a brief inhibitory input. We constructed the consecutive firing-rate distributions of the climbing activity from the spike trains obtained in repeated trials. In Fig. 3C, we calculated the distributions in three different periods to show the time evolution of the firing rate. The distributions resembled the rate distribution for stepwise-changing Poisson spike trains, but differed remarkably from that of gradually increasing Poisson spike trains (Fig. 2, E and F). The distributions, especially the one in an intermediate period of the simulated delay period, showed clear bimodal peaks. We could construct a similar bimodal profile of the rate distributions using spike trains of multiple neurons in a single simulation trial (Fig. 3D). Thus both trial- and ensemble-averaging procedures yielded similar climbing profiles, high variability in spike trains, and consecutive rate distributions. It is noted that the individual neurons significantly changed the order of activation from trial to trial (Fig. 3E). In fact, there was no a priori preferable order of activation among different neurons in the stochastic dynamics of the present uniform network, which consists of identical neurons interconnected by recurrent synapses of equal weights (Fig. 3F). This explains why spike trains were highly variable from trial to trial in the present model.

The network model exhibited descending activity if excitatory neurons were raised to the elevated firing state by a brief external input. As in climbing activity, averaging either spike trains of a single neuron over repeated trials (Fig. 4A) or spike trains of a neuron ensemble in a single trial generated descending activity (Fig. 4B). This activity terminated automatically when the activities of all excitatory neurons returned to the baseline level. The consecutive firing-rate distributions of descending activity showed clear bimodal peaks that indicate bimodal state transitions (Fig. 4, C and D). The temporal order of deactivation of the neurons changed significantly from trial to trial without a preferable order (Fig. 4, E and F).

The temporal profile of graded activity depended on several parameters in the model. In particular, the intensity of background input and the weight of recurrent synapses exerted significant influences on the network performance. Climbing activity grew at a nearly constant rate within adequate ranges of the values of these parameters and the slope of the climbing activity was modifiable by the intensity of a depolarizing external input. These properties were preserved if the noise intensity and the synaptic weight were tuned roughly within ±5 and ±10% of the optimal values, respectively (Fig. 5A1). However, the growth of population activity rapidly decelerated with time if recurrent connections were too weak (Fig. 5A2) or noise intensity was too strong (Fig. 5A3), or both (Fig. 5A6). By contrast, graded activity was accelerated excessively when the synaptic weight was too strong (Fig. 5A4) or the noise intensity was too weak (Fig. 5A5), or both (Fig. 5A8). The network model could tolerate the tuning errors if the two errors were correlated (Fig. 5, A7 and A9). These results imply that the climbing activity in this network model requires an adequate balance between the noise intensity and the synaptic weight. We can show similar results in the descending activity, where the intensity of a hyperpolarizing input modifies its slope (Fig. 5B).

Graded delay-period activity of ACC neurons

We tested the validity of the temporal integration with bimodal neuronal firing states using the data recorded from the anterior cingulate cortex. Anterior cingulate neurons of monkeys performing multitrial reward schedule tasks exhibited multimodal firing-rate distributions across trials and mixtures of a few Poisson distributions well represented a substantial proportion of the neuronal responses (Shidara et al. 2005). These results encouraged us to examine graded activity of ACC neurons. As subsequently demonstrated, results of our data analysis exhibited a surprising similarity to those of the present computational model.

The data set previously recorded from the monkey ACC contained 20 neurons that showed obvious climbing (n = 7) or descending (n = 13) activity during a delay period (Isomura et al. 2003). In this study, we primarily analyzed the activity of these neurons. Among the 20 cells, five cells showing climbing activity and 11 cells showing descending activity displayed the consecutive firing-rate distributions with clear bimodal peaks. Figure 6 shows typical examples of graded activity and their consecutive firing-rate distributions. The PSTHs representing the trial-averaged activity of single neurons display steady increases (Fig. 6A) or decreases (Fig. 6, C and E) during a delay period. In all examples, the trial-to-trial variability of temporal spiking patterns is quite significant during the period of graded rate changes. The rate distributions of each neuron exhibit bimodality, which is manifested by a trough between their low- and higher-frequency portions (Fig. 6, B,D,F). The higher-frequency portion grows with time in climbing activity, whereas it decreases gradually in descending activity.

We tested whether Poisson spike trains with gradually changing firing rates (i.e., the hypothesis illustrated in Fig. 2A) might consistently account for the consecutive firing-rate distributions of ACC neurons. To this end, we constructed the firing-rate distribution from the PSTH of graded activity of each ACC neuron. Examples of such rate distributions are displayed in Fig. 6, GI for the three ACC neurons shown previously. If the rate change in a single trial has a profile similar to that of the PSTH—that is, if the above hypothesis is true—the consecutive firing-rate distribution will exhibit a profile with a single peak. Therefore we examined whether the actual bin count represents a statistically significant deviation from the count obtained from the PSTH at the location of the peak (see methods). We note that the peaks were generally located near the troughs of the consecutive firing-rate distributions. Results of our statistical tests confirmed the statistical significance of the deviations. In total, 16 cells exhibited clear bimodal rate distributions. Among them, the deviations were statistically significant in all five cells showing climbing activity in both z-score and t-test (P < 0.05). The remaining 11 cells displayed descending activity and the troughs were statistically significant in nine or ten cells in z-score or t-test, respectively (P < 0.05). Thus, graded activity of ACC neurons gives spike trains, PSTHs, and the firing-rate distributions that are consistent with those obtained by simulations of the present network model.

We noticed that climbing activity of ACC neurons often displayed modulations of higher-frequency components at the late stage of a delay period. Figure 6J shows such an example for the climbing activity presented in Fig. 6A. The frequency components >20 Hz reflected the accelerated rate increase at times >3.4 s in the PSTH of this neuron. Similarly, Fig. 6K displays the firing-rate distribution during the initial epoch of the descending activity shown in Fig. 6C. The frequency components >40 Hz originated from the epoch of a steeper rate decrease at times <3 s in the PSTH. These results may suggest that the rate of the elevated firing state undergoes continuous modulations or that the ACC neurons might have more than two modes of firing states.


Slow temporal integration by stochastic bimodal neurons

Graded neuronal activity, which gradually increases or decreases typically during a delay period of cognitive task, is generally considered to represent temporal integration of information relevant to behaviors. We have proposed a novel form of temporal integration in a uniform recurrent network of stochastic bimodal neurons. In this network, stochastic neural dynamics recruits an ensemble of neurons for integrating an external input in single trials. The mechanism does not require N-methyl-d-aspartate (NMDA)–receptor-mediated slow recurrent synapses because the stochastic transitions between the bimodal neuronal states are slow enough to integrate input in a range of hundred milliseconds to seconds. A formal mathematical theory for this stochastic mechanism was developed in our previous paper (Okamoto and Fukai 2001; Sakai et al. 2006). Interestingly, the proposed mechanism gives no prescribed temporal order of activation to each neuron.

Our analysis of previously recorded data (Isomura et al. 2003) revealed that graded delay-period activity of ACC neurons exhibits bimodal neuronal responses analogous to those predicted by the network model. It was reported that activity of ACC neurons displayed bimodal firing-rate distributions across different trials of cognitive task (Shidara et al. 2005). Bimodal or, more generally, multimodal firing state might be the method by which ACC neurons encode information. However, the experimental data available for the present analysis were limited, so the results need to be confirmed by further experiments. By contrast, graded activity might constitute truly graded activity of individual neurons, the hypothesis used in several previous models of graded activity (Durstewitz 2003; Miller et al. 2003; Wang 2002). Prefrontal neurons display graded delay-period activity in decision making (Barraclough et al. 2004; Constantinidis et al. 2002; Takeda and Funahashi 2002), in which an internal signal representing subject's confidence level may start growing until it reaches a certain criterion for decision (Hanes and Shall 1996; Maimon and Assad 2006; Reddi and Carpenter 2000; Usher and McClelland 2001). The prefrontal cortex of the monkey performing a parametric working memory task (Brody et al. 2003) and the posterior parietal cortex of the monkey performing a duration discrimination task (Leon and Shadlen 2003) also exhibit graded activity. It seems intriguing to investigate whether and how the mechanism of temporal integration may differ in different cortical regions using the multiunit recording technique.

Physiological origin of bimodal neuronal responses

Our model used the spike afterdepolarizing potential induced by voltage-dependent Ca2+ current and Ca2+-dependent cation current to generate bimodal neuronal responses. Several other mechanisms can generate a similar bimodal response (Gruber et al. 2003; Loewenstein et al. 2005). For instance, the NMDA-receptor–mediated synaptic current provides a mechanism of bistability (Lisman et al. 1998). Spontaneous membrane potential fluctuations might be another possible source of bimodality in cortical neurons. In vivo and in vitro cortical neurons display membrane potential fluctuations between a depolarizing UP and a resting DOWN state (Anderson et al. 2000; Bazhenov et al. 2002; Cossart et al. 2003; Ikegaya et al. 2004; Peterson et al. 2003; Shu et al. 2003; Steriade et al. 2001; Stern et al. 1997). Although these fluctuations exhibit irregular oscillating patterns in anesthetized animals, they might temporally be organized in behaving states. However, whether the two-state fluctuations appear in an awake or a behaving animal remains controversial (Kitano et al. 2002; Peterson et al. 2003; Steriade et al. 2001).

Relation to previous models

Several models of neural integrator networks were previously proposed in the literature (Cannon et al. 1983; Galiana and Outerbridge 1984) for temporal integration. Seung et al. (2000) and Miller et al. (2003) showed that a special class of finely tuned recurrent networks produces a continuous attractor state useful for temporal integration. Later, this fine-tuning was relaxed by using bistable neurons or bistable dendritic compartments (Goldman et al. 2003; Koulakov et al. 2002). Because of the specific wiring structures or the neuron-specific activity thresholds, these models will activate neurons in an approximately fixed order with small trial-to-trial jitters arising from biological noise. By contrast, our model has essentially no preferable order of activation among neurons. We also note that our model does not require the perfect stability of bimodal states because stochastic noise destabilizes such states.

Models with slow synaptic transmissions also produce gradually climbing or descending neuronal activity (Brody et al. 2003; Wang 2002). Because spike trains in these models are highly variable, these models may display spike raster similar to that of our model. However, the inspection of firing-rate distributions will enable us to discriminate the two types of network models (Fig. 1). Single neurons can also perform temporal integration of stimuli (Durstewitz 2003; Egorov et al. 2002; Fransen et al. 2006; Loewenstein and Sompolinsky 2003; Teramae and Fukai 2005). Depending on the detailed cellular mechanisms, such neurons presumably produce a continuous or a multi-stepwise-graded rate change.

Mongillo et al. (2003) proposed another type of recurrent neural network model that may be consistent with the results of the present analysis of graded activity in the ACC. They considered two-state transitions in synaptic activity triggered by stochastic noise. Unlike the present model, their model generates bimodal responses of the entire network through recurrent excitation, so all neurons in the network exhibit coherent two-state transitions. Nevertheless, the trial average of single-cell activity displays a graded profile because the time of the coherent transition changes at random across trials. In addition, the model neurons exhibit bimodal rate distributions similar to those analyzed in this study. Although the model by Mongillo et al. does not serve as a neural integrator in single trials, the bimodal firing-rate distributions revealed in the ACC neurons do not exclude their model. Further discriminations between their model and ours require multiunit recordings and single-trial–basis analyses of graded activity.


The present work was partially supported by Scientific Research on Priority Areas—Higher-Order Brain Functions Grant-in-Aid 17022036 and The Ministry of Education, Culture, Sports, Science and Technology Grants-in-Aid for Scientific Research 16300096 and 16500190.


We thank S. Amari and Y. Sakai for fruitful discussion of the statistical tests in experimental data analysis. We thank T. Takekawa and S. Kang for technical assistance.


  • The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.


View Abstract