Fast Inference of Interactions in Assemblies of Stochastic Integrate-and-Fire Neurons from Spike Recordings

F ast Inference of In teractions in Assem blies of Sto c hastic In tegrate-and-Fire Neurons from Spik e Recor dings R. Monasson 1 , 3 , S. Co cco 2 , 3 Abstract W e present t wo Ba yesian p rocedures to infer the interactio ns and ex t ernal currents in an assembly of sto c hastic integrate-and-ﬁre neurons from the recording of their spiking activit y . The ﬁrst procedu re is based on the exact calculatio n of the most lik ely time courses of the neuron membrane p oten tials cond itioned by the recorded spikes, and is exact for a v anishing noise v ariance and for an instan taneous sy nap- tic integration. The second pro cedure t ak es into account the presence of ﬂuctuations around the most lik ely t ime courses of the p otentials , and can deal with mod erate noise levels . The running time of b oth pro cedures is prop ortional to the number S of spikes multiplied by th e squared num b er N of neurons. The algorithms are va lidated on synthetic data generated by netw orks with k no wn couplings and currents. W e also reanalyze previously pu blished recordings of the activity of the salamander retina (in- cluding fr om 32 to 40 neurons, and from 65 , 000 t o 170 , 000 spikes ). W e stu d y the dep endence of the inferred interactions on the membrane leaking time; th e d iﬀerences and similarities with th e classical cross-correlation analysis are discussed. 1 Intr o duction Over th e p ast decades, m ulti-electro de recordings ( T aketani and Baudry , 2006) hav e unv eiled th e nature of th e activity of p opulations of neural cells in v arious sy stems, such as the vertebrate retina (Schnitzer and Meister, 2003), cortical cultures (T ang et al, 200 8), or the pref rontal co rtex (P eyrache et al, 2009). The observ ation of substantial correlations in th e ﬁring activities of neurons has raised fun damen tal issues on their functional role (Romo, Hernand ez, Zai nos and Salinas, 2003 ; Averbeck, Latham and P ouget, 2006). F rom a stru ctural p oin t of view, a challenging problem is to infer the netw ork and the strengths of the functional interac tions b etw een the neural cells from 1 Laboratoire de Ph ysique Th´ eorique de l’ ENS, CNRS & UPMC, 24 rue Lhomond, 75005 Pa ris, F rance 2 Laboratoire de Physique Statistique de l’ENS, CNRS & UPMC, 24 rue Lhomond, 75005 Pa ris, F rance 3 The Simons Cen ter for Systems Biology , Institute for Adv anced Study , Einstein Dri ve, Prince- ton NJ 08540, USA 2 the spiking activit y (Fig. 1A) . Po w erful inference pro cedures are needed, capable to handle massive data sets, with millions of spikes emitted by tens or hundreds of neurons. A classical approach to infer funct ional neural connectivity is t hrough the study of pairwise cross-correlations ( Perkel , Gerstein and Moore, 1967; Aersten and Gers tein, 1985). The approach was applied in a v ariet y of neural systems, including the audi- tory midbrain of the grassfrog ( Ep ping and Eggermon t, 1987), the salamander retina (Briv anlou, W arland and Meister, 1998), th e primate and rat prefron tal cortex (Con- stanti dinidis, F rano wicz and Goldman-Raking, 2001; F ujisa wa, Amarasingham, Har- rison and Buzsaki, 2008). Other app roac hes, capable of taking into account netw ork- mediated eﬀects, were proposed based on concepts issued from statistics and graph the- ory (Seth and Edelman, 2007; Dahlhaus, Eichler and Sandk ¨ uhler, 1997; Sameshima and Baccal´ a, 1999; Jung, N am and Lee, 2010), information theory (Bettencourt, Stephens, Ham and Gross, 200 7), or statistica l physics (S c hneidman, Be rry , Segev and Bialek, 2006; Sh lens et al, 2006). An alternative approach is to as sume a particular dyn amica l mo del for the spike generation. The generalized linear mo del, which represen ts th e generation of spikes as a Po isson process with a time-d ependent rate is a p opular framework (Bro wn, N guy en, F rank, Wilson and Solo, 2001; T ruccolo et al, 2005; Pillo w et al, 2008). The Integrate- and-Fire (IF) model, where spik es are emitted according to the d ynamics of th e mem- brane p otentia l is anoth er n atural candidate (Gerstner and Kistler, 2002; Joliv et, Lewis and Gerstner, 2004). The problem of estimating the mod el parameters ( external cur- rent, va riance of the noise, capacitance and conductance of the membrane) of a single stochastic IF neuron from the observ ation of a spike train has received a lot of at- tention (Paninski, Pillo w and Simoncelli, 2004; Pill ow et al. , 2005; Mullowney and Iyengar, 2008; Lansky and Ditlevsen, 2008). F ew stu dies hav e focused on the inference of interactions in an assembly of IF neurons (Mak aro v, P anetsos and de F eo, 2005). Recently , we prop osed a Ba y esian algorithm to infer the interac tions in a netw ork of stochastic p erfect integrators when the synaptic integra tion is instantaneous and t he noise is v anishingly small (Co cco, Leibler an d Monasson, 2009). In t he present work we introduce a Bay esian algorithm to infer the coup lings and the external currents in an assembly of leaky IF neurons, and in presence of moderate input noise (Fig. 1A). The computational time grows as the pro duct of the num b er of recorded spikes, and the square of the number of n eu rons. W e val idate th e algorithm on synthetic d ata, and apply it to real recordings of the ganglion cell activit y in the salamander retina, presented with natural v isual stimuli, an d in the absence of stimulus (sp on taneous activity). 2 Materials and Metho ds 2.1 Deﬁn ition of the Leaky Integrate-and-Fire model In the Leaky Integrate-and-Fire (LIF) model, the membrane p otential V i ( t ) of neuron i at t ime t ob ey s the ﬁrst-order diﬀerential equation, C dV i dt ( t ) = − g V i ( t ) + I syn i ( t ) + I i + η i ( t ) (1) where C and g are, resp ectively , th e capacitance and conductance of the mem brane. The ratio τ = C /g is t h e membrane leaking t ime. I syn i ( t ) is th e synaptic current coming 3 time i j j i J ij J ji I I i j ? neuron A a b c neuron 2 neuron 1 B Fig. 1 A. Extra-cellular r eco rdings give access, through spike-sorting, to the times t i,k of the spike s emitted by a p opulation of neu rons. W e w an t to infer the v alues of the in teractions J ij and external inputs I i of the netw ork most l i k ely to hav e generat ed the recorded activit y . B. Example of ﬁring activit y of N = 2 neurons. The top panel sho ws three spike s emitted by neuron 1. Pane ls a, b, c show p ossible activities of neuron 2, with equal av erage ﬁring rates but wi th diﬀerent timings. from the other neurons and entering the neuron i at time t : I syn i ( t ) = X j ( 6 = i ) J ij X k δ ( t − t j,k ) (2) where J ij is th e strength of the connection from neuron j onto n euron i (Figure 1A); t j,k is the time at whic h n euron j ﬁres its k th spike. W e assume that synaptic inputs are instantaneo usly integrated, i.e. that th e synaptic integration time is muc h smaller than all the other time scales, including τ . O ur metho d for inferring t h e interactions relies on this assumpt ion, and should b e mod iﬁed in the presence of synaptic integration kernels with temp oral ﬁ ltering. I i is a constant external curren t ﬂowing into neuron i (Fig. 1A), and η i ( t ) is a ﬂuctu ating current, mo deled as a Gaussian noise p ro cess: h η i ( t ) i = 0, h η i ( t ) η j ( t ′ ) i = σ 2 δ ij δ ( t − t ′ ). The n oise standard dev iatio n, σ , has here the dimension of a current times the square ro ot of a time. A n alternative d eﬁnition w ould consist in rescaling σ with a time dep endent factor, e.g. √ τ ; our d eﬁnition allo ws us to reach the p erfect in tegrator limit ( τ → ∞ ) while keeping σ ﬁx ed . The neuron i remains silent as long as V i remains b elo w a threshold p oten tial V th . If the threshold is crossed at some time t 0 , i.e. V i ( t 0 ) = V th , then a spike is emitted , and th e p oten tial is reset to its resting v alue: V ( t + 0 ) = 0. The d ynamics then resumes follo wing ( 1). 2.2 Likelihoo d of th e spiking times for given interactions and currents Let J = { J ij } and I = { I i } denote the sets of, respectively , t he intera ctions and currents. Let t i,k ∈ [0; T ] b e the time at which n eu ron i emits its k th spike; T is the duration of the recording. H o w can w e infer the in teractions and currents from the observ ation of the spiking activity? Consider the raster plots in Fig. 1 B. In pattern a , the timings of the spik es of neuron 1 do n ot seem t o be correlated to t he activity of neuron 2. H ence, we ma y guess that there is no interactio n from neuron 2 to neuron 1 ( J 12 = 0). In pattern b , a spike of n euron 1 is lik ely to follo w a spike of neuron 2, whic h suggests that the interaction J 12 is p ositiv e. Conv ersely , in p attern c , it seems that the ﬁring of neuron 2 hinders the ﬁring of neuron 1, which indicates that J 12 has a negative v alue. 4 This crude reaso ning can be made mathematically rigorous in the framew ork of statistical inference (Cov er and Thomas, 200 6). Let us deﬁne the likelihood P ( T |J , I ) of a set of spiking times, T = { t i,k } , giv en J and I . According to Ba yes rule the m ost lik ely couplings and currents, ˆ J and ˆ I , given the set of spiking t imes T can b e inferred through the maximization of P ( T |J , I ) 1 . Due to the statistical in dependen ce of the noises η i from neuron to neuron, the likelihoo d P of T giv en J , I can b e written as the pro duct of First-Passage Time (FPT) probabilities, P ( T |J , I ) = Y i,k p F P T ( t i,k +1 | t i,k , { t j,ℓ } , { J ij } , I i ) . (3) Here p F P T denotes th e probability that V i crosses V th for the ﬁ rst time at time t i,k +1 , starting from 0 at time t i,k and conditioned to the in puts from the other neurons at times t j,ℓ , with t i,k < t j,ℓ < t i,k +1 . As the synaptic integ ration is in stantaneous, an incoming spike from neuron j results in a (p ositiv e or negative) jump of the p otenti al V i by J ij /C ; p F P T can th erefore be in terpreted as the FPT density probability for a one- dimensional Ornstein-Uh len b ec k pro cess with a time-dep endent force. It is imp ortan t to stress th at the presence of the p roducts ov er th e spike in terv als in (3) d oes n ot entail that the spiking times are indep endent. Consider now the p otential V i ( t ) during the inter-spike interv al (ISI) [ t i,k ; t i,k +1 ]. The b oundary conditions are V i ( t + i,k ) = 0 (reset of t he p otential right after a spike), and V i ( t − i,k +1 ) = V th (condition for ﬁring). At intermediate times, the p otentia l can take an y v alue smaller than V th . The logarithm of the p robabilit y of a dynamical path (time course) of t he p otentia l ov er the k th ISI of n eu ron i is, after m ultiplication by the v ariance σ 2 of the noise, L [ V i ( t ); k , T , J , I ] = − 1 2 Z t i,k +1 t i,k dt η i ( t ) 2 (4) = − 1 2 Z t i,k +1 t i,k dt  C dV i dt ( t ) + g V i ( t ) − I syn i ( t ) − I i  2 , according to the Gaussian nature of th e n oise η i ( t ) and to the dynamical equ ation of the LIF (1). 2.3 Dyn amica l eq uations for the optimal p oten tial and noise While n o exact expression is k no wn for p F P T , it can be analytically approximated by the contribution of t h e most probable d ynamical path for th e potential, V ∗ i ( t ) (Pa ninski, 2006). This approximation b ecomes exact when th e standard deviation σ of the n oise is small. The id ea is to replace the distribution of p ath s for the p oten tial V i ( t ) with a single, most likely p ath V ∗ i ( t ), which we call optimal. W e n o w ex plain h ow to derive V ∗ i ( t ) through the condition that the log–probabilit y L (4 ) is maximal. Let us assume ﬁrst that V ∗ i ( t ) < V th . Then, the deriv ativ e of L in (4) with resp ect to V ∗ i ( t ) must v anish, whic h gives δ L δ V i ( t )     V ∗ i = − C 2 d 2 V ∗ i dt 2 ( t ) + g 2 V ∗ i ( t ) + C dI syn i dt ( t ) − g I syn i ( t ) − g I i = 0 . (5) 1 W e consider here that the a priori measure o ve r the couplings and currents is ﬂat. 5 W e no w turn this second order diﬀeren tial equation for the optimal p otential into a ﬁrst order diﬀerential eq uation at the p rice of introducing a new function, η ∗ i ( t ), and a n ew ﬁrst order diﬀerential equation for this fun ction. It is straigh tforw ard to chec k that the solution of C dV ∗ i dt ( t ) = − g V ∗ i ( t ) + I syn i ( t ) + I i + η ∗ i ( t ) (6) is a solution of th e optimization equ ation (5) if η ∗ i ( t ) fulﬁlls dη ∗ i dt ( t ) = g C η ∗ i ( t ) = η ∗ i ( t ) τ , (7) where τ is t he membrane leaking t ime. The similarit y b et ween eqns (1) and (6) allo ws us to interpret η ∗ i ( t ) as a current n ois e. How ever, t his noise is n o longer sto c hastic, but rather it follo ws the deterministic path solution of (7). W e will, t h erefo re, in the follo wing refer to η ∗ i ( t ) as the optimal noise. η ∗ i ( t ) corresponds to the most lik ely v alue the n oise takes given the set of spiking times. S olving (7) sho ws that the optimal noise is an ex ponential function of t h e time: η ∗ i ( t ) = η exp(+ t/τ ) ( V ∗ i ( t ) < V th ) , (8) where η is a constant, which w e call noise co eﬃcien t. It may h appen that the optimal p oten tial only reaches the threshold without ac- tually crossing it at intermediate times. When this is th e case, the optimal p oten tial equals V ∗ i ( t ) = V th and its deriv ativ e with respect to the time v anishes. The expression for the op t imal noise can b e then read from (6 ) , and is given by η ∗ i ( t ) = g V th − I syn i ( t ) − I i ( V ∗ i ( t ) = V th ) . (9) Equation (9) ensu res that the potential d o es not cross the t h reshold va lue at a time t < t i,k +1 . Despite their apparen t simplicit y , eqns (6,8,9) are not easy to solve, due mainly to t he interpla y b etw een t h e tw o regimes, V ∗ < V th and V ∗ = V th , mentioned ab o ve. The d etermination of V ∗ i ( t ) was achieved numeric ally by P aninski for a single neuron (Pa ninski, 200 6). W e n o w sketch th e p ro cedure to determine V ∗ i ( t ) rapidly , even for tens of neurons. The pro cedure relies on th e searc h for con tacts, that is, times at which the optimal p oten tial touches the threshold. There are tw o t yp es of con tacts: contacts coinciding with a synaptic inpu t (the p oten tial touches the threshold at time t j,k ), and con tacts arising in b et w een tw o inputs. In the absence of leak age, only the former type of contacts matter, and a search procedure to lo cate those isolated-time contacts w as prop osed by Co cco, Leibler and Monasson (2009). In the presence of leak age, b oth types of contacts hav e to b e taken into accoun t. The searc h pro ced ure is more complex, and is explained b elo w. 2.4 Fixed Threshold pro cedure: optimal p aths for the p otentia l and the n oise W e assume in this Section that the couplings and curren ts are known. Consider neuron i at time t ∈ [ t i,k ; t i,k +1 ], where k is the index of t he ISI. The initial and ﬁnal con- ditions for the optimal p oten tial are: V ∗ i ( t + i,k ) = 0 and V ∗ i ( t − i,k +1 ) = V th . In b etw een, V ∗ i ( t ) obeys the LIF evolution eq uation (6 ) with an optimal ‘noise’ η ∗ i ( t ). η ∗ i ( t ) can 6 0 1 V * / V th 0 1 t i,k t i,k+1 0 1 1 1.6 6 0 .25 -.25 η * /(gVth) 0 2 4 .25 -.25 0 0 V th V 0 0,5 p s A A P B P A C D V th M A A 0 1 2 3 4 5 time (sec) 0 0.2 0.4 0.6 0.8 1 V/V th E Fig. 2 A & B. Sket che s of the optimal potentials V ∗ (top, blac k curves) and noises η ∗ (bottom, black curves) for one neuron receiving sev eral inputs from tw o other neurons ( red and green im pul s es, middle, with J r ed = − J gr een = . 2 C V th ). T he membrane conductance is g = . 8 I /V th . The jump i n th e optimal noise consecutiv e to an activ e con tact is alwa ys posi tive. A. Illustration of Passive (P) and Active (A) con tacts. B. Comparison with numerical simulations, a v eraged o ver ∼ 5 , 000 samples, f or ¯ σ = . 07 (red), . 18 (purple), . 36 (blue); noise curv es are av eraged ov er th e time-windo w ∆t = . 015 τ . C. Dashed l ines r epresen t p ossible paths f or the p oten tial when the noise standard deviation, σ , does not v anish. T he amplitude of the ﬂuctuations of the p oten tial around V ∗ at the mid-p oint of the ISI i s symbolized by the double arrow l ine. D. Probability p s ( δt | V ) that an Ornstein-Uhlenbeck pro cess starting from V does not cross the threshold V th for a time δ t = 3 . 5 τ . Paramet ers are g V th /I = 1 . 2, ¯ σ = . 15. The tange nt line to p s in V = V th crosses the p s = 1 2 line in V M th . E. System of t w o IF neurons, with g V th /I 1 = 1 . 5 , g V th /I 2 = 2 ., J 12 / ( C V th ) = . 1 , J 21 = 0 , ¯ σ = . 25. The dashed and full blac k curves represent the optimal poten tials for neuron 1 calculated by , r esp ectiv ely , the Fixed and Mo ving Thr eshold pro cedures; for the l att er, V M th is shown in red. One random realization of the membrane p oten tial (av eraged o ver a 10 msec time–window) is shown for comparison (blue curve) . b e interpreted as a non-sto c hastic, external, t ime-depend ent cu rren t to b e fed into th e neuron in order to drive its potential from 0 to V th , given th e sy naptic couplings. The expressions for the optimal ‘noise’ are give n by (8) when V ∗ i ( t ) < V th , and ( 9) when the optimal p oten tial V ∗ i ( t ) is equ al to the threshold v alue. When V ∗ i ( t ) reac hes the threshold at a time coinciding wi th an incoming spike, the co eﬃcien t η in (8) ma y abruptly change through an active c ontact ; the notion of active contact is illustrated in the simple case of a neuron receiving a single spik e in App endix A .1. The p otenti al V ∗ i ( t ) may also touch the threshold without crossing it, and the noise ma y remain constant ov er some time in terv al; we call suc h an event p assive c ontact . That th e p oten tial can brush, or remain at the threshold level without prod ucing a spik e is made possible by th e σ → 0 limit. W e wil l discuss later on the v alidit y of th is calculation, and how to mo dify it when t he noise standard deviation, σ , do es not v anish. Both typ es of con tacts are shown in Fig. 2A . Let us explain how the positions of active and passive contacts can b e determined. Let t 1 < t 2 < . . . < t M b e the emission times of the spikes arriving from th e neurons intera cting with i du rin g the t ime interv al [ t 0 ≡ t i,k ; t M +1 ≡ t i,k +1 ], and J 1 , J 2 , . . . , J M the corresp onding synaptic strengths 2 . Let V 0 = 0 b e the initial v alue of the potential, and m 0 = 1 b e the index of the ﬁrst input spike. If t h e time is small enough the optimal p oten tial is surely b elo w th e threshold va lue. A ccording to (8) the optimal noise is an 2 Due to the limited temp oral resolution of the measurement tw o inputs of amplitudes J and J ′ can apparen tly arriv e at the same tim e; if so, we consider, based on mo del (1) and (2 ), that a single input of amplitude J + J ′ en ters the neuron. 7 exp onen tial with noise co eﬃcien t η , and the optimal p otentia l is obtained by solving (6) with the result, V i ( η , t ) = V 0 e − ( t − t 0 ) /τ + M X m = m 0 J m C e − ( t − t m ) /τ θ ( t − t m ) + I i g (1 − e − ( t − t 0 ) /τ ) + η g sinh  t − t 0 τ  (10) where θ is the Heaviside function. It is tempting to look for the va lue of η such t hat a spike is emitted at time t M +1 , deﬁned by the implicit equ ation V i ( η , t M +1 ) = V th . How ev er, the corresp onding potential might not b e b elo w threshold at all intermediate times t 0 < t < t M +1 . Instead, we look for the smallest noise capable of driving the p oten tial from its initial value V 0 into contact with th e threshold: η ∗ = min  η : max t 0 0), and neuron 2 is indep endent from the activity of n euron 1 ( J 21 = 0). In presence of a strong noise, the optimal p oten tial calculated from the Fixed Threshold proced ure quickly reaches a stationary val ue close to V th , while the random p otential obtained from sim ulations ﬂ uctuates around a muc h low er level (Fig. 2E). The presence of a strong noise biases the membrane p otential to low er v alues to prevent early sp iking. A heuristic approac h to repro duce this bias consists in decreasing th e threshold from V th to a time- and con text-dep endent v alue, V M th . W e now exp lai n ho w this moving threshold, V M th , is determined . Consider ﬁrst a neuron with no synaptic inpu t, fed with an external current I , during the inter-spike interv al [ t i,k ; t i,k +1 ]. W e call p s ( δ t | V ) the probability that the p oten tial, taking v alue V at time t i,k +1 − δ t , remains b elo w the threshold at any larger time t , with t i,k +1 − δ t < t < t i,k +1 . This probabilit y dep ends on the current I , and can b e expressed for an arbitrary level of noise, σ , as a series of parab olic cylinder functions (Alili, Patie and P edersen, 2005). Figure 2D show p s as a function of V for some chara cteristic v alues of th e parameters. The probabilit y of surviva l, p s , sharply decreases to zero when V gets close to the th reshold, V = V th . W e mo del this outcome 10 by th e fol low ing app ro ximation, which in volv es a new, eﬀective threshold V M th : w e consider that the processes starting from a v alue of the p otentia l V > V M th will not survive for a time delay δ t . In other w ords, the true threshold, V th , is changed into a ’mo ving’ threshold, whic h is a function of the current I , the time δ t , and t h e parameters g , C, σ . A simple w a y to d eﬁne V M th is to lo ok at t he intersection of the t an gent line to p s in V = V th with, say , the p s = 1 2 line 3 ; the resulting expression for V M th is given in App endix E. Figure 2E shows the output of th e Moving Threshold pro cedure on the simple 2-neuron system d escribed above. The optimal p oten tial, ’pushed’ do wn b y the mo ving threshold V M th is much lo wer th an in the Fixed Threshold approach and in muc h b etter agreemen t with the random realization of t h e mem brane p otential . More details are give n in Section 3.1.4. T o account fo r the existence of synapt ic inputs, we ma y choose the p arame ter I en tering th e calculation of p s and V M th to b e the v alue of th e eﬀective cu rren t I e i = I i + X j ( 6 = i ) J ij f j , rather than the external current I i itself. Here, f j is the a v- erage ﬁ ring rate, d eﬁned as th e num b er of spikes ﬁred by neuron j d ivided b y the duration T . Contrary to the external curren t I i , th e eﬀectiv e current I e i takes into account th e (a vera ge) inp ut cu rren t coming from other n eurons. This c hoice was done in the numerical exp erimen ts rep orted in the Results section. T o further sp eed up the calculations, w e deriv e the value of V M th for d iscrete-va lues dela ys δ t only; in a discrete interv al, V M th is kept to a constant v alue. Alternative h euristic app roac hes to deal with the p resence of mo derate noise can b e prop osed. In App endix E w e introduce a cost-funct ion for the eﬀective curren t, whose eﬀect is also to decrease th e optimal potential. These approaches are eﬀective when the optimal p otentia l calculated by the Fixed Threshold pro cedure quickly saturates to a level close to V th . More p recisel y , w e ex pect the Mo ving Threshold pro cedure to b e eﬃcien t if th e m embrane leaking time is smaller or comparable to the ISI, and the leaking current, ≃ g V th , is larger or equ al to th e external current, I . 2.7 Maximization of th e log-lik elihoo d to infer the interactions and currents The Fixed or Moving Threshold procedures allo w us to calculate the optimal paths for the p oten tial and th e noise, given the coup lings and currents. Knowl edge of those paths gives us also access to th e log arithm of the likelihood P in the σ → 0 limit, L ∗ ( T |J , I ) = lim σ → 0 σ 2 log P ( T |J , I ) = X i,k L [ V ∗ i ( t ); k , T , J , I ] = − 1 2 X i,k Z t i,k +1 t i,k dt η ∗ i ( t ) 2 (16) Since L ∗ in (16 ) inv olves t he sum o ver diﬀerent neurons, the maximization o ver the couplings J ij and the current I i of neuron i can be done independently of the other couplings J i ′ j and currents I i ′ ( i ′ 6 = i ). F orma lly , we are left with N indep endent inferences of the most like ly couplings and current for a single neu ron, in presence of the spikes emitted b y the N − 1 oth er neurons. A s a consequence neurons ‘decouple’ in 3 This c hoice i s arbitrar y; other v alues, ranging from 1 4 to 1 hav e b een tried, do not quali- tativ ely aﬀect the results presen ted later in this article. 11 the invers e problem: the couplings J ij and th e current I i of neuron i can b e inferred indep endently of the other couplings J i ′ j and currents I i ′ ( i ′ 6 = i ). L deﬁn ed in (4) is a negative-semideﬁnite qu adratic function of its arguments V i ( t ) , J ij , I i . It is thus a concave function of the couplings an d the currents. This prop- erty h olds for L ∗ (16) (Bo yd and V and en b erghe, 2004). I n order to infer the most likely current I i and couplings J ij , we start from an arbitrary initial v alue e.g. I i = J ij = 0. The full path of the optimal noise, η ∗ i ( t ), ov er all the inter-spike interv als k of neuron i , is calculated follo wing the ab o ve pro cedure. W e then up date th e coup lings and th e current using the Newton-Raphson metho d to maximize log P , i.e. to minimize th e integ ral of the squared optimal n ois e, see (16). Con vergence follo ws from the conca vity prop ert y stated ab o ve. The pro cedure requires the expressions for the gradient and the Hessian matrix of log P with respect to the couplings J ij and the current I i , whic h can b e calculated exactly from (16) and (10 ). Note that log P is p iecew ise continuousl y tw ice-diﬀerentiable; while the gradient is con tinuous for al l J ij and I i , th e Hessian matrix is boun d ed and n ega tive, and ma y discontin uously jump due to a change of the conta ct p oin ts. Knowl edge of the H essi an matrix is also imp ortan t t o determine how reliable are th e v alues of the inferred parameters. 2.8 Accuracy on the inferred parameters When the va riance of the noise, σ 2 , v anishes the inferred p arameters cann ot dev iate from their most likely v alues. Ho wev er, for small but non zero σ , deviations are p ossi- ble 4 . The probability for such d eviations can b e estimated from the expansion of L ∗ around its maximum. W e introduce for eac h neuron i , th e N -dimensional vector v i whose comp onents are: v ( i ) i = I i τ , and v ( i ) j = J ij for j 6 = i . The multiplication of the current by the membrane leaking t ime ensures that all comp onen ts can b e ex pressed in units of a coupling. Similarly we call ˆ v ( i ) the vector obtained when the current and couplings take their most likel y va lues, th at is, maximize L ∗ . Let us call H ( i ) j,j ′ = − 1 σ 2 ∂ 2 L ∗ ∂ v ( i ) j ∂ v ( i ) j ′ ( T | ˆ J , ˆ I ) . (17) the Hessian matrix of L ∗ . The parameters v ( i ) j are normally distributed around their most likely v alues, with a cov ariance matrix given b y h  v ( i ) j − ˆ v ( i ) j  v ( i ) j ′ − ˆ v ( i ) j ′  i =  H ( i )  − 1 j,j ′ . (18) In particular, the error bars on th e inferred parameters are given b y the diagonal elemen ts of the invers e of H ( i ) . Note that, if t he val ue of σ is not known, formulas (17) and (18) can still b e u sed to compare the error bars b etw een eac h other. As the en tries of H ( i ) scale linearly with the duration T of the recording, or, more precisely , the num b er S of recorded spikes the uncertaint y on the inferred parameters will d ecrea se as S − 1 / 2 . A d etaile d sp ectral analysis of σ 2 H ( i ) /S in t h e case of weak 4 Note that the inferred parameters might b e less sensitiv e than the time course of the potential to the noise lev el σ . The reason is that the corrections to the log-like liho od L ∗ , to the low est order in the noise v ariance σ 2 , do not dep end on the current and inte ractions (Appendix D). 12 couplings, reported in App endix C, sho ws that the largest eigen v alue, λ max , is related to the ﬂ u ctuations of the eﬀective current, I e i = I i + X j ( 6 = i ) J ij f i,τ j , (19) where f i,τ j = 1 T X k,ℓ : t i,k 0). 15 0 0.1 0.2 0.3 0.4 0.5 0.6 fraction p of connections 0 0.05 0.1 0.15 ε s (J)/J 0 gV th /I = .8 gV th /I = .5 gV th /I = .1 A -0.2 0 0.2 -0.2 0 0.2 inferred J ij -0.2 0 0.2 true couplings J ij gV th /I=.2 gV th /I=.8 r=.03 r=.03 C-1 (.89) (.41) FT FT C-2 5 10 15 20 t/ τ 0 1 V * / V th D 0 0.25 0.5 0.75 1 1.25 gV th / I 0 0.05 0.1 0.15 0.2 0.25 ε s (J) / J 0 r =.15 r =.03 r =.15 B -0.2 0 0.2 -0.4 -0.2 0 0.2 inferred J ij -0.2 0 0.2 true couplings J ij r=0.15 gV th /I=1.2 r=0.15 gV th /I=1.2 FT MT (.10) C-3 C-4 5 10 15 20 t/ τ 0 0.5 1 V, V * , V th M / V th E Fig. 4 Results from the Fixed (empty squares) and Mov ing (full squares) Threshold algo- rithms on a r andom net work of N = 40 coupled neurons; the maximal amplitude of synapses is J 0 = . 2 C V th . The error on the couplings, ǫ s ( J ) /J 0 , is plotted as a function of the f raction p of connections ( A ) and conductance ov er curr en t r atio, g V th /I ( B ). Each simulated data set con tains S = 5 10 5 spike s, which is larger than the cross-ov er size S c.o. ; the symbol size correspond to the ﬂuctuations estimated from ten diﬀerent data sets f or the same net work of in- teractions. C. Inferred inte ractions vs. true v alues of J ij for v arious v alues of g V th /I and r , and for one random net w ork w i th a f raction p = . 2 of connections. Dashed li nes hav e slop e uni ty . Pa nels C- 1, C-2, C-3 sho w the results of the Fixed Threshold (FT) pro cedure; the slop es of the best linear ﬁts (full lines) are i ndicat ed b et w een parent hesis. Panel C-4 sho ws the outc ome of the Moving Thr eshold (MT) pr ocedure; ev en if m ultiplied by 10, the FT couplings of Panel C- 3 are in muc h worse agreemen t wi th the true interactions than the MT couplings. D. Optimal potentials V ∗ obtained wi th the Fixed Threshold pro cedure f or g = . 1 I /V th (dashed curve) and g = 1 . 2 I /V th (full curve) , and f or one arbitrari ly chosen neuron among the N = 40 neural cells; the noise r atio is r = . 15. E. Comparison of a random realization of the poten tial V (red) with the optimal p ote ntial V ∗ (blac k) obta ined with the M o ving Threshold V M th (green) procedure. The net work of int eractions, the s pi king times, and the ar bitrarily c hosen neuron are the same as the ones in D for g = 1 . 2 I /V th . The time-av erage of V M th is ≃ . 93 V th . connection j → i is preserved. A t the end of the construction process, the av erage num b er of outgoing (or incoming) n eigh b ors of a neuron is p ( N − 1). Each existing connection is then assigned a syn aptic weigh t, uniformly at random o ver the interv al [ − J 0 ; J 0 ]. All n eu rons receiv e the same external current I . I n addition, the membrane conductance, g , is now diﬀerent from zero. The v alues of p, J 0 , I , g , and σ are chosen so that the netw ork remains b elow saturation. W e h ave also p erformed sim ulations where the interaction graph is drawn as ab o ve, but each neuron i is chosen to b e eith er excitatory or inhibitory with eq ual probabilities. The outgoing interactions from i hav e al l the same sign, and random amplitudes in [0; J 0 ]. The p erformance of our inference algorithms are qualitatively similar for b oth mod els. Figure 4A sho ws the error on the couplings inferred with the Fixed Threshold algorithm, ǫ s ( J ), as a function of the fraction p of connections, for three va lues of the membrane conductance ove r curren t ratio. The error roughly increases as √ p , that is, the number of conn ections in the netw ork. This scaling suggests that muc h of the 16 inference error is du e to non-zero coup lings. This ﬁ nding agrees with Fig. 3C, which sho wed th at th e inferred interactions b et ween uncoupled neurons was very small in the g = 0 case. T o better understand the performance of the alg orithm, we compare in Fig. 4C the inferred interactions J ij with their true val ues for the 1560 oriented pairs j → i of a randomly drawn netw ork of N = 40 neurons, with p = . 2 an d J 0 = . 2 C V th . When the ratio g V th /I is small compared to un it y , the quality of the inference is very goo d (Fig. 4C-1). F or larger ratios g V th /I the inferred coup lings are still strongly correlated with their t rue v alues, but are ap p ro ximately rescaled by an ov erall factor < 1, corresponding to th e av erage slop e of th e linear regression in Fig. 4C-2. As g V th /I increases, this factor decreases and the inference error grows (Fig. 4A ) . Figure 4B show s that the inference error on the interactions increases not only with g V th /I but also with the noise ratio r . F or large v alues of r , th e n etw ork can sustain activit y even when g V th > I , and the inference error can take large v alues (up per cu rv e in Fig. 4B). I n th is regime, the couplings found by the Fixed Threshold algorithm b e- come small, and the inferred cu rren t I i gets close to g V th . The corresp onding p otential V ∗ i ( t ) rises sharply , in a time τ , to a v alue slightly b elow threshold, I i /g , with small ﬂuctuations due to the synaptic inputs. This phenomenon can be seen in Fig. 4D, whic h compares the optimal p otential of a neu ron for tw o diﬀerent v alues of mem brane conductance. As d iscussed in the Metho ds section, this b eha vior is a consequence of the σ → 0 limit taken in the calculation of the optimal p otentia l; when σ , or r , is not small, the p otential is unlikely to stay close to the threshold for a long time without prod ucing a spike, see Fig. 2B. In t he n ext paragraph, w e analyze the results of t h e Mo ving Threshold inference pro cedure. As a conclusion, zero couplings are p erfectl y inferred, while the amplitude of large (p ositiv e or negative) interactions can b e und erestima ted by the Fixed Threshold al - gorithm, esp ecially so when the noise is strong. How ev er, th e relative ordering of the intera ctions is essentiall y preserved by th e inference pro cedure. 3.1.4 Infer enc e err or with the Moving Thr eshold pr o c e dur e The Mo ving Threshold pro cedure was t ested in Fig. 2E on an asymmetric system of tw o IF neurons ( J 12 / ( C V th ) = . 1 , J 21 = 0) in the presence of a strong noise, see d escri ption in caption and Section 2.6. While the Fix ed Threshold pro cedure erroneously inferred that both interactions v anish, the Mo ving Threshold correctly inferred the sign and the order of the magnitude of the coup ling: J inf err ed 12 / ( C V th ) = . 2 ± . 1. The inferred currents were within 10% of their t rue v alues. These results were obtained from a large num b er S of spikes to a void ﬁnite- S eﬀects. The syn thetic data used in Fig. 4B were generated with tw o diﬀerent v alues of the noise ratio, r . W e estimate the relative ﬂuct uations of the p oten tial around the optimal path, av eraged o ver all the inter-spike interv als in th e d ata set, using form ula ( 14 ), and ﬁnd p h ( V i − V ∗ i ) 2 i V th ≃  . 028 for r = . 03 . 138 for r = . 15 (25) for all v alues of g V th /I comprised betw een .9 and 1.25. Hen ce, t he relative ﬂuctu atio ns cannot b e neglected when r = . 15. Figure 4B show s the inference error obtained from the Moving Threshold algorithm as a function of the membrane condu ctance for that v alue of the noise ratio . Not su rp risingl y , the Mo ving Threshold pro cedure is more accurate than the Fixed Threshold algorithm. 17 In the Mo ving Threshold algorithm, the optimal p otentia l is constrained to remain b elo w a certain threshold, V M th , which dep ends on the time preceding th e next spike and on t he eﬀective cu rrent I e i . Figure 4 E shows the v al ues of the mo ving threshold V M th and of the optimal p otential V ∗ i for a few spike interv als of the same neuron as in Fig. 4D. As exp ected, the v alue of V ∗ i ( t ) lie s substantia lly further aw a y from the threshold V th than in the Fixed Threshold proced ure. In addition, Fig. 4E show s a random realiz ation of the potential V i ( t ), ob t ained throu gh numerical integration of the LIF diﬀerential equ atio n (1), for the same neuron i . Although V i is stochastic, the comparison of severa l inter-spike interv als ind icates that V ∗ i ( t ) and V i ( t ) are in fair statistical agreement. T o inv estigate in more details the origin of the inference error on the couplings for large v alues of r and g V th /I , we p lot in Fig. 4 C the inferred val ues of the interaction J ij vs. t he true v alue for every pairs j → i of a randomly draw n netw ork of N = 40 neurons. The in teractions inferred b y the Fixed Threshold algorithm are about ten times smaller than their t ru e v alues (Fig. 4C-3). The use of t h e Mo ving Threshold proced ure leads to a sp ectacular impro vemen t for p osi tive-v alued coup lings (Fig. 4C-4). While p ositiv e couplings are accurately inferred, th e magnitude of negative couplings is often ov erestimated. These negative couplings are resp onsible for most of the error ǫ s in Fig. 4B. F rom th e Ba yes ian p oin t of view, when τ is small er than the av erage ISI , negative -va lued couplings are indeed intrinsicall y harder to infer th an p ositiv e-v alued ones. A p ositiv e inp u t driv es the p oten tial closer to th e threshold, whic h strongly reduces the ISI . Convers ely , a negative inp ut drives the p otenti al down, and a spike is unlikely to occur b efore the p oten tial ﬁrst relaxes t o its av erage level I /g after a time of th e order of τ . H ence, the inﬂuence of a negative inp ut is hardly seen in t h e increase of the IS I when τ is smaller than the a verage IS I . W e present an analytical calculation supp orting this argument in Section 3.2.2. 3.2 App lica tions t o multi-electrode recording data W e now apply our algori thm to m ulti-electro de recordings of th e ganglion cell activity of the salamander retina. Two data sets were considered. The ﬁrst one, hereafter re- ferred to as Dark (data courtesy of M. Meister), rep orts the spontaneous activity of 32 neurons for 2,000 seconds, and consists of 65 , 525 spikes (Schnitzer and Meister, 2003). In the second exp erimen t, referred to as Natural Movie (data courtesy of M. Berry), a retina w as presen ted a 26.5 second-long mo vie, repeated 120 times, and the activity of 40 neurons w as regis tered for the whole duration of 3,180 seconds (Schneidman, Berry , Segev and Bialek, 2006). Natu ral Mo vie includ es 172 , 521 spikes. The ﬁring rates, av- eraged o v er the p opulation of recorded n eurons, have similar v alues in the tw o data sets: f ≃ 1 . 02 spikes/sec in Dark, f ≃ 1 . 35 spikes/sec in Natural Movie. These t wo d ata sets we re analyzed in a prev ious work (Cocco, Leibler and Monas- son, 2009) with th e p erfect integrator mo del ( g = 0) and the Fixed Threshold algo - rithm. In this section w e extend the analysis to the case of the LIF mod el and use b oth the Fixed and Mo ving Threshold approac hes. I n particular w e sho w that the LIF mod el is capable of inferri ng th e asymmetry of th e interactions, which is see n in t h e cross-correlog rams but wa s n ot obtained with the p erfect integrator mo del. Moreov er w e discuss error bars on th e inferred couplings and the fact that strong negative in - teractions are more diﬃcult to infer than positive-v alued couplings. W e stress that the couplings we infer a priori d epend on the stim ulus. Co cco, L eibler and Monasson (2009) 18 hav e studied how the in teractions in ferred with the p erfect in tegrator model dep ended on th e stimulus based on the analysis of tw o recordings on th e same retina, namely the sp on taneous activity and random ﬂickeri ng squares. An alternative approac h to disentangl e stim ulus-indu ced and structural contributions to the couplings w ould b e to consider a time- and stimulus-depend en t external current I ( t ) (Section 4.4). The v alue of the mem brane leaking t ime τ strongly aﬀects the num b er of contacts and the run n ing time of th e alg orithm. It takes about 40 seconds to infer the cu rren ts and the interactions from either Dark or N atural Mo vie when τ ≃ 1 sec with one core of a 2.8 GH z I ntel Core 2 Quad deskt op compu ter, and ab out 10 times longer w hen τ = 100 msec. The number of p assive contacts of the opt imal p otential computed by the Fixed Threshold p ro cedure quickly decreases as τ increases. It is divided by ≃ 20 when th e membrane leaking time increases from 100 msec to 10 sec for b oth data sets. In comparison, the number of activ e contacts is less sensitive to the v alue of τ . W e ﬁ nd that th e ratio of th e number of conta cts p er neuron and p er second ov er th e avera ge ﬁring rate takes similar v alues for b oth d ata sets. F or τ = 1 msec, this ratio is ≃ 2 . 00 for Dark, and ≃ 2 . 04 for Natural Movie. The number of passive conta cts is smaller with the Mo ving Threshold algorithm, while the num b er of active contacts remains rath er unchanged compared to its v alue with the Fixed Threshold p rocedure. On the ov erall, the running time of the Moving Threshold proced ure is higher du e to t h e calculation of the time-dep endent threshold V M th . Knowl edge of the v ariance of the noise is requ ired for the Mo ving Threshold algo- rithm. The v alue of σ could, in principle, b e determined from exp erimental measures of the ﬂ uctuations of the synaptic curren t, b ut is u nkno wn for the tw o recorded data sets a v ailable to us. W e c ho ose σ so that the relativ e ﬂuctuations of th e p otential around the optimal path V ∗ i are less than 10%. W e compute these ﬂuctuations by av eraging (14) ov er all IS I and all neurons i in the p opulation. The corresponding v alue of the dimensionless stand ard deviation of the noise (13) are: for Dark, ¯ σ = . 13 , . 12 , . 11 for, respectively , τ = 200 , 100 , 20 msec; for Natural Movie, ¯ σ = . 15 , . 14 , . 12 for, resp ectiv ely , τ = 200 , 100 , 20 msec. 3.2.1 Amplitudes of the inferr e d inter act ions and curr ents Figure 5 A sho ws the av erage v alue of th e current and of the interactio n strength as a function of th e membrane leaking t ime. As exp ected with the Fixed Threshold inference proced ure, we ﬁnd that the av erage v alue of the couplings decreases as τ gets small. This eﬀect v aries from neuron to n euron: the closer I i is to g V th , the smaller are the couplings J ij . T o compare the matrices of couplings J, J ′ inferred with th e Fixed Threshold algorithm for diﬀerent val ues of τ , we use th e correlation co eﬃcien t (Hu bert and Baker, 1979) R ( J, J ′ ) = co v( J, J ′ ) p co v( J, J ) cov( J ′ , J ′ ) , (26) where co v( J, J ′ ) = N ( N − 1) X i 6 = j J ij J ′ ij −  X i 6 = j J ij  X i 6 = j J ′ ij  . (27) Identical matrices corresp ond to R = 1, and uncorrelated matrices giv e R = 0. R is indep endent of the scale of th e coupling matrices J and J ′ , i.e. R ( aJ, a ′ J ′ ) = R ( J, J ′ ) for any a, a ′ > 0; th eref ore, R is sensitive to the relative amplitudes of the couplings J ′ and J and not to t h eir absolute d iﬀerences. W e choose J to b e th e coupling matrix 19 0.01 0.1 1 10 τ (sec) 0 0.2 0.4 0.6 0.8 1 I/(gV th ) 0.01 0.1 1 10 τ (sec) 0 0.05 0.1 0.15 0.2 0.25 J/(CV th ) A FT MT MT FT 0.01 0.1 1 10 100 τ (sec) 0.4 0.6 0.8 1 R Dark Natural Movie B 0 1 2 3 4 5 latency(i,j)/ τ 0.1 1 10 100 -J ij τ = 200 msec τ = 100 msec τ = 20 msec D -0.6 -0.3 0 0.3 0.6 J ij (Moving Threshold) -0.6 -0.3 0 0.3 0.6 J ij (Fixed Threshold) -0.6 -0.3 0 0.3 0.6 C-1 C-2 0 0.05 0.1 0.15 0.2 0.25 latency (sec) 0 30 60 number of pairs 0 0.05 0.1 0.15 0.2 0.25 0 30 60 Dark Natural Movie E Fig. 5 Amplitudes of the i n teractions and currents in Dark (full circles) and Natural M ovie (empt y di amonds) . A. Average v alue of the current (left) and ro ot mean square v alue of the coupling (right) as a function of the membrane leaking time τ . P oints corresp onding to the Fixed Threshold (FT) pro ced ure are joined by full li nes, while dashed lines indicate the results from the Moving Thr eshold (MT) algori thm. Note that the currents are larger in Dark than i n Natural Movie. B. Correlation co eﬃcien t R (26) b et wee n the couplings at l eaking tim e τ and with no leak age. C . Comparison betw een the interactions J ij found with the Mo ving (x-axis) and the Fixed (y-axis) Threshold pr ocedures, for Dark (C-1) and Natural Mov ie (C-2). The dashed li ne is the x = y line, and τ = 20 msec. D. Strongly negat ive couplings J ij vs. latency o ver the m em brane leaking time τ for three v alues of τ . Couplings were obtained usi ng the Moving Threshold pro cedure, and corresp ond to the N atu ral Movie data set. Only interact ions J ij < − . 1 ar e considered; there are, r espectiv ely , 16, 28, and 60 such couplings f or τ = 200, 100, and 20 msec. The v alue of the slop e of the b est l inear ﬁt l og( − J ij ) = α latency( i, j ) /τ + β , sho wn b y the dashed l i ne, is α = 0 . 95. E. Distributions of the latencies (28) b et we en neurons in Dark (top) and Natural M o vie (bottom). Only latencies larger than 5 ms ec are tak en into accoun t in the histograms. in the absence of leak age and J ′ to b e the coupling matrix for a give n τ . The va lue of R as a function of τ is shown in Fig. 5B. Even for τ = 20 msec, the coupling matrix is substantial ly similar to th e one obtained with the p erfect integrator mo del ( R = . 6 for Dark, R = . 5 for Natu ral Movie). Despite the ov erall change in the amplitude of the inferred couplings, the relative ordering of the coup lings with t he pair indices ( i, j ) is largely indepen den t of τ , especially so for Dark. How ev er, for sp eciﬁc pairs of neurons, the intera ctions m ay strongly dep end on τ . Such a d ep endence eﬀect will b e illustrated in Section 3.2.3, and can b e related to the temp oral structure of the corresponding cross-correlog rams. The a verage v alue of the in teractions cal culated by the Moving Threshold algo- rithm does not decrease when τ gets smaller, and is larger than the one obtained from the Fixed Threshold p ro cedure (Fig. 5A). T o b etter un d erstand this discrepancy , we compare in Fig. 5C th e in teractions infe rred with b oth algorithms for every pairs of neurons in t h e Dark and Natural Movie data sets when τ = 20 msec. The agreemen t b e- tw een b oth pro cedures is very goo d for p ositiv e and strong couplings. Couplings which are slightl y p ositiv e with the Fixed Threshold pro cedure generally h a ve a larger v alue 20 with t he Moving Threshold pro cedure. This oﬀset is resp onsible for the diﬀerences in the a verage v alues of the interactions found in Fig. 5 A. I n add ition, in N atural Movie, negative -va lued couplings often have a stronger amplitude with the Mo ving Threshold proced ure. W e ﬁn d, in both approac hes, a few negative and v ery strong couplings. The amplitude of th ose extreme coup lings increases very quic kly as the membrane leaking time decreases. The emergence of strong negative intera ctions with the lo wering of τ can b e related to the presence of long latencies betw een the emission of spikes. W e deﬁne the latency of neuron i with resp ect to n euron j as the smallest d ela y b etw een a spike emitted by j and a later spike ﬁred by i , latency( i, j ) = min k,ℓ : t i,k / τ B -0.4 -0.2 0 0.2 0.4 J/(CV th ) 0 10 20 30 40 50 60 d/dJ x CV th / τ C Fig. 6 A. Marginal log-likelihoo ds L c ( I 1 ) (top panel), and L c ( J 1 , 27 ), L c ( J 1 , 20 ), L c ( J 1 , 4 ) (from l ef t to right in the b ott om panel) f or Natural Movie, and τ = 200 msec. Dashed lines correspond to the best ﬁts wi th a single quadratic function. The most li kely v alue for the current is ˆ I 1 = 1 . 14 g V th . The most li k ely v alues for the interactions are: ˆ J 1 , 27 = − . 11, ˆ J 1 , 20 = . 01, and ˆ J 1 , 4 = . 22, in units of C V th . The oﬀset on the ve rtical axis has been ch osen so that all maxima are at heigh t L c = 0. B. Average v alue of the ﬁrst-passage tim e t F P T after a s ynap tic entrance of amplitude J . C. Deriv ative of t F P T with resp ect to J . The parameters of the neuron are: g V th /I = 1 . 2, r = . 15, τ = 85 msec (full l i ne) and 20 ms ec (dashed line). The deriv ativ e i s maximal ar ound J opt / ( C V th ) = 1 − I / ( gV th ) ≃ . 167. The same pro cedure can obviously be used to obtain the error bar on the current I i . W e now illustrate this approach on the N atural Mo vie data set, and one arbitrarily chos en neuron, i = 1. Three interactions, representative of, resp ectiv ely , p ositiv e, weak, and negative couplings, were singled out among the 39 couplings incoming onto neu- ron 1. Figure 6 A shows the marginal log-likelihoo ds L c ( J 1 , 4 ), L c ( J 1 , 20 ), and L c ( J 1 , 27 ), in addition to L c ( I 1 ). F or all four parameters, the marginal likelihoo ds can be approxi- mated with parab olas in the vicinit y of t heir maxima. Estimating th e second d eriva tives from those b est qu adratic ﬁts and u sing (30 ), w e obtain ∆I 1 g V th ≃ . 020 ¯ σ , ∆J 1 , 27 C V th ≃ . 023 ¯ σ , ∆J 1 , 20 C V th ≃ . 021 ¯ σ , ∆J 1 , 4 C V th ≃ . 022 ¯ σ . (31) where ¯ σ is th e dimensionless noise level deﬁ n ed in (13) . Hence, the error bars on th e couplings and curren ts hav e very similar v alues. This common v alue dep ends on the noise level, ¯ σ . As discussed in the nex t Section 3.2.1, ¯ σ is exp ected to b e close to, or smaller than unity when τ = 200 msec. Consequen tly , the v alue for J 1 , 20 is compatible with zero, while th e interactio ns J 1 , 27 and J 1 , 4 are non zero, with 99.9999% conﬁdence. A closer insp ection of Fig. 6 A shows that the q ualit y of the quadratic ﬁ t of L c is excellen t for J 1 , 27 and J 1 , 4 , but less so for I 1 and J 1 , 20 . F or the la tter p arameters , it seems that the curv ature of L c takes t w o diﬀerent v alues, dep ending on whether the maximum is app roac hed from the left of from the right. This phen omenon results from the piece-wise structure of the L ∗ function, see Methods section. A practical consequence is that the errors I 1 − I ∗ 1 and J 1 , 20 − ˆ J 1 , 20 are not ev enly distributed around zero; for instance J 1 , 20 is more likel y to b e larger than ˆ J 1 , 20 than it is to b e smaller. Note that strong, negativ e interactions may b e harder t o infer than p ositiv e-v alued couplings, a ph enomenon already underlined by Aersten and Gerstein (1985). The underlying intuition is that th e d uration of the IS I is less aﬀected by an inhibitory input 22 -200 -100 0 100 200 delay t (msec) 0 5 10 15 20 25 H 5,17 -100 -50 0 50 100 delay t (msec) 0 5 10 15 20 H 1,22 A 0.01 0.1 1 10 100 membrane leaking time τ (sec) 0 0.5 1 1.5 symmetry ratio ρ pair 1,22 pair 5,17 B Fig. 7 A. Cross- correlograms H ( t ) f or pairs (5 , 17) and (1 , 22) i n Dark. The cross- corr elograms are normalized such that H ( t ) → 1 for large delays | t | . B. Ratios J ij /J j i of the interactions betw een the neurons 5,17 (top) and 1 , 22 (bottom) as a function of τ . than b y an excitatory input when the mem brane leaking time, τ , is small compared to the av erage val ue of the IS I . W e now p resen t an analytical argument su p porting this intuition. Consider a neuron, fe d w ith an external current I and with nois e va riance equal to σ 2 . Assume a synaptic inpu t of amplitude J is receiv ed at time t = 0. W e call t F P T the av erage v alue of the time at which the neuron will emit a spik e; the calculation of t F P T can b e done using a series of parab olic cylinder functions (Alili, P atie and Pedersen 2005). Figures 6B&C shows th at the dep endence of t F P T on J is muc h weak er for negati ve-v alued J than for p ositive couplings. A s th e set of spikin g times is th e only information we h a ve at our d isposal, the diﬃculty in inferring negative couplings is intrinsic to the Bay esian approac h, and cannot b e circum ven ted by any particular algorithm. 3.2.3 Symmetry of the inter actions and cr oss-c orr elo gr ams The dep endence of th e symmetry of coup lings up on the membrane leaking time τ can b e understo od, to some ext ent, from the structure of the cross-correlograms, that is, the histograms H ij ( t ) of the dela ys t = t i,k − t j,ℓ b et w een the times of the spike s ﬁred by th e t w o neurons i, j in each pair. T o do so, w e consider tw o pairs of neurons in Dark, called pairs (5 , 17) and (1 , 22). Figure 7A sho ws the cross-correlograms H 5 , 17 and H 1 , 22 . P air (5 , 17) is c haracterized by a positive p eak in H , centered in t = 0, and of width ≃ 20 msec. Pair (1 , 22) exhibits a p osi tive p eak of correlations, of the same width, but centered around t ≃ 20 msec. W e plot in Fig. 7B the symmetry rati os of the in teractions in the pairs , ρ 5 , 17 = J 5 , 17 /J 17 , 5 and ρ 1 , 22 = J 1 , 22 /J 22 , 1 . W e ﬁnd that ρ 5 , 17 is, to a large ex ten t, indep endent of τ . Conve rsely , ρ 1 , 22 sharply decreases with decreasing τ and is close to zero when τ = 20 msec, whic h coincides with the typical dela y in the cross-correlogram H 1 , 22 sho wn in Fig. 7A. W e conclude t hat the inference pro cedure is capable of capturing the directionality of the intera ction b etw een t h e neu rons 1 and 22, if τ is small en ough. This results shed some light on the corresp ondence b et ween the intera ctions inferred within the LIF mo del and within the Ising mo del (Schneidman, Berry , S egev and Bialek, 2006; Shlens et al, 2006). Couplings inferred with the p erfect integra tor mo del for Dark are in go od agreemen t with the Ising interactio ns, when th e time is binned into windo ws of width ∆t = 20 msec (Cocco, Leibler and Monasson, 2009). By construction, the I sing model pro duces symmetric interacti ons from the pair-wise correlations of the activities, av eraged of the binning windo w. In the absence of leak age, the Integrate-and- Fire inference algorithm hardly distinguishes b et we en a p ost-synaptic and pre-synaptic 23 ﬁring pattern, and prod uces rath er symmetric coup lings. But as τ decreases, th e LIF couplings may b ecome strongly asymmetric (Fig. 7 B). In this case, the correspond ence b et w een the Ising and LIF couplings b reaks down. The same phenomenon was observ ed in Natural Movie, where d ela ys in th e cross-corre lograms are even stronger. 4 Discussi on In this article, we ha ve p resen ted a pro cedure to in fer th e interactio ns and currents in a netw ork of Leaky Integrate-and-Fire neurons from their spiking activity . The v alidit y of the pro ced ure w as established through numerical tests on syn thetic data generated from netw orks with k no wn couplings. W e hav e also applied our al gorithm to real recordings of the activity of tens of ganglion n eurons in the salamander retina. Though our algorithm is limited to mo derate noise levels and instan taneous sy naptic integ ration, it is fast and can, to our kno wledge, hand le much bigger data sets than the existing inference methods for the sto chas tic IF mo del. It is our in tention to make this algorithm av ailable to the n eurobiolo gy communit y in a near futu re. 4.1 Comparison with previous studies Cross-correla tion analysis ( P erk el, Gerstein and Mo ore, 196 7; Aersten and Gerstein, 1985) consists in study ing the distribution of delays b etw een the spikes of neurons in a p air. This approach has b een u sed to characterize the conn ections b et we en neurons (amplitude, time-scale, d ep endence on distance), or their dynamical ev olution (F u ji- sa wa , Amarasingham, Harrison and Buzsaki, 2008). The analysis d o not require an y com binatorial pro cessing of t he activity of a large part of the neu ral assembly . A s a re- sult, the approac h is not limited to small n et w orks. H ow ever, cross-correlatio n analysis ma y ﬁn d diﬃcult to separate direct correlations from indirect correlations modu lated through interactio ns wi th neurons in the surrounding netw ork (Osto jic, Brunel and Hakim, 2009 ; Cocco, Leibler and Monass on, 2009), or due to common inp u ts (Con- stanti dinidis, F rano wicz and Goldman-Rakic, 2001; T rong an d Rieke, 2008). In statistical approac hes a widely-used concept is the one of causality (Seth and Edelman, 2007). A causal interaction exists from neuron i to n euron j if the kn o wledge of the activit y of i h elps predict the ﬁring of j b ey ond what can be ac hiev ed f rom the activity of j alone. In practice, causal relationships are detected through linear multiv ariate statistical regressions (Sameshima and Baccal´ a, 1999), and may o verlook non-linear d ep endencies. Causal analysis h a ve also diﬃculties in ev aluating t h e strength of the intera ctions. Maxim um en tropy models, whic h deduce in teractions fro m pairwis e correlations only , hav e b een sho wn to accurately repro duce higher-order correlations b et wee n neu - rons in t h e vertebrate retina (Schneidman, Berry , Segev and Bialek, 2006; S hlens et al, 2006; Cocco, Leibler an d Mo nasson, 2009). These mo dels, how ev er, suﬀer from some limitations. Intera ctions are constrained to b e symmetric, and temporal correlations are partially discarded (Marre, El Boustani, F r´ egnac and Destexhe, 2009 ). I n add ition obtaining the interactions from the correlations may b e computationally ve ry hard for large netw orks, though eﬃcient approximate algori thms hav e recently b een developed (Cocco and Monasson, 2010). 24 Generalized linear models (GLM), whic h represents the generation of spikes as a P oisson pro cess wi th a time-depend ent rate, hav e b een applied to v ari ous neural systems (Brow n, Nguyen, F rank, Wilson and Solo, 2001; T ruccolo et al, 2005; Pillo w et al, 2008). The inference of parameters in the GLM framew ork is apparently easier to solv e than for IF models, whic h has made the GLM framew ork very attractiv e. Whether GLM are b etter th an I F mo dels to accoun t for real neural activity , regardless of the computational complexity of b oth inference framew ork, is an imp ortan t issue (Gertsner and Naud, 2009). W e hop e that our wo rk, which makes p ossible to apply the IF mo del to large data sets, will help to answer this qu estion. Approaches to infer model parameters in the IF framework hav e b een so far capable of pro cessing a very limited number of neurons or of spikes. Pillo w et al. ( 2005) inferred the mo del parameters of one sto c hastic IF neuron based on a 50 second-long recording with a pro cedure t olerating an y level of noise; Mak arov, Panetsos and de F eo (20 05) inferred th e connections betw een 5 deterministic IF neurons from a 60 second-long synthetic spike train. I n comparison we hav e analyzed a 3180 -second long recording of the activity of 40 neurons. The ru n ning time of ou r pro cedure increases as N 2 S ∼ N 3 T f , where T is the du- ration of the recording and f is the a vera ge ﬁring rate. R ecen tly , Koy ama and P aninski (2009) hav e prop osed a numerical procedu re for calculating the op t imal potential and inferring t h e interactio ns. In their approac h, the time is discretized in to many time- bins of sma ll duration ∆ , and the v alues of th e optimal p oten tials at those discrete times can b e found by means of the interi or-p oin t metho d for discrete constrained opti- mization problems. The runn ing time of the proced ure, O ( N 3 T / ∆ ), is approxima tely 1 / ( f ∆ ) times larger th an ours. In p ractice, f is of the order of 1 t o 10 Hz, while th e discretization time, ∆ , is of the order of 1 msec; h en ce, 1 / ( f ∆ ) ranges from 100 to 1000. How ev er, th is order of magnitude do es not take into account the existence of multi- plicativ e constants; a comparativ e test of the tw o app roac hes on the same synthetic or real d ata w ould be useful to accurately estimate their runn ing times. F urthermore, the algorithm introduced by Koy ama and Paninski can easily incorp orate the presence of temp oral ﬁltering in the intera ctions. Our pro cedure is, in its present form, v ali d when the integrati on kernel is instantaneous only; considering other syn aptic kernels would require ad ho c mod iﬁ catio ns to the expressions of the opt imal noise and p otential and to the search pro cedure for contacts. 4.2 How to include a ﬁnite integration time One of th e ma jor assumptions in our ap p roac h is that the syn aptic integratio n time, τ s , is v anishingly small. In practice, τ s does not v anish, but migh t often be smaller than the mem brane leaking time, τ , and the av erage ISI. Assume that neu ron i , whose p oten tial V i is close to the threshold V th , receives a spike at time t from another neuron, j , through a strongly excitatory connection J ij > V th − V i . Then , n euron i will reach the threshold level after ha ving received a c harge ∆q = C ( V th − V i ), smaller than J ij . As a consequence, large p ositive in teractions can b e underestimated when the latency of neuron i from neu ron j (28 ) is smaller th an τ s . T o comp ensate for this eﬀect we could introduce a time-dep endent va lue fo r the intera ction, J ij ( t, t i,k +1 ) = J ij min  t i,k +1 − t τ s , 1  , (32) 25 where t i,k +1 ( > t ) is the closest ﬁ ring time of n euron i . Hence the eﬀective interaction J ij ( t, t i,k +1 ) is equal to its nominal v alue J ij only if th e synaptic current has enough time to enter the neuron j , and is a fraction of J ij otherwise. The mo diﬁed pro ced ure will b e correct as long as τ s < τ . If the syn aptic and membrane time-scales are com- parable, one needs to take into accoun t th e complete shap e of the syn aptic integration kernel, K ( t ). Cho osing simple enough integratio n synaptic kernel, such as the piece- wise linear function K ( t ) = 0 if t < 0 or t > τ s , K ( t ) = 2 min( t, τ s − t ) /τ 2 s if 0 ≤ t ≤ τ s , could lead to tractable d ynamical equ ations for the optimal p otential and n ois e. The resolution of t h ose equations is left for future wo rk. 4.3 T ow ards a more realistic in fere nce mod el The inference procedu re that w e ha ve introdu ced here can be extend ed to include realistic features such as a refractory p eriod, τ R . T o do so, we restrict the sum in (10) to t h e sp ikes m entering the neuron i at times larger t h an t 0 + τ R . W e hav e run the mod iﬁ ed inference pro cedure on the recordings of th e retin al activity , for v alues of τ R ranging from 2 to 5 milliseconds. The couplings d id not change much with resp ect to the val ues found with τ R = 0. N ote that th e introduction of a propagation d ela y τ D in the synaptic interactio n is straightfo rward, as long as the integration kernel remains a Dirac distribution (centered in τ D ). Bounds on th e val ues of the couplings and cu rren ts e.g. to prev ent the exponen- tial gro wth of negative in teractions with the leaking conductance can naturally be introduced through a prior distribution. As an example, assume t hat the interactions J ij take v alues in [ J − , J + ]. Then, one could maximize L ∗ − X i,j W ( J ij ) instead of th e log-lik elihoo d L ∗ alone, where W ( J ) = w 2 ( J − J − ) 2 if J < J − , 0 if J − < J < J + , w 2 ( J − J + ) 2 if J > J + and w is a large p ositiv e co eﬃcien t. W e ha ve assumed, throughout this work, th at the v alues of g and V th w ere know n. In practical situations, while the orders of magnitudes are known, the p recise v alues of these parameters should b e inferred, and could dep end on the neuron i . The inference proced ure could be modiﬁed to up date th e va lues of g i and ( V th ) i at th e same time as the synapt ic couplings J ij and the curren t I i . The num b er of parameters to infer (p er neuron) would simply increase from N to N + 2, and the running time should not increase to o muc h. 4.4 Inference from a limited neural p opulation and in the p resence of a stimulus Now ada ys, multi-ele ctro de ex periments can record a few t ens, or hundreds of neurons. T o whic h extent do the in teractions inferred from this sub- population coincide with the interactio ns one would ﬁnd from th e knowledge of the whole p opulation activity? The q uestion do es not arise in cross-correlation analysis: the correlation b etw een th e ﬁring activities of tw o neurons is obviously ind ependent of whether a third neu ron is recorded or not. H o w ever the issue must b e addressed as so on as a collective mo del for generating th e activity is assumed, such as the coup led LIF mo dels studied here. A detailed analysis suggests that the interaction b et ween a pair of neurons is n ot aﬀected by t h e activit y of other neurons distan t by more than ℓ = 300 µ m in the case of sp ontaneous activity (Cocco, Leibler and Monasson, 2009). The electro de arra y 26 should b e at least twice longer and wider than ℓ , and should b e dense enough t o capture al l the neu rons on t he recorded area . It is estimated th at about 10% of the ganglion cells are registered in the Dark experiment, compared to more t h an 80% with the denser b ut smaller electrode array u sed in the N atu ral Movie ex periment (Segev, Puchal la and Berry , 2005). It would thus b e very interesting to repeat our stu dy on other multi-electrode recordings, with suﬃciently large and d ense arra ys. T aking into account th e stimulus S in the infer ence process w ould also b e inter- esting. T o do so, w e could add a stim ulus-indu ced current, I s ( t | S ), to (1). A simple expression for this current would b e I s ( t | S ) = R t 0 dt ′ K s i ( t − t ′ ) S i ( t ′ ), where K s i is a kernel similar to the one used in generalized linear mo dels (Pillo w et al, 2008 ). The ex- pression of the current-depend en t term in the p oten tial V ( η , t ) (10) should be mo diﬁed accordingly , while the noise-dep endent t erm w ould remain unchanged. It is imp ortan t to note that the search pro cedure for contacts presented in Section 2.4 w ould remain v alid. H o w ever, the expressions of the noise co eﬃcien t, th e contact time and the du- ration of a passive contact giv en in App endix B for the cas e of a constan t current I should b e rederived and would depend on the precise temp oral stru cture of the stim ulus-indu ced curren t I s ( t | S ). Ackno wledgment: This w ork originates from a collab oration with S . Leibler, whom w e th an k for numerous and fruitful discussions. W e thank C. Barbieri for a critical reading of the manuscript. W e ac kno wledge th e hospitalit y of The R ock efeller Un i- versi ty , where th is w ork was initiated. P artial funding w as p ro vided by the Agence Nationale de la Recherc he under contract 06-JCJC-051. A Active contacts In this App endix, we justif y the prescri ptions in the search for activ e cont acts presented in Section 2.4. F or the sak e of si mplicit y we restrict to the g = 0 case (no membrane leak age); the extension to non-zero g is brieﬂy discussed in App endix B. W e consider a neuron i , and call M the num b er of spikes receiv ed by this neuron duri ng its k th int er-spike interv al [ t 0 ≡ t i,k ; t M +1 ≡ t i,k +1 ]. The arriv al times are t 1 < t 2 < . . . < t M , and the corresponding synaptic strengths ar e J 1 , J 2 , . . . , J M . T o lighten notations we hereafter omit the index i of the neuron. A.1 Case of M = 0 or 1 input spike T o understand the key notion of contact , we ﬁrst consider the si m ple case of a neuron receiving no spike duri ng the inte r-spi k e inte rv al [ t i,k = 0; t i,k +1 = T ]. The optimal noise is constant according to (7). Equation (6) then shows that the optimal potential is a li nea r function of the time, which is fully determined from the b ound ary conditions V ∗ (0) = 0 , V ∗ ( T ) = V th . W e obtain V ∗ ( t ) = V th t T and η ∗ ( t ) = C V th T − I . (33) This solution is correct since the p oten tial remains below the threshold at al l times 0 < t < T . Let us no w assume now that the neuron receives one i nput fr om another neuron, of s tr ength J 1 at time t 1 ∈ ]0; T [. The eﬀect of the i nput is a discon tinuous j ump of the p ote ntial at time t 1 and of size J 1 C , shown in Fig. 8. Rep eating the calculation ab o v e, we obtain the foll o wing expressions for the optimal p oten tial and noise V ∗ A ( t ) =  V th − J 1 C  t T + J 1 C θ ( t − t 1 ) and η ∗ A = C V th − J 1 T − I (case A) , (34) 27 0 t 1 T 0 1 V * / V th 0 t 1 T 0 1 0 t 1 T 0 1 0 t 1 T 0 1 2 η * x T/(CV th ) 0 t 1 T 0 1 2 time A B C 0 t 1 T 0 1 2 Fig. 8 Sk etc hes of the optimal p oten tials V ∗ (top) and noises η ∗ (bottom) for one neuron receiving one w eak ( A ), one strong negativ e ( B ), and one strong positive ( C ) input. The jump in the optimal noise consecutive to an active con tact is alwa ys p ositiv e. V alues of the parameters used for the ﬁgure are: I = 0, t 1 = T / 2, J 1 / ( C V th ) = . 2 ( A ), - 1.2 ( B ), 1.2 ( C ). where θ is the Hea viside function: θ ( x ) = 1 i f x > 0, 0 otherwise. This solution is sk etc hed in Fig. 8A. It is v alid when the p oten tial V ∗ A remains b elo w the threshold at all times. W e call this situation case A . As V ∗ A is a piece-wise l inear f unct ion we only need to chec k that V ∗ A ( t − 1 ) and V A ( t + 1 ) are b oth smaller than V th . The tw o conditions are fulﬁlled provided that J − ≡ − C V th T − t 1 t 1 < J 1 < J + ≡ C V th . (35) What happens when the ab o ve condition is vi olate d? Let us consider ﬁrst J 1 < J − (referred to as case B hereafter). Then V ∗ A exceeds the threshold V th befor e the input enters the neuron. T o preven t the p oten tial fr om crossing the threshold at time t 1 , the true optimal noise, η ∗ B , should be smaller than η ∗ A . But, if η ∗ B < η ∗ A , the p oten tial could not reach V th when the neuron emits its s pi k e at time T acc ording to the v ery deﬁnition of η ∗ A ! The only wa y out is that η ∗ B tak es tw o di ﬀerent v alues corresp onding to the tw o sub-interv als [0; t 1 [ and ] t 1 ; T ], which w e call, respective ly , η ∗ B, − and η ∗ B, + . W e expect η ∗ B, − < η ∗ A < η ∗ B, + . The noise can c hange v alue in t = t 1 through (9) only if the potential reac hes the threshold in t 1 . W e ﬁnd that V ∗ B ( t ) = V th t t 1 and η ∗ B, − = C V th t 1 − I (case B , 0 < t < t 1 ) , (36) from the b ounda ry conditions V ∗ (0) = 0 , V ∗ ( t − 1 ) = V th , and V ∗ B ( t ) = V th + J 1 C T − t T − t 1 and η ∗ B, + = − J 1 T − t 1 − I (case B , t 1 < t < T ) , (37) from the boundary conditions V ∗ ( t + 1 ) = V th + J 1 C , V ∗ ( T ) = V th . This solution is drawn in Fig. 8B. It is i mportant to stress that the abov e solution is based on the capabilit y of the noise to abruptly cha nge its v alue when the p oten tial touc hes the threshold in t = t 1 . A detailed study of the behavior of the noise close to such ‘conta ct points’ proving that this i s indeed the case i s p ostponed to A ppendix A.2. Finally , we turn to case C corresp onding to J 1 > J + . In this case the input is so excitatory that the noise has to b e negative to preven t the neuron from emitting a spike at a time t < t 1 . As in case B, the potential r eac hes the threshold i n t = t 1 to allow the noise to c hange its v alue after the input has en tered the neuron. W e ﬁnd V ∗ C ( t ) =  V th − J 1 C  t t 1 and η ∗ C, − = C V th − J 1 t 1 − I (case C , 0 < t < t 1 ) , (38) according to the b oundary conditions V ∗ C (0) = 0 , V ∗ C ( t − 1 ) = V th − J 1 C . Right after the spike has been receiv ed, the p oten tial has reached its threshold v alue, and will keep to this v alue un til a 28 spike is emitted at time T , hence V ∗ C ( t ) = V th and η ∗ C, + = − I (case C , t 1 < t < T ) . (39) This s ol ution is dra wn in Fig. 8C. W e now give the v alues of log-like liho ods L ∗ corresponding to the cases li sted ab o ve. The v alue of L ∗ can b e calculated fr om the kno wledge of the optimal noise η ∗ through (16). In the case of M = 0 spik e, we ﬁnd, using (33) with T = t 1 − t 0 , L ∗ ( t 0 , t 1 | I ) = − ( C V th − I ( t 1 − t 0 )) 2 2( t 1 − t 0 ) . (40) The optimal curr en t i s then inferred by maximizing L ∗ ( I ) with the result ˆ I = 1 t 1 − t 0 , which corresponds to a v anishing v alue f or the optimal noise, as exp ect ed. When M = 1 spike is receiv ed b y the neuron, the log-li k elihoo d L ∗ has three distinct expressions corresp onding to the case A, B, C discussed i n Section A.1. The resulting expression is (with t 2 = T ): L ∗ ( t 0 , t 1 , t 2 | J 1 , I ) =          − ( C V th − J 1 − I ( t 2 − t 0 )) 2 2( t 2 − t 0 ) if J − < J 1 < J + (case A) − ( C V th − I ( t 1 − t 0 )) 2 2( t 1 − t 0 ) − ( J 1 + I ( t 2 − t 1 )) 2 2( t 2 − t 1 ) if J 1 < J − (case B) − ( C V th − J 1 − I ( t 1 − t 0 )) 2 2( t 1 − t 0 ) − I 2 2 ( t 2 − t 1 ) i f J 1 > J + (case C) . (41) The l og-lik eliho od L ∗ is a contin uous and con vex function of i ts argument. The ﬁrst deriv ativ es of L ∗ are con tin uous in J − , J + , but the second deri v atives are not. A.2 Stud y of the optimal noise close to an active con tact p oint The nois e coeﬃcient η in (8) are constan t o ver the time interv al separating tw o active conta cts. The v alue of η m a y how eve r change upon the crossing of an activ e cont act. The scope of this section i s to sho w that the noise r igh t after the cont act can take any value larger than the noise immediately befor e the cont act. This m onot onicity prop ert y justiﬁes the search for the minimal noise co eﬃcient done i n (11), see App endix A.3. T o show that the noise alwa ys increases through an active contac t, we consider that the synaptic int egration is not instant aneous, but tak es place ov er a ﬁnite alb eit small time, τ s . W e th us replace the expression for the current I sy n i in (2) with I sy n i ( t ) = X j ( 6 = i ) J ij X k K ( t − t j,k ) (42) where J ij is the str eng th of the connection from neuron j onto neuron i , and K ( τ ) i s is the memory k ernel of the integration of synaptic entries (top panel in Fig. 9). W e assume that K ( τ ) v anishes f or τ < 0 and for τ > τ s where the integrat ion tim e τ s is indep end ent of the pair ( i, j ). In addition, K is p ositiv e, and its in tegral o ver the interv al [0; τ s ] is equal to unit y . W e consider the case of a single incoming spike, as in Section A.1. W e wan t to show that, in the τ s → 0 limit, the only constraint linking the v alues η ∗ − and η ∗ + of the optimal nois e, resp ectively , before and after a spik e entering at t 1 , i s η ∗ + > η ∗ − , as we ha ve f ound for a single incoming i nput in cases B and C. T o do so, we assume that the time of synaptic inte gration τ s is small but ﬁnite , and consider case B. The dynamics of V ∗ and η ∗ can be divi ded i n several steps, whose num bered are reported on Fig. 9: 1. Prior to the input, i. e. at times < t 1 , the optimal noise η ∗ − is constan t and the optimal potential V ∗ is a linear function of the tim e, with sl ope ( I + η ∗ − ) /C , as shown in Fi g. 9(left). 2. A s tr ongly negativ e input of amplitude J 1 ( < J − ) is then receiv ed by the neuron betw een times t 1 and t 1 + τ s . The deriv ative of the poten tial now decrease s with the time until it v anishes at time t c deﬁned through K ( t c − t 1 ) = k − where k − ≡ η ∗ − + I − J 1 . (43) 29 t + τ 1 s s τ η * − η * − * η + * η + t + τ 1 s 1 t t 1 t K 1 η * t t’ c c k + k − ~1/ k + k − c c V * K 1 B C t’ t c t’’ c t’’ η * V=V V=V th th C C J J t 1 3 V * 4 2 6 5 Fig. 9 Behaviors of the optimal poten tial V ∗ (middle) and noise η ∗ (bottom) close to a con tact p oin t, compared to the memory k ernel K (top). An input of total amplitude J 1 en ters the neuron during the time int erv al t 1 < t < t 1 + τ s . Left: J 1 is strongly negativ e as in Fi g. 8B; italic num ber s r efer to the steps listed in the main text. Right: J 1 is strongly posi tiv e as in Fig. 8C. See text for a detailed descri ption of the curves, of the constan ts k − , k + (43,44), and of the times t 1 , t c , t ′ c , t ′′ c , τ s . 3. If the v alue of η ∗ − is cho sen so that V ∗ ( t c ) = V th , the p oten tial tangen tially r eac hes the threshold at t c (con tact p oin t). Then, the p ote ntial r emains constan t and equal to V th . The noise obeys eqn. (9) and, therefore, increases un til i t reac hes the prescrib ed v alue, η ∗ + , at time t ′ c suc h that K ( t ′ c − t 1 ) = k + where k + ≡ η ∗ + + I − J 1 , (44) see b ottom panel in Fig. 9(left). 4. Then the p ote ntial starts decreasing fr om its threshold v alue through eqn. (6), and reache s a mini mum i n t ′′ c , s ol ution of the same equation (44) as t ′ c , see Fig. 9(top left). 5. At l ate r times the deriv ative of the p oten tial is p ositiv e from eqn. (6), and increases un til time t 1 + τ s , coinciding with the end of the synaptic integrat ion. 6. At times larger than t 1 + τ s , the p oten tial ke eps growing with a constant sl ope equal to ( I + η ∗ + ) /C . In the τ s → 0 l imit, all times t c , t ′ c , t ′′ c tend to the same v alue, that is, the time t 1 . More precisely , as the s l ope of K is of the order of τ − 2 s (in absolute v alue), and η ∗ − , η ∗ + , V ∗ ( t 1 ) are ﬁnite (= O (1)), then for τ s → 0, t 1 , t c , t ′ c diﬀer fr om eac h other by O ( τ 2 s ). Hence the cha nge in the p oten tial V ∗ betw een t ′ c and t 1 + τ s equals J 1 C + O ( τ s ). W e conclude that, for τ s → 0 the p oten tial b ecome s a discon tinuou s function of time wi th a discon tinuit y J 1 C . In addition, the noise η ∗ can also jump abruptly fr om its v alue η ∗ − at t − 1 to an y lar ger v alue η ∗ + at time t + 1 since the maximum of K tends to inﬁnit y when τ s → 0. Note that the dr a wing of Fig. 9(left) tacitly assumes that k + > k − . A hypothetic scenario wo uld be that the noise exactly compensates the synaptic input for a longer time inte rv al (including the top of K ), while the potential would remain equal to V th . In this case, the p eak v alue of the noise would b e O (1 /τ s ). T he co ntribution to the integral (16) would b e of the order of 1 /τ s and w ould di v erge in the τ s → 0 limit. Hence this p ossibility is precluded. 30 The abov e discussion is straight forwardly exte nded to case C. The optimal p ote ntial and noise are sketc hed in Fig. 9(right) . Note that the con tact i n terv al spr ead s beyond [ t ′ c , t c ] in this case. In the generic case of more than one incoming spikes, the con tact i n terv al is restricted to [ t ′ c ; t c ] as in case B. The noise can also discontin uously c hange from its v al ue η ∗ − < − I b efore the con tact to any larger v alue, η ∗ + , af ter the con tact. A.3 Case of M ≥ 2 incoming spikes W e no w consider the general case of M input spikes. Let V 0 = 0 , m 0 = 1 be, resp ectiv ely , the initial v alue of the p oten tial and the i ndex of the ﬁrst input spike. W e deﬁne the piece-wis e linear f unct ion solution of (6) for a constan t noise η , V ( η , t, t 0 ) = V 0 + I + η C ( t − t 0 ) + M X m = m 0 J m C θ ( t − t m ) . (45) W e are l ooking for the small est v alue of the noise co eﬃcien t η capable of bri nging the p oten tial V ( η , t, t 0 ) f rom its initial v alue V ( η , t 0 , t 0 ) = 0 to the threshold. The cont act tim e, t c , coincides with an entering spike , i. e . t c = t m ∗ for some m ∗ ≥ 1. If m ∗ = M + 1 then the optimal potent ial is V ( η ∗ , t, t 0 ) throughout the i nter-spike int erv al [ t 0 ; t M +1 ], and the problem is solv ed. If m ∗ ≤ M , t m ∗ is the ﬁrst activ e cont act p oin t of the p oten tial. η ∗ and V ( η ∗ , t ) are, resp ectiv ely , the optimal noise and p oten tial on the interv al [ t 0 , t m ∗ ]. The correctness of the abov e statemen t can be established using a proof by con tradiction. – assume that the optimal noise, η opt , is small er than η ∗ on some sub-i n terv al of [ t 0 ; t c ]. Remark that the potent ial V i n (45) is an increasing function of the noise, η ′ > η = ⇒ V ( η ′ , t, t 0 ) > V ( η , t, t 0 ) , (46) for all t > t 0 . By virtue of (46) and the deﬁnition of η ∗ , V ( η opt , t, t 0 ) cannot touch the threshold at any time so the noise is constan t throughout the interv al [ t 0 ; t c ]. Hence no activ e contact can take place at tim e t c . As η ∗ is the mini mal v alue of the noise which can drive the potential i n to con tact with the threshold ov er [ t 0 ; t M +1 ], we conclude that V ( η opt , t, t 0 ) cannot cr oss the threshold at any time ≤ t M +1 . The neuron can therefore not spike at time t M +1 . – con versely , suppose that the optimal noise is equal to η α > η ∗ on the inte rv al [ t 0 ; t m α ] with 1 ≤ m α < m ∗ , and tak es another v alue on the interv al [ t m α ; t c ] 7 . As the change of noise can tak e place only through an active contact , and the c hange is necessaril y p ositiv e (Section A. 2), we ha ve η β > η α . A ppl yi ng (46) to the interv al [ t m α ; t c ], w e ha v e V ( η β , t c , t m α ) > V ( η α , t c , t m α ) . (47) Adding the v alue of the opt imal potent ial in t m α to b oth members of the previous in- equalit y , we ﬁnd V ∗ ( t c ) = V ( η β , t c , t m α ) + V ( η α , t m α , t 0 ) > V ( η α , t c , t m α ) + V ( η α , t m α , t 0 ) = V ( η α , t c , t 0 ) > V ( η ∗ , t c , t 0 ) (48) where the last inequality comes from (46). But, by deﬁnition of η ∗ , V ( η ∗ , t c , t 0 ) = V th . Hence, we ﬁnd that the optimal p ote ntial in t c is abov e threshold, which cannot be true. The optimal noise and p oten tial on the remaini ng part [ t c ; t M +1 ] of the inte r-s pi k e interv al can b e determined iteratively . W e replace t 0 with t m ∗ and V 0 with V th if J m ∗ > 0 or V th + J m ∗ C if J m ∗ < 0 in (45), and lo ok for the lo west noise pr oducing a ne w con tact point ov er the int erv al [ t m ∗ , t M +1 ]. The procedure is rep eated until the whole interv al is exhausted. This wa y an increasing sequence of noise v alues is obtained, eac h corresponding to the slop e of the optimal poten tial b et w een tw o successive contact points. 7 The case of three or a higher num b er of v alues for the noise can be handled exactly i n the same w ay . 31 t t t t 0 2 1 V th V ∆ ∆ c c (2) (1) t c Fig. 10 Sketc h of the optimal p ote ntial close to a passi v e contac t sta rting at tim e t c . The duration of the passiv e con tact is ∆ c ( ℓ ), where ℓ is the index of time t ℓ corresponding to the next active contac t. The p otentials corr esponding to the hy p othesis ℓ = 1 and ℓ = 2 are sho wn with the dashed and f ull curves resp ectiv ely . B Passiv e contacts When the m embrane leaking conductance is non zero, some change hav e to b e brought to the abov e calculation of the optimal noise and potential. Fir st, in the absence of inputs, the noise is no longer constan t, but rather it i s an exponen tially increasing (in absolute v alue) function of the time (8 ). Sim ilarly , the p oten tial V ∗ itself is not a linear function of the time as in (48), but i s a linear combination of exp( ± t/τ ) with appropriate co eﬃcien ts, see (10). The m ain conclusion of App endix A still holds: the diﬀerence betw een the noise v alues just after and b efore an active cont act p oin t, coinciding wi th a synaptic input, is alwa ys p ositiv e (Fig. 2A). Consequen tly , the pro cedure of Section 2.4, i.e. the iterative searc h for the active con tact p oin ts and the mi ni mal noise coeﬃcient η ∗ , deﬁned through (11), remains unchang ed. Note that some care must be taken to translate the statemen t ab out the growth of the noise to the v alues of the noise co eﬃcien ts. Consider f or instance tw o successive con tact times, t and t ′ , and call η , η ′ the corresp onding noise coeﬃcients. That the noise is larger at time t ′ than at time t implies that η × exp(( t ′ − t ) /τ ) < η ′ , but do es not imply that η ′ is larger than η 8 . There exists, ho wev er, a ma jor diﬀerence b et we en th e g = 0 and g 6 = 0 cases. When g > 0, the optimal p oten tial i s not guaran teed to b e a monotonous function of the time, as sho wn in Fig. 10. F or given v alues of g , I , and of the times and the amplitudes of the synaptic inputs, the optimal p oten tial V ∗ may touc h the threshold at an in termediate time, t c . Such a situation is called passiv e con tact. It is imp ortan t to note that the v alue of the optimal nois e during a passiv e contact is, according to eqn. (9), equal to g V th − I . As the optimal noise is a monotonous f unct ion of the time betw een tw o activ e con tact s, see eqn (8), the v alue g V th − I can b e crossed at m ost once: there is at most one passive conta ct in betw een tw o successiv e activ e ones. T o be more precise, there are at most A + 1 passive contact s in an inter-spik e int erv al with A activ e con tacts. T o decide the existence of a passive contact in an int erv al [ t 0 ; t M +1 ], we l ook for a solution of th e t wo coupled equat ions expressing that the optimal p oten tial touc hes the threshold without crossing it, V ∗ ( η p , t c ) = V th and ∂ V ∗ ∂ t ( η, t c ) = 0 . (49) The solutions of these equations give the noise co eﬃcien t η p and the contac t time t c at which the optimal p oten tial reaches the threshold v alue (Fig. 10). The solution can be calculated analytically , wi th the f ollo wing result. Let us call V 0 the v alue of the p oten tial of the neuron at time t + 0 . F or each m ≤ M w e deﬁne V m = V 0 + X ℓ : t 0

Original Paper

Loading high-quality paper...

Comments & Academic Discussion

Loading comments...

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment