Statistical-Mechanical Measure of Stochastic Spiking Coherence in A Population of Inhibitory Subthreshold Neurons
aa r X i v : . [ phy s i c s . b i o - ph ] O c t Noname manuscript No. (will be inserted by the editor)
Statistical-Mechanical Measure of Stochastic SpikingCoherence in A Population of Inhibitory SubthresholdNeurons
Woochang Lim · Sang-Yoon Kim
Received: date / Accepted: date
Abstract
By varying the noise intensity, we study stochastic spiking coherence (i.e., col-lective coherence between noise-induced neural spikings) in an inhibitory population ofsubthreshold neurons (which cannot fire spontaneously without noise). This stochastic spik-ing coherence may be well visualized in the raster plot of neural spikes. For a coherent case,partially-occupied “stripes” (composed of spikes and indicating collective coherence) areformed in the raster plot. This partial occupation occurs due to “stochastic spike skipping”which is well shown in the multi-peaked interspike interval histogram. The main purposeof our work is to quantitatively measure the degree of stochastic spiking coherence seen inthe raster plot. We introduce a new spike-based coherence measure M s by considering theoccupation pattern and the pacing pattern of spikes in the stripes. In particular, the pacingdegree between spikes is determined in a statistical-mechanical way by quantifying the aver-age contribution of (microscopic) individual spikes to the (macroscopic) ensemble-averagedglobal potential. This “statistical-mechanical” measure M s is in contrast to the conventionalmeasures such as the “thermodynamic” order parameter (which concerns the time-averagedfluctuations of the macroscopic global potential), the “microscopic” correlation-based mea-sure (based on the cross-correlation between the microscopic individual potentials), and themeasures of precise spike timing (based on the peri-stimulus time histogram). In terms of M s , we quantitatively characterize the stochastic spiking coherence, and find that M s re-flects the degree of collective spiking coherence seen in the raster plot very well. Hence,the “statistical-mechanical” spike-based measure M s may be used usefully to quantify thedegree of stochastic spiking coherence in a statistical-mechanical way. Keywords
Inhibitory Subthreshold Neurons · Stochastic Spiking Coherence · Statistical-Mechanical Measure
Recently, brain rhythms have attracted much attention (Buzs´aki 2006). Synchronous oscil-lations in neural systems may be used for efficient sensory processing (e.g., binding of the
Woochang Lim · Sang-Yoon Kim (corresponding author)Department of Physics, Kangwon National University, Chunchon, Kangwon-Do 200-701, KoreaE-mail: [email protected] integrated whole image in the visual cortex is accomplished via synchronization of neuralfirings) (Gray 1994). In addition to such neural encoding of sensory stimuli, neural syn-chronization is also correlated with pathological rhythms associated with neural diseases(e.g., epileptic seizures and tremors in the Parkinson’s disease) (Hauschildt et al. 2006;Grosse et al. 2002; Tass et al. 1998). Here, we are interested in these synchronized neu-ral oscillations.A neural circuit in the major parts of the brain such as thalamus, hippocampus, andcortex consists of a few types of excitatory principal cells and diverse types of inhibitory in-terneurons. Functional diversity of interneurons increases the computational power of prin-cipal cells (Buzs´aki 2006; Buzs´aki et al. 2004). To understand the mechanisms of syn-chronous brain rhythms, neural systems composed of excitatory neurons and/or inhibitoryneurons have been investigated, and thus three types of synchronization mechanisms havebeen found (Wang 2003): recurrent excitation between principal cells, mutual inhibitionbetween interneurons, and feedback between excitatory and inhibitory neurons. Perfect syn-chronization occurs in the network of pulse-coupled excitatory neurons (Kuramoto 1991;Mirollo & Strogatz 1990). However, this synchronization via mutual excitation cannot bealways stable in networks with slowly decaying synaptic couplings. Stability of synchro-nization is shown to depend on both the time course of synaptic interaction and the responseof neurons to small depolarization (Vreeswijk et al. 1994; Hansel et al. 1995). When thedecay time of the synaptic interaction is enough long, mutual inhibition (rather than exci-tation) may synchronize neural firing activities. By providing a coherent oscillatory outputto the principal cells, interneuronal networks play the role of the backbones (i.e., synchro-nizers or pacemakers) of many brain oscillations such as the 10-Hz thalamocortical spin-dle rhythms (Wang & Rinzel 1992; Golomb & Rinzel 1994) and the 40-Hz fast gammarhythms in the hippocampus and the neocortex (Wang & Buzs´aki 1996; White et al. 1998;Whittington et al. 2000; Tiesinga et al. 2001). When the feedback between the excitatoryand inhibitory populations is strong, neural synchrony appears via the “cross-talk” betweenthe two populations (Whittington et al. 2000; Tiesinga et al. 2001; Hansel & Mato 2003;B¨orgers & Kopell 2003, 2005).Most past studies exploring mechanisms of neural synchronization were done in neuralsystems composed of spontaneously firing (i.e., self-oscillating) suprathreshold neurons. Forthis case, neural coherence occurs via cooperation of regular firings of suprathreshold self-firing neurons. In contrast, neural systems composed of subthreshold neurons have receivedlittle attention. Unlike the suprathreshold case, each subthreshold neuron in the absence ofcoupling cannot fire spontaneously without noise; it can fire only with the help of noise.Recently, stochastic spiking coherence (i.e., collective coherence emerging via cooperationof noise-induced spikings) was observed in an excitatory population of pulse-coupled sub-threshold neurons (Wang et al. 2000; Lim & Kim 2007, 2009). This kind of works may bethought to correspond to a “subthreshold version” of neural synchronization through mutualexcitation. Due to the stochastic spiking coherence, synaptic current, injected into each in-dividual neuron, becomes temporally coherent. Hence, temporal coherence resonance of anindividual subthreshold neuron in the network may be enhanced. This enhancement of co-herence resonance in an excitatory network of subthreshold Hodgkin-Huxley neurons wascharacterized in terms of the coherence factor b , representing the degree of sharpness ofthe peak in the power spectrum of an individual neuron (Wang et al. 2000). In this way, themeasure b is used for characterization of the temporal coherence of an individual neuron.However, we note that b is not a measure for directly measuring the degree of collectivecoherence in the whole population. In this paper, we are concerned about the “subthreshold version” of neural synchro-nization via mutual inhibition. In Section 2, we describe the biological conductance-basedMorris-Lecar (ML) neuron model with voltage-gated ion channels (Morris & Lecar 1981;Rinzel & Ermentrout 1998; Tsumoto et al. 2006). The ML neurons (used in our study)exhibit the type-II excitability (i.e., the firing frequency begins to increase from a non-zero value when the stimulus exceeds a threshold value), and they interact via inhibitoryGABAergic synapses whose activity increases fast and decays slowly. In Section 3, wecharacterize stochastic spiking coherence in a large population of inhibitory subthresholdML neurons by varying the noise intensity for a fixed coupling strength. Weakly coherentstates with oscillating ensemble-averaged global membrane potential V G are thus found toappear in a range of intermediate noise intensity [i.e., regular global oscillation (with re-duced amplitude and increased frequency) emerges via cooperation of irregular individualoscillations]. Emergence of collective coherence may be well described in terms of the con-ventional “thermodynamic” order parameter which concerns the time-averaged fluctuationsof the macroscopic global potential V G . We note that this stochastic spiking coherence maybe well visualized in the raster plot of neural spikes (i.e., a spatiotemporal plot of neuralspikes) which is directly obtained in experiments. For the coherent case, “stripes” (com-posed of spikes and indicating collective coherence) are found to be formed in the rasterplot. Due to coherent contribution of spikes, local maxima of the global potential V G appearat the centers of stripes. However, these stripes are partially occupied. Individual inhibitoryneurons exhibit intermittent spikings phase-locked to V G at random multiples of the period of V G . This “stochastic phase locking” leading to “stochastic spike skipping” is well shown inthe interspike interval (ISI) histogram with multiple peaks. The multi-peaked ISI histogramshows some indication of weak collective spiking coherence.The main purpose of our work is to quantitatively measure the degree of stochastic spik-ing coherence seen in the raster plot. In Section 3, we introduce a new type of spike-basedcoherence measure M s by taking into consideration the occupation pattern and the pacingpattern of spikes in the stripes of the raster plot. In particular, the pacing degree betweenspikes is determined in a statistical-mechanical way by quantifying the average contributionof (microscopic) individual spikes to the (macroscopic) global potential V G . This “statistical-mechanical” measure M s is in contrast to the conventional measures such as the “thermo-dynamic” order parameter (Golomb & Rinzel 1994; Hansel & Mato 2003), the “micro-scopic” correlation-based measure (based on the cross-correlations between the microscopicindividual potentials) (Wang & Buzs´aki 1996; White et al. 1998), and the measures of pre-cise spike timing based on the peri-stimulus time histogram (PSTH) (Mainen & Sejnowski1995; Schreiber et al. 2003). The “thermodynamic” order parameter and the “microscopic”measure concern just the the macroscopic global potential V G and the microscopic individualpotentials, respectively without considering any quantitative relation between V G and the mi-croscopic individual potentials. (The auto-correlation of the global activity used in the workof Brunel and Hakim (1999) may also be regarded as a kind of “thermodynamic” measure.)For the PSTH-based measure “events,” corresponding to peaks of the instantaneous popula-tion firing rate, are selected through setting a threshold. Then, the measures for the reliabilityand the precision of spike timing concern only the spikes within the events, in contrast to thecase of the “statistical-mechanical” measure where all spikes are considered (without select-ing events). A main difference between the conventional and the new spike-based measureslies in determining the pacing degree of spikes. The precision of spike timing for the con-ventional case is given by just the standard deviation of (microscopic) individual spike timeswithin an event without considering the quantitative contribution of (microscopic) individ-ual spikes to the (macroscopic) global activity. Hence, the PSTH-based measure is not a statistical-mechanical measure. However, if we take the instantaneous population firing rateas a global activity and exactly define the global cycles [see Fig. 4(a)] and the global phases[see Eqs. (15) and (16)] like our case, then the conventional PSTH-based measure mayalso develop into a similar statistical-mechanical measure. By varying the noise intensity,we quantitatively characterize the stochastic spiking coherence in terms of the “statistical-mechanical” measure M s , and find that M s reflects the degree of collective spiking coherenceseen in the raster plot very well. We also expect that M s may be implemented for character-izing the degree of collective coherence in the experimentally-obtained raster plot of neuralspikes. Finally, a summary is given in Section 4. In this section, we describe the biological neuron model used in our computational study.We consider an ensemble of N globally coupled neurons. As an element in our neural sys-tem, we choose the conductance-based ML neuron model with voltage-gated ion channels,originally proposed to describe the time-evolution pattern of the membrane potential forthe giant muscle fibers of barnacles (Morris & Lecar 1981; Rinzel & Ermentrout 1998;Tsumoto et al. 2006). The population dynamics in this neural network is governed by a setof the following differential equations: C dv i dt = − I ion , i + I DC + D x i − I syn , i , (1) dw i dt = f ( w ¥ ( v i ) − w i ) t R ( v i ) , (2) ds i dt = a s ¥ ( v i )( − s i ) − b s i , i = , · · · , N , (3)where I ion , i = I Ca , i + I K , i + I L , i (4) = g Ca m ¥ ( v i )( v i − V Ca ) + g K w i ( v i − V K ) + g L ( v i − V L ) , (5) I syn , i = JN − N (cid:229) j ( = i ) s j ( t )( v i − V syn ) , (6) m ¥ ( v ) = . [ + tanh { ( v − V ) / V } ] , (7) w ¥ ( v ) = . [ + tanh { ( v − V ) / V } ] , (8) t R ( v ) = / cosh { ( v − V ) / ( V ) } , (9) s ¥ ( v i ) = / [ + e − ( v i − v ∗ ) / d ] . (10)Here, the state of the i th neuron at a time t (measured in units of ms) is characterized bythree state variables: the membrane potential v i (measured in units of mV), the slow recov-ery variable w i representing the activation of the K + current (i.e., the fraction of open K + channels), and the synaptic gate variable s i denoting the fraction of open synaptic ion chan-nels. In Eq. (1), C represents the capacitance of the membrane of each neuron, and the timeevolution of v i is governed by four kinds of source currents.The total ionic current I ion , i of the i th neuron consists of the calcium current I Ca , i , thepotassium current I K , i , and the leakage current I L , i . Each ionic current obeys Ohm’s law. Theconstants g Ca , g K , and g L are the maximum conductances for the ion and leakage channels, and the constants V Ca , V K , and V L are the reversal potentials at which each current is bal-anced by the ionic concentration difference across the membrane. Since the calcium current I Ca , i changes much faster than the potassium current I K , i , the gate variable m i for the Ca + channel is assumed to always take its saturation value m ¥ ( v i ) . On the other hand, the activa-tion variable w i for the K + channel approaches its saturation value w ¥ ( v i ) with a relaxationtime t R ( v i ) / f , where t R has a dimension of ms and f is a (dimensionless) temperature-liketime scale factor.Each ML neuron is also stimulated by the common DC current I DC and an indepen-dent Gaussian white noise x [see the 2nd and 3rd terms in Eq. (1)] satisfying h x i ( t ) i = h x i ( t ) x j ( t ′ ) i = d i j d ( t − t ′ ) , where h· · ·i denotes the ensemble average. The noise x is a para-metric one which randomly perturbs the strength of the applied current I DC , and its intensityis controlled by the parameter D . The last term in Eq. (1) represents the coupling of thenetwork. Each neuron is connected to all the other ones through global synaptic couplings. I syn , i of Eq. (6) represents such synaptic current injected into the i th neuron. Here the cou-pling strength is controlled by the parameter J and V syn is the synaptic reversal potential. Weuse V syn = −
80 mV for the inhibitory synapse and V syn = s obeys the 1st order kinetics of Eq. (3) (Golomb & Rinzel 1994;Wang & Buzs´aki 1996). Here, the normalized concentration of synaptic transmitters, acti-vating the synapse, is assumed to be an instantaneous sigmoidal function of the membranepotential with a threshold v ∗ in Eq. (10), where we set v ∗ = d = v is largerthan v ∗ ). The synaptic channel opening rate, corresponding to the inverse of the synaptic risetime t r , is a =
10 ms − , which assures a fast rise of I syn (B¨orgers & Kopell 2003, 2005). Onthe other hand, the synaptic closing rate b , which is the inverse of the synaptic decay time t d , depends on the type of the synaptic receptors. For the inhibitory GABAergic synapse(involving the GABA A receptors) we set b = . − , and for the excitatory glutamatesynapse (involving the AMPA receptors), we use b = . − (B¨orgers & Kopell 2003,2005). Thus, I syn decays slowly.The ML neuron may exhibit either type-I or type-II excitability, depending on the sys-tem parameters (Rinzel & Ermentrout 1998). For the case of type-I (type-II), the firing fre-quency begins to increase from zero (a non-zero finite value) when I DC passes a threshold.Throughout this paper, we consider the case of type-II excitability where g Ca = . / cm , g K = / cm , g L = / cm , V Ca =
120 mV , V K = −
84 mV , V L = −
60 mV , C = m F / cm , f = . , V = − . , V =
18 mV , V = , and V =
30 mV. As I DC passes a thresh-old in the absence of noise, each single type-II ML neuron begins to fire with a nonzero fre-quency that is relatively insensitive to the change in I DC (Hodgkin 1948; Izhikevich 2000).Numerical integration of Eqs. (1-3) is done using the Heun method (San Miguel & Toral2000) (with the time step D t = .
01 ms) similar to the second-order Runge-Kutta method,and data for ( v i , w i , s i ) ( i = , . . ., N ) are obtained with the sampling time interval D t = [ v i ( ) , w i ( ) , s i ( )] for the i th ( i = , . . ., N ) neuron with uniform probability in the range of v i ( ) ∈ ( − , ) , w i ( ) ∈ ( . , . ) , and s i ( ) ∈ ( . , . ) . (a) t (ms) v ( m V ) (b) ISI (ms) P r obab ili t y (c) t (ms) v ( m V ) Fig. 1
Intermittent noise-induced firings of a single subthreshold ML neuron: (a) time series of the membranepotential v and (b) interspike interval (ISI) histogram for I DC = m A / cm and D = m A · ms / / cm ; theISI histogram is made of 5 × ISIs and the bin size for the histogram is 5 ms. Regular firings of a singlesuprathreshold ML neuron: (c) time series of v for I DC = m A / cm and D = In this section, we study stochastic spiking coherence in a population of inhibitory sub-threshold ML neurons. The stochastic spiking coherence is quantitatively characterized interms of a new “statistical-mechanical” spike-based measure.We first consider the case of a single ML neuron. An isolated subthreshold neuron cannotfire spontaneously without noise; it may generate firings only with the aid of noise. Figure1(a) shows a time series of the membrane potential v of a subthreshold neuron for I DC = m A / cm and D = m A · ms / / cm . Noise-induced spikings appear intermittently. Forthis subthreshold case, the ISI histogram is shown in Fig. 1(b). The most probable value ofthe ISIs (corresponding to the main highest peak) is 97 . v is in contrast to the regular oscillation of v for a suprathreshold neuron (whichfires spontaneously without noise) when I DC = m A / cm and D = N globally coupled subthreshold ML neuronsfor I DC = m A / cm and set the coupling strength as J = / cm . (Hereafter, for con-venience we omit the dimensions of I DC , D , and J .) By varying the noise intensity D , weinvestigate the stochastic spiking coherence. Emergence of collective spiking coherence maybe well described by the (population-averaged) global potential, V G ( t ) = N N (cid:229) i = v i ( t ) . (11)In the thermodynamic limit ( N → ¥ ) , a collective state becomes coherent if D V G ( t ) (= V G ( t ) − V G ( t )) is non-stationary (i.e., an oscillating global potential V G appears for a co-herent case), where the overbar represents the time average. Otherwise (i.e., when D V G isstationary), it becomes incoherent. Thus, the mean square deviation of the global potential (b) i (d) i t (ms) (c) V G ( m V ) t (ms) N s t (ms) N s t (ms) (e) V G ( m V ) (a) N=10 N=10 N=10 l og O log D Fig. 2
Order parameter and partial occupation in N (= ) globally coupled inhibitory subthreshold MLneurons: (a) plots of the order parameter O versus the noise intensity D . (b) partially-occupied raster plotof spikings a ( i : neuron index) and plot of number of spikes ( N s ) versus time ( t ) and (c) time series of theensemble-averaged global potential V G . Full occupation in N (= ) globally coupled excitatory ML neurons:(d) fully-occupied raster plot of spikings and plot of N s versus t and (e) time series of V G . For both inhibitoryand excitatory cases, I DC = m A / cm , D = m A · ms / / cm , J = / cm , and the bin size for N s in(b) and (d) is 1 ms V G (i.e., time-averaged fluctuations of V G ), O ≡ ( V G ( t ) − V G ( t )) , (12)plays the role of an order parameter used for describing the coherence-incoherence tran-sition (Manrubia et al. 2004). For the coherent (incoherent) state, the order parameter O approaches a non-zero (zero) limit value as N goes to the infinity. Figure 2(a) shows a plotof the order parameter versus the noise intensity. For D < D ∗ l ( ≃ . O tends to zero as N → ¥ . As D passes the lower threshold D ∗ l ,a coherent transition occurs because of a constructive role of noise to stimulate coherencebetween noise-induced spikings. However, for large D > D ∗ h ( ≃ .
4) such coherent statesdisappear (i.e., a transition to an incoherent state occurs when D passes the higher threshold D ∗ h ) due to a destructive role of noise to spoil the collective spiking coherence. As an example, we consider a coherent state for D =
20. For this case stochastic spikingcoherence may be well visualized in terms of the raster plot of neural spikes. The upperpanel in Fig. 2(b) shows the raster plot for N = . We note that stripes (composed of spikesand indicating collective coherence) appear successively in the raster plot. The number ofspikes ( N s ) in each stripe is given in the lower panel. About 100 spikes constitute a stripe[i.e., only a fraction (about 1/10) of the total neurons fire in each stripe]. In this way partialoccupation occurs in the stripes. A regularly oscillating global potential V G emerges throughcooperation of spikes in partially-occupied stripes. A time series of V G is shown in Fig. 2(c).The amplitude of V G is much smaller than that in Fig. 1(a) for the isolated single case.This reduction of amplitude occurs mainly because of partial occupation. Local maximaof V G appear at the centers of stripes where the spiking densities are locally highest. Theglobal period T G of V G (corresponding to the average intermax interval of V G or equivalentlycorresponding to the average interstripe interval in the raster plot) is 54.2 ms. Hence theglobal frequency f G (= . V G is roughly as twice as the most probable frequency (= . v for the isolated single case. Thus, a regular global oscillation with reducedamplitude and increased frequency occurs for the partially-occupied case. For comparison,we also consider the case of full occupation which occurs in a population of excitatoryneurons for the same set of parameters I DC , D , and J . [We emphasize that partial occupationmay also occur even for the excitatory case of smaller J (e.g., J = V G for the fully-occupied case, respectively.Fully-occupied stripes appear successively at nearly regular time intervals D t ( = . V G with large amplitudeappears via cooperation of spikes in the fully-occupied stripes, and its global frequency f G (= . f G for the fully-occupied case is nearly equal to the most probable frequency of v for the isolated singlecase).To further understand the full and the partial occupations in the raster plots, we examinethe local and the global output signals for both the fully-occupied and the partially-occupiedcases of D =
20. Since our neural network is globally coupled, any local neuron may bea representative one. Unlike the isolated single case (shown in Fig. 1), each local neuronwithin the network is stimulated by a synaptic current which is directly related to the aver-age of firing events of all the other neurons. First, we consider the fully-occupied excitatorycase, and investigate the time series of the local potential v of the 1st neuron and the globalpotential V G which are shown in Fig. 3(a). The 1st neuron fires spikes phase-locked to V G atevery cycle of V G (i.e., the local potential v exhibits a nearly regular phase-locked oscilla-tion). To confirm this 1:1 phase locking, we collect 5 × ISIs from all local neurons, andobtain the ISI histogram [see Fig. 3(b)]. A single peak with a small width is located aroundthe average value ( = . T G of V G . This 1:1phase locking is expected to occur for the fully-occupied case because the global frequency f G (= . V G is nearly equal to the most probable frequency (= . v of the 1st neuron and theglobal potential V G . The 1st neuron exhibits intermittent spikings phase-locked to V G at ran-dom multiples of the period T G (=54.2 ms) of V G . In addition to these intermittent spikingphases, hopping phases (exhibiting small subthreshold oscillations) appear in most of globalcycles. We also note that after occurrence of a spiking, recovery from a hyperpolarized toa resting state is made during the next global cycle. Hence, a “preparatory” phase withoutspiking and hopping (for preparing for generation of the next spike or hopping) follows t (ms) V G ( m V ) t (ms) V G ( m V ) -4030 (c) v ( m V ) -4030 (a) v ( m V ) (b) ISI (ms) P r obab ili t y (d) ISI (ms) P r obab ili t y Fig. 3 V G at every cycle of V G ) in N (= ) globally coupled excitatory ML neurons: (a) time series of the local potential v of the 1st neuron andthe ensemble-averaged global potential V G and (b) interspike interval (ISI) histogram with a single peak.Stochastic phase locking leading to stochastic spike skipping (intermittent spikings phase-locked to V G atrandom multiples of the period of V G ) in N (= ) globally coupled inhibitory subthreshold ML neurons: (c)time series of v and V G , and (d) ISI histogram with multiple peaks. For both inhibitory and excitatory cases, I DC = m A / cm , D = m A · ms / / cm , and J = / cm . Vertical dashed lines in (a) and (c) representthe times at which local minima of V G occur. Each ISI histogram in (b) and (d) is made of 5 × ISIs andthe bin size for the histogram is 5 ms. Vertical dotted lines in (d) denote integer multiples of the global period T G (=54.2 ms) of V G . each spiking phase (see the gray parts). In this way, the local potential exhibits an irregularfiring pattern consisting of randomly phase-locked spiking and hopping phases and prepara-tory phases. We note that a regular fast global oscillation with a small amplitude emergesfrom these irregular individual oscillations. To confirm the stochastic spike skipping (arisingfromn stochastic phase locking) in the local potential, we collect 5 × ISIs from all localneurons and get the ISI histogram. Multiple peaks appear at multiples of the period T G ofthe global potential V G . However, due to appearance of preparatory cycles, the 1st peak ofthe histogram (which is the highest one) appears at 2 T G (not T G ). Hence, local neurons firemostly in alternate global cycles. Due to this stochastic spike skipping partial occupationoccurs in the raster plot. Similar skipping phenomena of spikings (characterized with multi-peaked ISI histograms) were found in single noisy neuron models driven by a weak periodicexternal force for the stochastic resonance (Longtin 1995, 2000; Gammaitoni et al. 1998).Unlike this single case, stochastic spike skipping in networks of inhibitory subthresholdneurons is a collective effect because it occurs due to a driving by a coherent ensemble-averaged synaptic current. As discussed in details in Section 4, similar results on skippingwere also obtained in simplified networks of IF neurons (Brunel and Hakim 1999; Brunel i M i i M i (c) O i (d) O i G G (a) t (ms) V G ( m V ) - G G (b) t (ms) P i P i Fig. 4 “Statistical-mechanical” spike-based coherence measure in N (= ) globally coupled ML neuronsfor I DC = m A / cm , D = m A · ms / / cm , and J = / cm . Partially-occupied inhibitory case: (a)time series of the global potential V G , (b) plot of the global phase ( F ) versus time ( t ), and (c) plots of O i (occupation degree of spikes in the i th stripe), P i (pacing degree of spikes in the ith stripe), and M i (spikingmeasure in the ith stripe) versus i (stripe). In (a) and (b), vertical dashed and solid lines represent the timesat which local minima and maxima (denoted by open and solid circles) of V G occur, respectively and G i ( i = ,
2) denotes the ith global cycle. Fully-occupied excitatory case: (d) plots of O i , P i , and M i versus i . D =
20, we introduce a new“statistical-mechanical” spike-based coherence measure M s by considering the occupationpattern and the pacing pattern of spikes in the stripes of the raster plot. Particularly, thepacing degree between spikes is determined in a statistical-mechanical way by quantifyingthe average contribution of microscopic individual spikes to the macroscopic global poten-tial V G . The spiking coherence measure M i of the i th stripe is defined by the product of theoccupation degree O i of spikes (representing the density of the i th stripe) and the pacingdegree P i of spikes (denoting the smearing of the i th stripe): M i = O i · P i . (13)The occupation degree O i in the i th stripe is given by the fraction of spiking neurons: O i = N ( s ) i N , (14) where N ( s ) i is the number of spiking neurons in the i th stripe. For the full occupation, O i = O i <
1. The pacing degree P i of each microscopic spike inthe i th stripe can be determined in a statistical-mechanical way by taking into considerationits contribution to the macroscopic global potential V G . Figure 4(a) shows a time series ofthe global potential V G ; local maxima and minima are denoted by solid and open circles, re-spectively. Central maxima of V G between neighboring left and right minima of V G coincidewith centers of stripes in the raster plot. The global cycle starting from the left minimum of V G which appears first after the transient time (= ms) is regarded as the 1st one, whichis denoted by G . The 2nd global cycle G begins from the next following right minimumof G , and so on. Then, we introduce an instantaneous global phase F ( t ) of V G via linearinterpolation in the two successive subregions forming a global cycle (Freund et al. 2003).The global phase F ( t ) between the left minimum (corresponding to the beginning point ofthe i th global cycle) and the central maximum is given by: F ( t ) = p ( i − / ) + p t − t ( min ) i t ( max ) i − t ( min ) i ! for t ( min ) i ≤ t < t ( max ) i ( i = , , , . . . ) , (15)and F ( t ) between the central maximum and the right minimum (corresponding to the be-ginning point of the ( i + ) th cycle) is given by: F ( t ) = p ( i − ) + p t − t ( max ) i t ( min ) i + − t ( max ) i ! for t ( max ) i ≤ t < t ( min ) i + ( i = , , , . . . ) , (16)where t ( min ) i is the beginning time of the i th global cycle (i.e., the time at which the leftminimum of V G appears in the i th global cycle) and t ( max ) i is the time at which the maximumof V G appears in the i th global cycle. The global phase F in the first two global cyclesis shown in Fig. 4(b). Then, the contribution of the k th microscopic spike in the i th stripeoccurring at the time t ( s ) k to V G is given by cos F k , where F k is the global phase at the k th spiking time [i.e., F k ≡ F ( t ( s ) k ) ]. A microscopic spike makes the most constructive (in-phase) contribution to V G when the corresponding global phase F k is 2 p n ( n = , , , . . . ),while it makes the most destructive (anti-phase) contribution to V G when F i is 2 p ( n − / ) .By averaging the contributions of all microscopic spikes in the i th stripe to V G , we obtainthe pacing degree of spikes in the i th stripe, P i = S i S i (cid:229) k = cos F k , (17)where S i is the total number of microscopic spikes in the i th stripe. By averaging M i ofEq. (13) over a sufficiently large number N s of stripes, we obtain the ‘statistical-mechanical”spike-based coherence measure M s : M s = N s N s (cid:229) i = M i . (18)We follow 3 × stripes and get O i , P i , and M i in each i th stripe for the partially-occupied inhibitory case of Fig. 2(b). The results are shown in Fig. 4(c). For comparison,we also measure the degree of stochastic spiking coherence seen in the raster plot of Fig. 2(d)for the fully-occupied excitatory case and give the results in Fig. 4(d). We note a distinct dif-ference in the average occupation degree h O i i , where h· · ·i denotes the average over stripes. (e) M s log D (a1) i (a2) (a3) (a4) (b1) P r obab ili t y (b2) (b3) ISI (ms) (b4) (c) < O i > log D (d) < P i > log D N s t (ms) Fig. 5
Stochastic spiking coherence for various values of D in N (= ) globally coupled inhibitory sub-threshold ML neurons when I DC = m A / cm and J = / cm . (a) Raster plots of spikings ( i : neuron in-dex) and plots of number of spikes ( N s ) versus time ( t ) and (b) interspike interval (ISI) histograms for (a1) and(b1) D = m A · ms / / cm , (a2) and (b2) D = m A · ms / / cm , (a3) and (b3) D = m A · ms / / cm ,and (a4) and (b4) D = m A · ms / / cm . The bin size for N s in (a) is 1 ms. Each ISI histogram in (b) ismade of 5 × ISIs and the bin size for the histogram is 5 ms. Vertical dotted lines in (b) represent theinteger multiples of the global period T G of V G ; T G = (b1) 67.4 ms, (b2) 54.2 ms, (b3) 48.6 ms, and (b4) 47.7ms. (c) Plot of h O i i (average occupation degree of spikes) versus log D . (d) Plot of h P i i (average pacingdegree of spikes) versus log D . (e) Plot of M s (“statistical-mechanical” spike-based coherence measure)versus log D . To obtain h O i i , h P i i , and M s in (c)-(e), we follow the 3 × stripes for each D . Open circlesin (c)-(e) denote the data for D = m A · ms / / cm . For the inhibitory case, partial occupation with very small h O i i (= . ) occurs due tostochastic spike skipping, while for the excitatory case full occupation with h O i i = h P i i (= . ) is large in contrast to h O i i , although it is smaller than h P i i (= . ) for the fully-occupied case of excitatory coupling. As a result, the “statistical-mechanical” coherence measure M s (which represents the collective coherence seen in thewhole raster plot) is 0.081 for the partially-occupied inhibitory case which is very low whencompared with M s (= . ) for the fully-occupied excitatory case. The main reason for thelow degree of stochastic spiking coherence is mainly due to partial occupation.Finally, by varying D we characterize the stochastic spiking coherence in terms of the“statistical-mechanical” measure M s in a population of inhibitory subthreshold ML neuronswhen I DC =
87 and J =
3. Such stochastic spiking coherence may be well visualized in theraster plots of neural spikes. Figures 5(a1)-5(a4) show the raster plots and the temporal plotsof the number of spikes ( N s ) for N = in the coherent region for D =
10, 20, 30, and
32, respectively. The corresponding ISI histograms are also given in Figs. 5(b1)-5(b4). Wemeasure the degree of stochastic spiking coherence in terms of h O i i (average occupationdegree), h P i i (average pacing degree), and M s for 20 values of D in the coherent regime, andthe results are shown in Figs. 5(c)-5(e). (Here, for comparison with other values of D weinclude the case of D =
20 which is studied in details above; open circles in Figs. 5(c)-5(e)represent the data for D = h O i i = .
106 and h P i i = . n T G ( T G : global period of V G and n =2,3,...). As the value of D is increased from 20,the average occupation degree h O i i increases only a little, as might be seen from the rasterplots and the plots of N s in Figs. 5(a3) and 5(a4) ( h O i i = 0.114 and 0.115 for D =
30 and32, respectively). This slow increase in h O i i is well shown in Fig. 5(c). However,the aver-age pacing degree h P i i for D >
20 decreases rapidly, as shown in Fig. 5(d). For example,stripes in the raster plots of Figs. 5(a3) and 5(a4) become more and more smeared, andhence the average pacing degree is decreased with increase in D . This smearing of stripescan be understood from the change in the structure of the ISI histograms. As D is increased,the heights of peaks are decreased, but their widths are widened [see Figs. 5(b3) and 5(b4)].Thus, peaks begin to merge. This merging of peaks results in the smearing of stripes. Wenote that for large D merged multiple peaks in the ISI histogram are directly associated withsmeared partially-occupied stripes in the raster plot. Thus, for D >
20 the degree of stochas-tic spiking coherence is rapidly decreased as shown in Fig. 5(e), mainly due to the rapiddecrease in h P i i (a little increase in h O i i has negligible effect). Eventually, when passing thehigher threshold D ∗ h ( ≃ . ) stripes no longer exist due to complete smearing, and thenincoherent states appear.By decreasing the value of D from 20, we also characterize the stochastic spiking co-herence. As an example see the raster plot of spikes in Fig. 5(a1) for D =
10. The averageoccupation degree is much decreased to h O i i = . N s . This rapid decrease in h O i i for D <
20 can be seen in Fig. 5(c). Contrary to h O i i theaverage pacing degree h P i i (= . ) for D =
10 is decreased a little when compared tothe value of h P i i (= . ) for D =
20; only a little more smearing occurs. For this case,the ISI histogram has multiple peaks without merging like the case of D=20, as shown inFig. 5(b1). However, when compared to the case of D =
20 the heights of peaks decrease,and the average ISI increases via appearance of long ISIs. Thus, the degree of stochasticspiking coherence for D = ( M s = . ) is much decreased mainly due to rapid decreasein the average occupation degree h O i i (small decrease in h P i i has only a little effect). In fact,as can be seen in Fig. 5(d) there is no noticeable change in h P i i for 10 < D <
20; as D isdecreased from 20 to a value of D ( ≃ ) , h P i i increases slowly, and then it begins to de-crease. Thus, for D <
10 both h O i i and h P i i decreases so rapidly, as shown in Figs. 5(c) and5(d). Hence, as D is decreased from 10 the degree of stochastic spiking coherence decreasesrapidly due to both effects of h O i i and h P i i . Eventually, when D is decreased through thelower threshold D ∗ l ( ≃ . ) , completely scattered sparse spikes appear without forming anystripes in the raster plot, and thus incoherent states exist for D < D ∗ l . In this way, we char-acterize the stochastic spiking coherence in terms of the “statistical-mechanical” measure M s in the whole coherent region, and find that M s reflects the degree of collective coherenceseen in the raster plot very well. When taking into consideration both the occupation degreeand the pacing degree of spikes in the raster plot, a maximal spiking coherence occurs near D ≃
20 [see Fig. 5(e)]. As discussed in details in Section 1, M s is in contrast to the conven-tional measures where any quantitative relation between (microscopic) individual spikes and the (macroscopic) global potential V G is not considered [e.g., the “thermodynamic” orderparameter (Golomb & Rinzel 1994; Hansel & Mato 2003), the “microscopic” correlation-based measure (Wang & Buzs´aki 1996; White et al. 1998), and the PSTH-based measuresof precise spike timing (Mainen & Sejnowski 1995; Schreiber et al. 2003)]. By changing the noise intensity, we have characterized stochastic spiking coherence in an in-hibitory population of subthreshold biological ML neurons. In a range of intermediate noiseintensity, a regularly oscillating global potential V G (with reduced amplitude and the in-creased frequency) emerges via cooperation of irregular individual firing activities. We notethat this stochastic spiking coherence may be well visualized in the raster plot of spikes. Forthe case of a coherent state, partially-occupied stripes appear in the raster plot. This partialoccupation occurs due to stochastic spike skipping, which is shown clearly in the multi-peaked ISI histogram. Recently, Brunel and Hakim (1999) have obtained similar results onspike skipping in a simplified network of integrate-and-fire (IF) neurons which are randomlyconnected via instantaneous inhibitory synapses (modeled by the delta function with a trans-mission delay). [Such skipping was also reported in the network of inhibitory and excitatoryIF neurons (Brunel 2000; Brunel & Wang 2003).] For the case of the IF neuron a spike istriggered instantaneously when the membrane potential reaches a threshold, in contrast tothe case of the biological ML neuron where dynamics of voltage-gated ion channels leads tospike generation. As the stimulus exceeds a threshold, the firing frequency of the IF neuronbegins to increase from zero (i.e., the IF neuron exhibits the type-I excitability), unlike thecase of the type-II ML neuron (used in our computational study). These type-I IF neurons aredriven by random excitatory inputs activated by independent Poisson processes. In this sim-plified network, it was shown that weakly coherent global activities (roughly correspondingto the global potential V G ) emerge from irregular firing patterns of neurons [refer to Figure 3in the paper of Brunel and Hakim (1999)]. The collective coherence was well shown in thetemporal auto-correlation function of the global activity. Although the basic results in thesimplified network are similar to ours in a realistic network, there are some differences in thespiking pattern of individual neurons. The stochastic spike emission of individual IF neuronswas well shown in the Poissonian ISI histogram showing little indication of the collective be-havior. Since no clear multiple peaks appear in the Poisson-like ISI histogram, the individualIF neurons exhibit more stochastic spike emission than the ML neurons, although weak col-lective coherence occurs in both cases. There are some additional differences in the effect ofthe external noise on the collective coherence. As the magnitude of the external noise is in-creased from zero and passes a threshold, the global activity was found to exhibit a transitionfrom an oscillatory to a stationary state (Brunel and Hakim 1999). Thus, incoherent statesappear for the case of strong noise. On the other hand, through competition between theconstructive role (stimulating coherence between noise-induced spikes) and the destructiverole (spoiling the collective coherence) of noise, stochastic spiking coherence (in our work)appears in a range of intermediate noise intensities. Hence, incoherent states appear for bothcases of sufficiently strong and weak noise. These differences in both the spiking patternof individual neurons and the noise effect on the collective coherence seem to arise mainlyfrom different types of drivings on individual neurons; IF neurons are driven by random ex-citatory inputs activated by Poisson processes, while ML neurons are driven by a Gaussianwhite noise in addition to a subthreshold DC stimulus. The main purpose of our work isto quantitatively measure the degree of stochastic spiking coherence seen in the raster plot. We have introduced a new type of spike-based coherence measure M s by taking into con-sideration the occupation degree and the pacing degree of spikes in the stripes. Particularly,the pacing degree between spikes is determined in a statistical-mechanical way by quantify-ing the average contribution of (microscopic) individual spikes to the (macroscopic) globalpotential V G . This “statistical-mechanical” measure M s is in contrast to the conventionalmeasures (e.g., the “thermodynamic” order parameter, the “microscopic” correlation-basedmeasure, and the PSTH-based measures of precise spike timing) where any quantitative re-lation between the microscopic individual spikes and the macroscopic global potential V G is not considered. In terms of M s , we have quantitatively characterized the stochastic spik-ing coherence, and found that M s reflects the degree of collective spiking coherence seen inthe raster plot very well. Finally, we also expect that M s may be implemented to character-ize the degree of collective spiking coherence in an experimentally-obtained raster plot ofneural spikes. Acknowledgements
This research was supported by the Basic Science Research Program through the Na-tional Research Foundation of Korea funded by the Ministry of Education, Science and Technology (2009-0070865). S.-Y. Kim thanks Dr. D.-G. Hong for his interest in our work.
References
B¨orgers, C. & Kopell, N. (2003). Synchronization in networks of excitatory and inhibitory neurons withsparse, random connectivity.
Neural Computation, 15,
Neural Computation, 17,
Neural Computation, 11,
Journal of Computational Neuroscience, 8,
Journal of Neurophysiology, 90,
Trends in Neuroscience, 27,
Rhythms of the Brain . New York: Oxford University Press.Freund, J., Schimansky-Geier, L., & H¨anggi, P. (2003). Frequency and phase synchronization in stochasticsystems.
Chaos, 13,
Reviews of ModernPhysics, 70,
Physica D, 72,
Journal ofComputational Neuroscience, 1,
Clinical Neurophysiology, 113,
Neural Computation,7,
Neural Computation, 15,
Physical Review E, 74,
Journal of Physiology, 107,
International Journal of Bifurcation andChaos, 10,
Physica D,50,
Journalof the Korean Physical Society, 51,
International Journal of Modern Physics B, 23,
NuovoCimento D, 17,
Stochastic Dynamics and Pattern Formation in Biological and Complex Systems (pp.219-239). New York: AIP.Mainen, Z. F. & Sejnowski T. J. (1995). Reliability of spike timing in neocortical neurons.
Science, 268,
Emergence of Dynamical Order . Singapore:World Scientific.Mirollo, R. E. & Strogatz, S. H. (1990). Synchronization of pulse-coupled biological oscillators.
SIAM Jour-nal on Applied Mathematics, 50,
Biophysical Journal,35,
Methods in Neural Modeling: from Ions to Networks (pp. 251-292). Cambridge: MIT Press.San Miguel, M. & Toral, R. (2000). Stochastic effect in physical systems. In Martinez, J., Tiemann, R., &Tirapegui, E. (Eds.),
Instabilities and Nonequilibrium Structures VI (pp. 35-130) Dordrecht: KluwerAcademic Publisher.Schreiber, S., Fellous, J. M., Whitmer, D., Tiesinga, P., & Sejnowski, T. J. (2003), A new correlation-basedmeasure of spike timing reliability.
Neurocomputing, 52-54 , 925-931.Tass, P., Rosenblum, M. G., Weule, J., Kurths, J., Pikovsky, A., Volkmann, J., Schnitzler, A., & Freud, H.J. (1998). Detection of n : m phase locking from noisy data: application to magnetoencephalography. Physical Review Letters, 81,
Hippocampus, 11,
Neurocomputing, 69,
Journal of Computational Neuroscience, 1,
Neural Computation, 4,
The Journal of Neuroscience, 16,
Encyclopedia of Cognitive Science (pp. 272-280).London: MacMillan.Wang, Y., Chik, D. T. W., & Wang, Z. D. (2000). Coherence resonance and noise-induced synchronization inglobally coupled Hodgkin-Huxley neurons.
Physical Review E, 61,
Journal of Computational Neuroscience, 5,