Adaptation Reduces Variability of the Neuronal Population Code

Adaptation Reduces V ariabilit y of the Neuronal P opulation Co de 1 F ar zad F ar khooi, 1 , ∗ Eilif Muller , 2 and Martin P . Nawrot 1 2 1 Neur oi nformatics & The or etic al Neur oscienc e, F r eie Uni vers it¨ at Berlin and BCCN-Berlin, Germany 3 2 L ab or atory of Comput ational Neur oscienc e, EPFL, L ausanne, Switzerland 4 (Dated: August 8, 2018) 5 Sequences of ev ents in noise-driven excitable systems with slow v ariables often show serial correla- 6 tions among their in terv als of even ts. H ere, w e employ a master eq uation fo r generalized non-renewal 7 processes t o calculate th e interv al and count statistics of sup erimposed pro cesses go verned by a slow 8 adaptation v ariable. F or an ensemble of spike-frequency adapting neuron s, this results in the reg- 9 ularization of th e p opulation activity and an enhanced p ost-synaptic signal deco ding. W e confirm 10 our theoretical results in a p opulation of cortical neurons recorded in vivo . 11 P ACS n umbers: 87.19.ls, 05.40.-a, 87.19.lo, 87.19.lc 12 Statistical mo dels of ev ents assuming the renewal prop- 13 erty , tha t the insta n taneous probability for the o ccur- 14 rence of an even t depe nds uniquely on the time s ince 15 the last even t, enjo ys a long history of in terest and ap- 16 plications in ph ysic s. Howev er, many even t pro cesses in 17 nature violate the renew al prop erty . F or instance, it is 18 known that photon emissio n in multilev el quantum sys- 19 tems constitutes a non-r e ne wal pr o cess [1]. Likewise, the 20 time s e r ies of earthquakes typically exhibits a memo ry of 21 previous s hocks [2], as do the times of activ ated escape 22 from a meta s table state, as encountered in v ar io us sci- 23 ent ific fields such as chemical, bio logical, a nd solid state 24 ph ysic s [3]. Often, the departure from the renewal prop- 25 erty arises when the pr ocess under study is mo dulated 26 by so me slow v ariable, whic h r esults in ser ial correlatio ns 27 among the in terv als be tw een success iv e ev ents. In par- 28 ticular, the ma jority of spiking neurons in the nervous 29 systems of different species show a s erial dependence b e- 30 t ween inter-even t interv als (ISI) due to the fact that their 31 spiking activity is mo dulated by an in trins ic slow v ari- 32 able of s elf-inhibition, a pheno menon known a s spike- 33 frequency adaptation [4]. 34 In this letter, w e present a non-renewal formalism 35 based on a po pulation density tre a tmen t that enable s 36 us to quantitativ ely study ensemble pro cesses augmented 37 with a slow noise v aria ble. W e formally der iv e genera l ex- 38 pressions for the higher -order interv al and count statistics 39 of single and superimp osed non-renewal pro cesses for ar- 40 bitrary observ ation times. In spiking neurons, intrinsic 41 mechanisms of ada ptation reduce o utput v ariability and 42 facilitate p opulation co ding in neur al ensembles. W e con- 43 firm our theor etical r esults in a set o f ex perimental in vivo 44 recording s and a nalyse their implications for the read-o ut 45 prop erties o f a p ostsynaptic neur al deco der. 46 Non-r enewal Master Equation. W e define the limiting 47 probability density fo r a n even t given the state v ariable 48 x by the so-c a lled hazard function h x ( x, t ) where t de- 49 notes explicit dep endence on time due to external input 50 following [5, 6]. Here, we assume x has a shot-no ise-lik e 51 dynamics, which is widely used as a mo del of spike in- 52 duced neuronal a daptation [6] 53 ˙ x := − x ( t ) /τ + q P k δ ( t − t k ) , (1) where δ is the Dirac delta function, t k is the time of the 54 k th even t, and q is the q ua n tile change in x a t each even t. 55 The dynamics of x deviates from standard treatments of 56 shot-noise (such as in [7]) in that the ra te of even ts has 57 a dep endence on x a s expressed by the haza r d function 58 h x ( x, t ). It is s traight for ward to show that the dis tr ibu- 59 tion of x in an ensemble, denoted by P r( x, t ), is governed 60 by 61 ∂ t Pr( x, t ) = ∂ x [ x τ Pr( x, t )] + h x ( x − q , t ) Pr( x − q , t ) − h x ( x, t ) Pr( x, t ) . (2) Much insight can be gained b y apply ing the method of 62 characteristics [8] to establish a link betw een the state 63 v a riable x and its time-like v ar iable t x . F or Eq. (1) 64 we define t x = η ( x ) := − τ ln( x/q ), where by d dt t x = 1. 65 When an e v ent o ccurs, t x 7→ ψ ( t x ), wher e ψ ( t x ) = 66 η ( η − 1 ( t x ) + q ) = − τ ln( e − t x /τ + 1) with its inverse 67 given by ψ ( t x ) − 1 = − τ ln( e − t x /τ − 1 ). Thu s, we de- 68 fine h ( t x , t ) := h x ( η − 1 ( t x ) , t ). This tr ansformation of 69 v a riables to t x elucidates the connection of the mo del to 70 renewal theory . Here, the r eset condition after ea c h even t 71 is not t x 7→ 0 (renewal) but t x 7→ η ( x + q ) [5]. Therefore, 72 the v ariable t x that we may call a ’pseudo age’ is a general 73 state v aria ble tha t no longer re presen ts the time since the 74 last even t (ag e). T ra nsforming v ariables in Eq. (2) from 75 x to t x yields in the s teady state 76 ∂ t x Pr( t x ) = − h ( t x ) P r( t x ) + (1 − Θ 0 ( t x ))[ h ( ψ − 1 ( t x ) P r( ψ − 1 ( t x ))] , (3) where Θ 0 ( t x ) is the Heaviside step function, and for con- 77 venience w e defined ψ − 1 ( t x ≥ 0 ) ≡ 0. An efficient al- 78 gorithm for so lving Eq. (3) is given in [6]. W e denote 79 this so lution b y Pr eq ( t x ). F ur ther , the time-like trans- 80 formation in Eq. (3) allo ws computation of the ISI b y 81 analogy to the re ne wal theory [6] and also per mits the 82 compariso n to the master eq ua tion for a renewal pr o - 83 cess a s g iv en in Eq. (6.43) in [9]. The distribution 84 2 of t x just prior to an even t is a quantit y o f interest 85 and it is deriv ed as Pr ∗ ( t x ) = h ( t x )Pr eq ( t x ) /r eq , where 86 r eq = R h ( t x )Pr eq ( t x ) dt x is a normalizing co nstan t and 87 also the pro cess intensit y or r ate o f the ensemb le. Simi- 88 larly , o ne can derive the distribution of t x just a fter the 89 even t, Pr † ( t x ) = Pr ∗ ( ψ − 1 ( t x )) d dt x ψ − 1 ( t x ) [6]. Then the 90 relationship betw een t x and the ordinary ISI distribution 91 can b e written as 92 ρ (∆) = Z + ∞ −∞ h ( t x + ∆)Ω( t x + ∆)Pr † ( t x ) dt x , (4) where Ω( t x + ∆) = e − R ∆ t x h ( t x + u ) du . Now the n th moment 93 µ n of the distribution and its co efficient of v aria tion C v 94 can b e numerically determined. 95 Counting Statistics . In order to der ive the coun t dis- 96 tribution, we gener alize the elegant approa c h for deriving 97 the moment g enerating function as intro duced in [10]: let 98 ρ n ( t n , t n x | t 0 x ) b e the joint probability density g iv en its ini- 99 tial state t 0 x , w her e t n stands for time to n th even t and t n x 100 is the cor respo nding adaptive state of the pr oces s . Ther e - 101 after, o ne ca n recur siv ely derive 102 ˜ ρ n +1 ( s, t n +1 x | t 0 x ) = Z ˜ ρ n ( s, t n x | t 0 x ) ˜ ρ ( s, t n +1 x | t n x ) dt n x , (5 ) where ˜ ρ n +1 ( s, t n +1 x | t 0 x ) = L [ ρ n +1 ( t n +1 , t n +1 x | t 0 x )] and L is 103 the Lapla ce transform with resep ect to time, assuming 104 ˜ ρ 1 ( s, t 1 x | t 0 x ) = ˜ ρ ( s, t 1 x | t 0 x ) [10]. Next, defining the opera- 105 tor P n ( s ) and applying Bra-Kat notation as suggested 106 in [10], le a ds to t he Laplace tra nsform o f n th even ts or- 107 dinary density 108 ˜ ρ n ( s ) = h 1 | P n ( s ) | P r † i = h 1 | [ P ( s )] n | Pr † i , (6) where the op erator P a ssocia ted with ˜ ρ ( s ), which inter- 109 estingly c o rresp onds to the moment g enerating function 110 of the sum of n non-indep enden t interv a ls ˜ f n ( s ) as de- 111 fined in [11]. Now, following E q s. (2.1 5 ) in [11 ] Lapla ce 112 transform o f count distribution denoted a s ˜ P ( n, s ). 113 The F a no facto r provides an index f or t he quan tifica- 114 tion of the co un t v aria bilit y . It is defined as J T = σ 2 T /µ T , 115 where σ 2 T and µ T are the v aria nce and the mean o f 116 the num b er of event s in a certain time windo w T . It 117 follows fro m the additive prop ert y of the exp ectation 118 that µ T = R T 0 r ( u ) du and assuming constant fir ing rate 119 µ T = r eq T . T o calculate the second moment of ˜ P ( n, s ), 120 we require ˜ A s = P k ˜ ρ k ( s ), thus 121 ˜ A s = h 1 | P ( s )( I − P ( s )) − 1 | Pr † i , (7) where I is the identit y op erator . Note, a ssuming a re- 122 newal interv a l dis tribution in Eq. (4) o ne obtains ˜ A r s = 123 ˜ ρ ( s ) / (1 − ˜ ρ ( s )) and L − 1 [ r eq ˜ A s ] = r eq A ( u ) is the joint 124 density of an even t at time t and another even t at time 125 t + u . Thu s, the a uto corr lation of even ts is A ( u ) = 126 r eq [ δ ( u ) + A ( u )]. Now, by using Eq. (7 ) and the Eq. 127 (3.3) in [11], the second momen t of the count statistics 128 can b e derived. Thus, we obtain the F ano factor 129 J T = 1 + (2 /T ) R T 0 ( T − u ) A ( u ) du − r eq T , (8) The asy mptotic prop ert y of F = lim T →∞ J T can be de- 130 rived from the result s tated in Eq . (7.8) in [11] as 131 lim s → 0 [ ˜ A s − 1 / ( µ 1 s )] = C 2 v [1 / 2 + P ∞ k =1 ξ k ] − 1 / 2 , (9) where ξ k is the linear cor relation coefficient betw een tw o 132 k lagg ed in terv als. P r o vided the limit ex its, w e find F = 133 C 2 v [1 + 2 P ∞ k =1 ξ k ] in [12]. 134 Sup erp osition. W e now generalize our results on the 135 counting statistics to the s uperp osition of indep enden t 136 po in t pro cesses. This is of pra ctical interest in all cases 137 where w e observe super imposed even ts that stem from 138 m ultiple independent pr ocess , e.g. in photon detec tio n 139 devices, o r in the case of a p ostsynaptic neuro n that re- 140 ceives conv erg ing inputs f ro m multiple lines. W e study 141 the sup erp osition o f k stationar y , orderly , a nd indep en- 142 dent pro cesses. The ensemble pro cess will hav e a ra te 143 ˇ r = P k i =1 r i and following Eq. (4.18) in [13] ˇ A ( u ) = 144 ˇ r + ˇ r − 1 P k i =1 r i [ A i ( u ) − r i ]. Here, for the sak e of simplic- 145 it y , we derive the desired relationship b etw een C 2 v and 146 the ensemble ˇ F for k identical pro cesses. T o this end, 147 we plug ˇ r and L [ ˇ A ( u )] in to the Eq. (9) a nd therefore 148 it b ecomes lim s → 0 [ ˜ A s − 1 / ( µ 1 s )] = CV 2 [1 / 2 + Ξ] − 1 / 2 , 149 where CV and Ξ = P ∞ i =1 Ξ i are the co efficient of v a ri- 150 ation and the in terv al correlations of the superimp osed 151 pro cess. Note that the left hand side of this equation and 152 Eq. (9) ar e sim ular . Thus, we obtian 153 CV 2 [1 + 2 Ξ] = C 2 v [1 + 2 P ∞ i =1 ξ i ] . (10) The left hand side o f Eq. (10) is indeed the F ano fa ctor 154 ˇ F of the ensemble pro cess as desired. No w, [13] s ug- 155 gests as k → ∞ , CV 2 → 1. Interestingly , if a ll individual 156 pro cesses fullfill the renewal conditio n, it follows from 157 Eq. (10) that ˇ F = C 2 v = [1 + 2Ξ], and therefore if C 2 v 6 = 1 158 the p opulation activity is non-renewal with Ξ < 0 (Ξ > 0) 159 for pro cesses with C 2 v > 1 ( C 2 v < 1). This imp ortant find- 160 ing explains the numerical observ ation in [14] rega rding 161 emergance of non-renewal pro cesses as the result of the 162 sup e rpositio n op eration. 163 A daptation in a Neu r onal En s emble. In [6] it has 164 bee n shown by an adia ba tic elimina tion of fa st v ari- 165 ables that the master eq uation desc r iption of a detailed 166 neuron mo del including voltage dy namics, conductance- 167 based s y napses, a nd spike-induced a daptation reduces to 168 a sto c hastic po in t pro cess sim ular to Eq. (3). The cor re- 169 sp onding hazar d function can be approximated as 170 ˆ h x ( x ) = a t exp( − b t x ) , (11) where a t and b t are determined b y the time dependent 171 statistics of inputs [5] and the equilibrium rate co nsis- 172 tency equa tion r eq ≈ ˆ h x ( r eq q τ ) [6] with the solution 173 r eq = W ( abq τ ) / ( bq τ ) , (12) 3 0.1 1 20 0 1 Expected count in T Fano factor Theory Renewal J T 0.1 1 20 Experiment Shuffled J T FIG. 1. Ad aptation reduces the F ano factor of the ensemble process. Left Magenta : J T for arbitrary observ ation time T according to Eqs. (8) and (11 ) with bq = 1 . 4, a = 5 . 0 an d τ = 400 ms . Blue : F ano factor for equiv alent renew al ensemble process with interv al distribution of Eq. (4). Squar e Dots : Numerically estimated F ano factor for sup erposition of the 5 realization runs of the full-detailed adaptive neuron model as in [6]. Dash-dotte d line : C 2 v . Righ t M agent a : Empirical ˆ J T estimated from the p ooled spike trains of 5 cortical neurons. Blue : F ano factor for th e p ool of shuffled spike trains. Dash- dotte d : A v erage C 2 v of the 5 individual spike trains. where W is the L am b ert function. In the case of v an- 174 ishing adaptation ( bq → 0 ) the pro cess apporaches the 175 Poisson pro cess with r eq → a . 176 W e show in [5] that the adaptation dynamics in Eq . (1) 177 pro duces negative ser ial correlatio ns ξ k < 0. The 178 strength of s erial co rrelation deca ys with increas ing lag 179 k and dep ends on the mean adaptation, E [ x ] = r eq q τ . 180 Such a v anishing of negative seria l interv a l co r relations 181 with increasing lag is well s upported by a large b ody of 182 exp erimen tal evidence [4]. The depa rture from the re- 183 newal prop erty induced by adaptation r educes the F ano 184 factor Eq. (8) fo r the single pro cess as well a s for the 185 po pulation mo del of sup erimpos ed p o cesse s . 186 W e v alidate our theo retical result of the r educed F ano 187 factor in a set o f exp erimental spike trains of N = 5 188 in vivo intracellular r ecorded neurons in the somatosen- 189 sory cortex of the rat. The sp ontaneous activity of each 190 of these neurons shows nega tive ser ial interv al correla- 191 tions [15] where the empir ic a l sum ov er co r relation co- 192 efficients a moun ts to an av era ge P 10 i =1 ξ i = − 0 . 28. W e 193 construct the p opulation activity b y sup erimpo s ing all 194 5 s pik e trains. Therea fter, w e e s timate the F ano factor 195 as a function of the obser v ation time a nd compare it to 196 the case where , prior to sup erp osition, r e new al statistics 197 is enfor ced for each individual neur o n through in terv al 198 shufflin g. Our exp erimen tal o bserv atio n in Fig. 1 (Right) 199 confirms the theore tica l prediction of a reduced F ano fac- 200 tor simular to individual neur ons [16] in the population 201 level. 202 Benefits for Neur al Co ding. W e provide thre e ar gu- 203 men ts that demonstrate how the mechanism of s pik e- 204 frequency ada ptation b enefits neura l pr o cessing and po p- 205 ulation coding. First, our result of a reduced F a no factor 206 ˇ F < C 2 v for t he po pulation a ctivit y of stationa ry adap- 207 tive pro cesses ( bq > 0) directly implies a reduction o f the 208 noise in the neurona l p opulation rate co de. Our ana lysis 209 of a set of co rtical data sug gests a reductio n of > 50% 210 for long o bserv a tio n times. The reduction of J T in Fig. 1 211 bec omes significant ev en for s mall observ ation times of 212 ≈ 2 a verage interv als, whic h is a r elev ant time scale for 213 the transmission of a p opulation r a te signal. This res ult 214 is reminiscent o f an effect t hat has previously been a c- 215 knowledged as noise sha ping and weak stimuli detection 216 expressed in the re ductio n of the low frequency pow er in 217 a spe ctral a nalysis of spik e trains with negative ser ial in- 218 terv al corr e la tions [17]. Our result confirms their findings 219 at the p opulation level. 220 Our seco nd a rgument is concerned w ith the trans mis- 221 sion of a p opulation rate signal. W e may define a func- 222 tional neural ensemble by the common po stsynaptic tar- 223 get neur on that receives the co n verging input of all en- 224 semble mem b ers. T o elucidate the pos ts y naptic effect 225 of adaptation we simplify the ensemble autocor relation 226 function A ( u ) following [18] with an exponential a ppro x- 227 imation 228 ˆ A ( u ) = r eq δ ( u ) + [( F − 1) / 2 τ c ] e x p( − u/τ c ) , (13) where the second term is the approximation of r eq A ( u ). 229 F or given observ ation time window u , and τ c the r educ- 230 tion of F implies that ˆ A r u < ˆ A u . Therefore , the p ostsy- 231 naptic neuron r e ceiv es inputs fro m a n adaptive ensemble 232 that expr esses an extended autoco rrelation structure as 233 compared to the inputs from a no n-adaptive ensemble. 234 F ollowing the theory on the effect of input autocor rela- 235 tion o n s ig nal transmission in spiking neurons as devel- 236 op ed in [1 8], a longer τ c reduces the input current fluc- 237 tuations and this facilitates a faster and more reliable 238 transmission of the modulated input rate sig na l by the 239 po stsynaptic targ et neuron. 240 Finally we a r gue that a p ostsynaptic neuron ca n b et- 241 ter deco de a small change in its input if the presynaptic 242 neurons are ada ptiv e. T o this end, w e compute the in- 243 formation gain o f the po stsynaptic activity , b et ween t wo 244 counting distributions o f an ada ptiv e pr esynaptic ensem- 245 ble when ˆ h x ( x ) is adia batically trans ferred to ˆ h x ( x − ǫ ) 246 with a small change ǫ in the input ensem ble. It is co n- 247 venien t to use ˜ ρ n ( s ) w hich as socia ted with counting dis- 248 tribution ˜ P ( n, s ). Thus, we apply the Kullback-Leibler 249 divergence to the Eq. (6) b efore and after the adiabatic 250 change in the input 251 D K L ( ˜ ρ ǫ n || ˜ ρ n ) = P i ˜ ρ ǫ i ( s ) ln( ˜ ρ ǫ i ( s ) / ˜ ρ i ( s )) . (14) Using E q . (7) w e obtain D K L ( ˜ ρ ǫ n || ˜ ρ n ) = A ǫ s [ln( A ǫ s / A s )]. 252 Due to Eqs. (1 ) and (12) , the mean adaptatio n after the 253 change is E [ x ǫ ] = τ qr ǫ eq . If ǫ > 0 it follows that r ǫ eq ≥ r eq . 254 Therefore the mean adaptation level incre a ses and the 255 adapted pro cess exhibits str onger negative serial corre- 256 lations and A ǫ s > A s . Thus, by Eq. (13), it is straigh t 257 forward to deduce that D K L > D r K L , for rene wal and 258 adaptive pro cesses with identical interv al distributions. 259 4 0 25 0 30 r eq (Hz) r eq ε (Hz) Poisson adapted r eq (Hz) D KL / κ 0 25 0 1.25 FIG. 2. Information gain p er spike due to adaptation. Left : T ransfer of equilibrium rate f or fixed ǫ change of the input in adaptive and Poi sson mod el. Righ t : Kullback-Leibler Di- verge nce p er extra spike as th e measure of information gain for n th even t densit y of adapt iv e and P oisson processes while u = 200ms and ǫ = 0 . 01nS with the same initial r eq and κ = r ǫ eq − r eq . W e now compute the informatio n gain of the adap- 260 tive ensem ble pr ocess relative to a matched Poisson rate 261 mo del. F or differen t in itial rate v alues r eq we as s ume a 262 small but fixed increase ǫ in the input that w e express in 263 parameter c hanges a ǫ and b ǫ in Eq. (11) a s outlined in 264 [5]. This leads to a n increa se κ = r ǫ eq − r eq in rate that 265 is effectively consta nt ov er a wide range of initial v a lue s 266 r eq (Fig. 2, Left). In the r ate mo del, assuming the same 267 initial v alue of r eq , the sa me input step leads to a higher 268 equilibrium rate increase κ P oisson > κ , whic h dep ends 269 on the inital rate (Fig. 2 , Left) b ecause the r ate mo del 270 lacks a mechanism of self-inhibition, which in the adap- 271 tive mo del counteracts the r ate incr e ase. T he r eafter, we 272 compute the Kullback-Leibler divergence for b oth mo d- 273 els and no rmalize it by the change in the output r ate κ . 274 The result in Fig. 2 (Right) shows that D K L /κ is larger 275 for the adaptive mo del than for the rate mo del across 276 the ra nge o f tested input ra tes. Thus, the informa tion 277 per extra spike is larger in the adaptiv e ensem ble than in 278 the rene wal ensemble, and a p ostsynaptic neur on can dis - 279 criminate small c hanges ǫ more efficiently , even though 280 the absolute change in fir ing rate is low er . 281 Discussion. Our r e sults p oint out a new asp ect of spike 282 frequency adaptation that b enefits the reliable transmis- 283 sion and p ostsynaptic deco ding o f the neural p opulation 284 co de. This a spect adds to the k no wn prop erties of co m- 285 pression and temp oral filtering of sensory input signa ls 286 [19] in spike frequency adapting neurons. The s pecific 287 result o f Eq. (10) is a ls o of pra ctical consequence for the 288 empirical analy sis o f the so-ca lled multi-unit activity . By 289 estimating F ano factor and serial co rrelations we readily 290 obtain an estimate of the average C v and serial correla - 291 tion of the individua l pro cesses. 292 W e developed a new formalism to tre a t even t emit- 293 ting systems that a re influenced by a s lo w state v aria ble, 294 and we provide a num b er of useful general results on 295 the hig her o rder even t statistics of sup erimp osed renewal 296 and non- renew al ev ent pro cesses, which are applicable 297 to a wide r ange of event -ba sed systems in nature [5]. 298 The der iv ation o f the state dependent haz ard a nd mas- 299 ter equation [6] assumes incoherent input fluctuation as 300 in the mea n-field theory , where common input is neg li- 301 gible. T r eating a net work with coheren t fluctuations as 302 encountered in finite size netw ork s requir es an alterna tiv e 303 deriv ation of the hazar d function [5]. 304 A cknow le dgment s . W e thank Carl v an V reeswijk and 305 Stefano Car danobile for v aluable c o mmen ts, and Clemens 306 Boucsein for making the da ta av a ilable for re-analy s is. 307 This work was funded by the German Ministry of E du- 308 cation and Rese a rc h to BCCN Berlin. 309 ∗ Corresponding author: farzad@zedat.fu-b erlin.de 310 [1] See Refs. [6-7] and [19-26] in F. C. S oler, F. J. Ro driguez, 311 and G. Zumofen, Phys. Rev. A 78 , 053 813 (2008). 312 [2] V. N. Livina, S. Ha vlin, and A. Bun de, 313 Phys. Rev. Lett. 95 , 20 8501 (2005). 314 [3] B. Lindner and T. S c hw alger, Phys. Rev. Lett. 98 , 210603 315 (2007). 316 [4] Y. H. Liu and X. J. W ang, J Comput Neurosci 10 , 2 5 317 (2001); F. F arkho oi, M. F. Strub e-Bloss, and M. P . 318 Naw rot, Phys. Rev. E 79 , 021905 (2009). 319 [5] Supplementa ry material at . 320 [6] E. Muller, L. Buesing, J. Schemmel, and K. Meier, 321 Neural Comp. 19 , 2958 (2007). 322 [7] C. W. Gardiner, Handb o ok of Sto chastic Metho ds , 3rd ed. 323 (Springer, 2004). 324 [8] M. de Kamps, Neural Comp. 15 , 2129 (2003). 325 [9] W. Gerstner and W. M. Kistler, Spiking Neur on Mo dels , 326 1st ed. (Cam bridge Universit y Press, 2002). 327 [10] C . va n V reeswijk, in Analysis of Par al lel Spike T r ains , 328 Springer Series in Comput. Neu rosci. (Sp rin ger, 2010). 329 [11] J . A. McF adden, J. of the Roy al Stat. So c. Series B 24 , 330 364 (1962). 331 [12] D. R. Cox and P . A . W. Lewis, The Statistic al Analysis 332 of Series of Events , Methuen’s Monographs on A p plied 333 Probabilit y and Statistics (Methuen, London, 1966). 334 [13] D. Cox and V. Isham, Point Pr o c esses , 1st ed. (Chapman 335 and Hall/CRC , 1980). 336 [14] A. J. La wrence, J. of the Ro yal Stat. Soc. Se- 337 ries B 35 , 306 (1973); H. Cˆ a teau and A. D. 338 Reyes, Phys. Rev. Lett. 96 , 058101 (2006); B. Lind ner, 339 Phys. Rev. E 73 , 022901 (2006). 340 [15] M . P . Naw rot, C. Boucsein, V. Rodriguez- 341 Molina, A. Aertsen, S. Gruen, and S. Rotter, 342 Neuro computing 70 , 1717 (2007). 343 [16] R. Ratnam and M. E. Nelson, 344 J. Neurosci. 20 , 6672 (2000). 345 [17] M . J. Chacron, B. Lindn er, and A . Longtin, Phys. 346 Rev. Lett. 9 2 , 080601 (200 4); I. F uwa p e and 347 A. B. Neiman, Ph ys. R ev . E 78 , 051922 ( 2008); M. J. 348 Chacron, A . Longtin, and L. Maler, J. Neurosci. 21(14) , 349 5328 (2001). 350 [18] R. Moreno-Bote, A. Renart, and N . Parga , 351 Neural Comput 20 , 1651 (2008); R. Moreno, 352 J. de la R ocha, A. R enart, and N . Parga , 353 Phys. Rev. Lett. 89 , 28 8101 (2002). 354 [19] B . N. Lundstrom, M. H. Higgs, W. J. S pain, and A . L. 355 F airhall, Nat Neurosci 11 , 1335 (2008). 356