Instantaneous Non-Linear Processing by Pulse-Coupled Threshold Units

Contemporary theory of spiking neuronal networks is based on the linear response of the integrate-and-fire neuron model derived in the diffusion limit. We find that for non-zero synaptic weights, the response to transient inputs differs qualitatively from this approximation. The response is instantaneous rather than exhibiting low-pass characteristics, non-linearly dependent on the input amplitude, asymmetric for excitation and inhibition, and is promoted by a characteristic level of synaptic background noise. We show that at threshold the probability density of the potential drops to zero within the range of one synaptic weight and explain how this shapes the response. The novel mechanism is exhibited on the network level and is a generic property of pulse-coupled networks of threshold units.

Dynamical systems driven by uctuations are ubiquitous models in solid state physics, quantum optics, chemical physics, circuit theory, neural networks and laser physics.Absorbing boundaries are especially interesting for diusion over a potential step, escape from a particle trap, or neuron models [1].Approximating uctuations by Gaussian white noise enables analytical solutions by Fokker-Planck theory [2,3].For non-Gaussian noise, however, the treatment of appropriate boundary conditions gains utmost importance [4].The advantage of a transport description to solve the mean rst passage time problem has previously been demonstrated [5].In this line of thought, here we present a novel hybrid theory that augments the classical diusion approximation by an approximate boundary condition for nite jump Poisson noise.Exact results have so far only been obtained for the case of exponentially distributed jumps amplitudes [6].We apply our theory to the leaky integrate-and-re neuron model [1], a noise-driven threshold system widely used to uncover the mechanisms governing the dynamics of recurrent neuronal networks.An incoming synaptic event causes a nite jump of the membrane potential which decays back exponentially.
The neuron res a spike if the membrane potential reaches a threshold.This simplicity renders the model analytically tractable, ecient to simulate with precise spike times [7], and yet it captures the gross features of neural dynamics [8].The commonly pursued approach is to linearize this non-linear input-output unit around a given level of background activity and to treat deviations in the input as small perturbations.This technique has been applied successfully to characterize the phase diagram of randomly connected recurrent networks by a mean-eld approach [9], to quantify the transmission of correlated inputs by pairs of neurons [10,11] and to understand the interplay of spike-timing dependent learning rules with neural dynamics [12].For Gaussian white noise input current the linear response kernel is known exactly [9]: It constitutes a low-pass lter for signals modulating the mean [13]; only modulations of the uctuations are transmitted immediately [14].In this Letter we show how our novel hybrid approach allows the analytical prediction of a genuinely instantaneous non-linear response never reported so far.Poisson noise with nite synaptic jumps even enhances this response.

CREMENTS
We consider a rst order stochastic dierential equation driven by point events from two Poisson processes with rates ν + and ν − .Each incoming event causes a nite jump J k = J + > 0 for an increasing event and J k = J − < 0 for a decreasing event (1) We follow the notation in [15] and employ the Kramers-Moyal expansion with the innitesi- The rst and second innitesimal moment evaluate to A 1 (x) = f (x) + µ and A 2 = σ 2 , where we introduced the The time evolution of the probability density p(x, t) is then governed by the Kramers-Moyal expansion, which we truncate after the second term to obtain the Fokker-Planck equation where S denotes the probability ux operator.In the presence of an absorbing boundary at θ, we need to determine the resulting boundary condition for the stationary solution of (2).Without loss of generality, we assume an absorbing boundary at θ at the right end of the domain.Stationarity implies a constant ux φ through the system.Rescaling the density by the ux as q(x) = φ −1 p(x) results in the linear inhomogeneous dierential equation of rst order ( The ux over the boundary has two contributions, the deterministic drift and the positive stochastic jumps crossing the boundary with [x] + = {x for x > 0; 0 else}.To evaluate the integral in (5), for small J + θ − x we develop q(x) into a Taylor series around θ. To this end, we solve (3) for q (x) = − 2 It is easy to see by induction, that all higher derivatives q (n) can be written in the form q (n) (x) = c n (x)+d n (x)q(x) whose coecients obey the recurrence relation Inserting the Taylor series into (5) and performing the integration results in which is the probability mass moved across threshold by a perturbation of size s and hence also quanties the instantaneous response of the system.Dividing (4) by φ we solve it for q(θ) to obtain the Dirichlet boundary condition If J + is small compared to the length scale on which the probability density function varies, the probability density near the threshold is well approximated by a Taylor polynomial of low degree; throughout this letter, we truncate (7) and (8) at n = 3.

APPLICATION TO THE LEAKY INTEGRATE-AND-FIRE NEURON
We now apply the theory to a leaky integrate-and-re neuron [1] with membrane time constant τ and resistance R receiving excitatory and inhibitory synaptic inputs, as they occur in balanced neural networks [16].We model the input current i(t) by point events t k ∈ {incoming spikes}, drawn from homogeneous Poisson processes with rates ν e and ν i , respectively.The membrane potential is governed by the dierential equation τ dV dt (t) = −V (t) + Ri(t).An excitatory spike causes a jump of the membrane potential by J k = w, an inhibitory spike by J k = −gw, so Ri(t) = τ J k δ(t − t k ) + Ri 0 , where i 0 is a constant background current.Whenever V reaches the threshold V θ , the neuron emits a spike and the membrane potential is reset to V r , where it remains clamped for the absolute refractory time τ r .With µ def = τ w(ν e − gν i ) and σ 2 def = τ w 2 (ν e + g 2 ν i ), we choose the natural units u = t/τ and y = (V − Ri 0 − µ)/σ to obtain A 1 (y) = −y and A 2 = 1.The probability ux operator (2) is then given as S = −y − 1 2 ∂ ∂y .For the stationary solution q(y) = ν −1 0 p(y) the probability ux between reset y r and threshold y θ must equal the ring rate ν 0 , and is zero else, so The equilibrium solution q(y) = Aq h (y) + q p (y) of ( 9) is a linear superposition of the homogeneous solution q h (y) = e −y 2 and the particular solution q p (y) = 2e −y 2 y θ max(yr,y) e u 2 du, chosen to be continuous at y r and to vanish at y θ .The constant A is determined from (8) as A = q(y θ )/q h (y θ ).We obtain the mean ring rate ν 0 from the normalization condition of the density 1 = ν 0 τ y θ −∞ q(y) dy + ν 0 τ r , where ν 0 τ r is the fraction of neurons which are currently refractory Figure 1 shows the equilibrium solution near the threshold obtained by direct simulation to agree much better with our analytic approximation than with the theory for Gaussian white noise input.Close to reset V r = 0, the oscillatory deviations with periodicity w are due to the higher occupation probability for voltages that are integer multiples of a synaptic jump away from reset.They wash out due to coupling of adjacent voltages by the deterministic drift as one moves away from reset.
We now proceed to obtain the response of the ring rate ν to an additional δ-shaped input current τ s R δ(t) fed into the neuron.This input causes a jump s of the membrane potential at t = 0 and (2) suggests to treat it as a time dependent perturbation of the mean input Since the perturbation has a at spectrum, up to linear order in s the spectrum of the excess rate is fs (z) = sτ H(z) + O(s 2 ), where H(z) is the linear transfer function with respect to perturbing µ at Laplace frequency z.In particular, P r (s) = fs (0).As H(0) is the DC susceptibility of the system, we can express it up to linear order as H(0) = ∂ν 0 ∂µ .Hence, FIG. 1: Finite synaptic potentials distort the stationary membrane potential density P (V ).A Black: Incoming spike rates ν e = 29800 Hz, ν i = 5950 Hz (corresponding to µ = 12 mV and σ = 5 mV).
Histogram binned with ∆V = 0.01 mV.Gray: novel approximation B Magnication of A around spike threshold.Light gray: solution in diusion limit of [9].C,D Density for supra-threshold current Ri 0 = 20mV and incoming rates ν e = 95050Hz, ν i = 22262.5Hz (corresponding to µ = 12 mV and σ = 9.5 mV).Other parameters and gray code as in A,B.
We also take into account the dependence of A on µ to calculate dν 0 dµ from (10) and obtain Figure 2D shows the integral response to be in good agreement with the linear approximation.The integral response in the diusion limit is almost identical.
The instantaneous response of the ring rate to an impulse-like perturbation can be quantied without further approximation.The perturbation shifts the probability density by s so that neurons with V ∈ [V θ − s, V θ ] immediately re.This results in the nite ring probability P inst (s) of the single neuron within innitesimal time (5), which is zero for s < 0. This instantaneous response has several interesting properties: For small s it can be approximated in terms of the value and the slope of the membrane potential distribution below the threshold (using (7) for n ≤ 2), so it has a linear and a quadratic contribution in    (7).Histograms binned with h = 0.1 ms.C Medium gray curve: instantaneous response P inst (7) as a function of s for nite weights w = 0.1 mV.Black dots: Direct simulation.Light gray curve: Diusion limit of (7).Medium gray dots: Direct simulation of diusion limit with temporal resolution 10 −4 ms.D Gray curve: integral response for nite weights (11).Black dots: direct simulation.Gray dots: direct simulation for Gaussian white noise background input.Simulated data averaged over 2.5 • 10 8 (s = 0.1 mV) . . .2.5 • 10 6 (s = 1.0 mV) perturbation events.Other parameters as in Figure 1A.
s. Figure 2A shows a typical response of the ring rate to a perturbation.The peak value for a positive perturbation agrees well with the analytic approximation (5) (Figure 2C).
Due to the expansive nature of the instantaneous response (Figure 2C) its relative contribution to the integral response increases with s.For realistic synaptic weights ≤ 1 mV the contribution reaches 30 percent.Replacing the background input by Gaussian white noise, and using the boundary condition q(y θ ) = 0 in (7) yields a smaller instantaneous response (Figure 2C) which for positive s still exhibits a quadratic, but no linear, dependence.The integral response, however, does not change signicantly (Figure 2D).An example network in which the linear non-instantaneous response cancels completely and the instantaneous response becomes dominant is shown in Figure 3A.At t = 0 two populations of neurons simultaneously receive a perturbation of size s and −s respectively.This activity may, for example, originate from a third pool of synchronous excitatory and inhibitory neurons.The   pooled linear ring rate response of the former two populations is zero.The instantaneous response, however, causes a very brief overshoot at t = 0 (Figure 3B). Figure 3C reveals that the response returns to baseline within 0.3ms.Figure 3D shows the non-linear dependence of peak height P net (s) = 1 2 P inst (s) on s.In this Letter we present an extension to the diusion approximation with nite but small increments.Our theory describes the probability density evolution as a diusion on the length scale σ determined by the uctuations, but we take the quantization of the noise near an absorbing boundary into account, leading to a non-zero Dirichlet condition.This hybrid approach enables us to nd analytical approximations hitherto unknown for pure jump models.In particular we accurately quantify the instantaneous contribution to the escape rate in response to a perturbation.There is a formal similarity to the ndings of [17] for large perturbations.Applied to the integrate-and-re neuron with Gaussian white noise input, we quantify a quadratic non-linear instantaneous ring rate response not captured by the existing linear theory [13,14].Finite jumps in the noise qualitatively change and en- We consider a dynamic variable governed by a dierential equation of rst order with Poissonian input causing nite jumps J k = J + > 0 for increasing events arriving with rate ν + and J k = J − < 0 for decreasing events of rate (1) We follow the notation in [4] and employ the Kramers Moyal expansion with the innitesimal moments The rst innitesimal moment evaluates to where κ + , and κ − denote the Poisson distributed number of positive and negative incoming events in the time interval of size h, respectively.We used lim h→0 1 h κ + = ν + , lim h→0 1 h κ − = ν − and introduced the shorthand where the rst and second term in the second last row vanish for h → 0 and for the independent individual Poisson distributed random numbers κ + and κ − it holds So the innitesimal variance is independent of x The time evolution of the density is then governed by the Kramers-Moyal expansion We truncate this series after the second term to obtain the Fokker-Planck equation with the probability ux operator

II. ABSORBING BOUNDARY CONDITION
Without loss of generality, we assume an absorbing bondary at θ that has to be crossed from below.We need to determine the proper boundary condition for the stationary solution L FP p(x) = 0 at the threshold.Stationarity implies a constant ux φ through the system.Rescaling the solution by the ux as q(x) = φ −1 p(x) results in the linear inhomogeneous dierential equation of rst order On the other hand, the ux crossing the boundary from left has two contributions, the deterministic drift and the positive stochastic jumps crossing the boundary with [x] + = {x for x > 0; 0 else}.For small J + θ − x we express q(x) near the threshold θ by a Taylor series.To this end, we solve (7) It is easy to see by induction, that all higher derivatives q (n) can be written in the form q (n) (x) = c n (x) + d n (x)q(x) which yields the recurrence relations Inserting the Taylor series of q around θ q into ( 9) and integrating results in Inserted into (8) and dividing by φ we solve for q(θ) to obtain the Dirichlet boundary condition We can make use of the fact, that typically J + is small compared to the length scale on which the probability density function varies.So the probability density near the threshold is well approximated by a Taylor polynomial of low degree; throughout this letter, we truncate ( 12) and ( 13) at n = 3.

III. LEAKY INTEGRATE-AND-FIRE NEURON
We now apply the theory to the leaky integrate-and-re neuron [6] with membrane time constant τ and resistance R receiving excitatory and inhibitory synaptic inputs, as they occur in balanced neural networks [7].We model the input current i(t) by point events t k ∈ {incoming spikes}, drawn from homogeneous Poisson processes with rates ν e and ν i , respectively.The membrane potential is governed by the dierential equation an excitatory spike causes a jump of the membrane potential by J k = w, an inhibitory by J k = −gw.Whenever V reaches the threshold V θ , the neuron emits a spike and the membrane potential is reset to V r where it remains clamped for the absolute refractory time τ r .i 0 is a constant input current.With µ def = wτ (ν e − gν i ) and σ def = w 2 τ (ν e + g 2 ν i ) we choose the natural units u = t/τ and y = V −Ri0−µ σ proposed by ( 2) and ( 3) to obtain the form ( 1) where the events u k = t k /τ arrive with rate τ ν e (excitatory) and τ ν i , respectively.Consequently, in these coordinates A 1 (y) = −y and A 2 = 1, so the probability ux operator (6) acting on the rescaled probability density q(y) = ν −1 is given as For the stationary solution q of (6), the probability ux between reset and threshold must equal the stationary ring rate The equilibrium solution of ( 15) is a linear superposition of the homogeneous solution q h satisfying the dierential equation Sq h = 0 given by ( 14), which is and the particular solution q p .The latter is chosen to be continuous at y r , to vanish at y θ , and to fulll the inhomogeneous dierential equation ( 15) q p = −2yq p − 2H(y − y r ) (with H denoting the Heaviside function) =2e −y 2 y θ max(yr,y) The solution which fullls the boundary conditions is the superposition where the constant A is determined as A = q(y θ )q −1 h (y θ ).Using the recurrence (10) we determine the sequence of the c n and d n as The Dirichlet value at the threshold of given by (13) as Since the explicit solution ( 16) of q h is known and since we have chosen q p (y θ ) = 0, it follows that d n (y θ )e −y 2 θ is the n-th derivative q (n) h (y θ ) so that we can replace the sum in the denominator by Hence the boundary value can be calculated as A comparison of this approximation to the full solution is shown in Figure 1C.Within the range of several synaptic weights w = 0.1 mV used here, the approximation is nearly perfect.

IV. RESPONSE PROPERTIES
In this section, we calculate the response of an integrate-and-re neuron to an additional δ-shaped input current.We follow two approaches here: First, we quantify the integral response in linear perturbation theory.Then we go beyond linear response and determine the non-linear instantaneous response to a δ-shaped input current.
Treating the neuron in linear response theory, it is sucient to determine the response of the neuron's ring rate with respect to an input δ-current at t = 0.Such a perturbation causes a jump of the membrane potential and can be treated as a time dependent perturbation of the mean input by µ(t) = µ + sτ δ(t) as can be seen from (1).We call the time-dependent response of the ring rateν s (t) and its excess over baseline f s (t) = ν s (t) − ν 0 .The integral response is dened as P r (s) := ∞ 0 f s (t) dt.We denote by fs (z) = ∞ 0 e zt f s (t) dt its Laplace transform.The spectrum of the perturbation is at, so to linear order in s the spectrum of the response is fs (z) = sτ H(z) + O(s 2 ), where H(z) is the linear transfer function of the system with respect to perturbing µ.In particular, P r (s) = fs (0).Since H(0) is the linear DC susceptibility of the system (changing µ innitesimally slowly), we can equate H(0) = ∂ν0 ∂µ .Hence we arrive at FIG. 1: A analytic ring rate compared to simulation depending on uctuations σ of the input.Black: Analytic solution (20), light gray: ring rate in diusion limit [1,5].Parameters τ = 20 ms, V θ = 15 mV, Vr = 0, w = 0.1 mV, g = 4, τr = 1 ms.νe, νi chosen to realize mean input µ = 12 mV and the uctuation σ as given by the abscissa.B Probability density near threshold.Black dots: direct simulation, solid gray line: analytic result P (V ) = ν0τ /σQ((V − µ)/σ) given by (19).Dark gray solid line: Polynomial of degree 3 the density around V θ using (11).Parameters as in A, but with νe = 29800 Hz, νi = 5950 Hz (corresponding to µ = 12 mV and σ = 5 mV).C Analytic ring rate compared to simulation depending on the size of synaptic jumps w for xed σ = 5 mV and µ = [7, 9.5, 12, 14] mV (from bottom to top).Gray code as in A. D Membrane potential distribution near threshold for µ = 12 mV and σ = 5 mV for dierent w = 0.05 mV (black) w = 0.1 mV (dark gray), w = 0.2 mV (gray), w = 0.5 mV (light gray).
We calculate dν0 dµ , from (20), where we need to take into account also the dependence of A on µ.With dν −1 dµ , we only need to calculate dν −1 0 dµ .Additionally, we note that dy dµ = −1 σ .So we nd To calculate ∂A ∂µ we use ( 8) The last term's integrand unfolds to .
Taken together the derivative unfolds to To calculate the instantaneous response of the neuron to an impulse-like perturbation we have to leave linear response theory.The incoming spike has a direct contribution on the ring of the neuron, due to the probability density just below threshold which is shifted over the threshold.This is illustrated in Figure 3A.This direct contribution manifests itself in a nite response probability given by P inst (s) (9), which is zero for s < 0. The instantaneous response has several interesting properties: For small strength s it can be expressed by the value and the rst derivative of the membrane potential distribution at the threshold.The expression has a linear and a quadratic term in the perturbation s, so the neuron will inherently respond non-linearly to a positive perturbation s > 0, whereas it will not respond to a negative perturbation s < 0, as shown in Figure 2B.
The integral response (21) as well as the instantaneous response (9) both exhibit stochastic resonance, an optimal level of synaptic background noise σ that alleviates the transient.Figure 2A shows this noise level to be at about σ = 3 mV for the integral response.The responses to positive and negative perturbations are symmetric and the maximum is relatively broad.The instantaneous response in Figure 2B displays a pronounced peak at a similar value for σ.This non-linear response only exists for positive perturbations and is zero for negative ones.Stochastic resonance has been reported for the linear response to sinusoidal periodic stimulation [3] and also for non-periodic signals that are slow compared to the neuron's dynamics an adiabatic approximation reveals stochastic resonance [2].In contrast to the latter work cited, the rate transient observed in our work is the instantaneous response to a very fast (Dirac δ) perturbation.
For those neurons, which after the perturbation have not red instantaneously, the membrane potential is shifted by s. Figure 3B shows a sketch of the situation just after the perturbation has been received.The emission rate ν + of the neuron in a small time bin (0, h] of size h τ following the perturbation can be approximated by similar arguments as P inst .Knowing the distribution J(γ) = ∞ i,j=0 P(i|ν e h)P(j|ν i h)δ w(i−gj),γ of the membrane potential The upper trace is the response to a positive perturbation of magnitude s = 0.5mV, the lower to a perturbation of s = −0.5mV.
Black dots: direct simulation.Gray solid line: analytic result using (24).B Instantaneous response depending on synaptic background noise σ.The upper trace is the response to a positive impulse of weight s = 0.5 mV, the lower one to a negative of s = −0.5 mV.Black dots: direct simulation.Gray solid line: analytic peak response Pinst using (12).Simulations averaged over N = 1000 neurons for T = 100 s.Incoming spike rates νe, νi chosen to realize µ = 11.5 mV and σ on the abscissa for synaptic weights w = 0.1 mV and g = 4.

A B
FIG. 3: A Sketch illustrating the instantaneous response due to a perturbation of size s: The probability mass in the shaded area is shifted above threshold with probability 1 and produces the instantaneous response.Due to the shape of the probability density near the threshold, this area has a linear and a quadratic increasing contribution in s.B For all neurons in the population which have not been carried over the threshold directly by the perturbation, the membrane potential is shifted by s.
During a small time interval h τ after the perturbation, the density remains shifted.This causes the ring rate to increase approximately proportional to −s ∂P ∂V (V θ ) due to jumps γ occurring with probability P(γ) by synaptic input.
jumps γ due to synaptic input in a time step h, where P(k|µ) = µ k k! e −µ is the Poisson distribution, the instantaneous rate of the neuron can be approximated by the probability per time moved above the threshold

2 ]FIG. 2 :
FIG.2: Instantaneous ring rate response to perturbation.A Black: Response to perturbation s = 0.5 mV at t = 0, gray: s = −0.5 mV.B Magnication of A. Black cross: analytic peak response P inst h(7).Histograms binned with h = 0.1 ms.C Medium gray curve: instantaneous response P inst(7)  as a function of s for nite weights w = 0.1 mV.Black dots: Direct simulation.Light gray

3 ]FIG. 3 :
FIG. 3: The non-linear response to perturbations is exhibited on the network level.A Two identical populations of N = 1000 neurons each, receive uncorrelated background input (light gray spikes).At t = 0 the neurons simultaneously receive an additional input of size s in the upper and −s in the lower population (symbolized by black single spike).B Pooled response of the populations normalized by the number of neurons.C Magnication of B. D Instantaneous response P net (s) (black dots: direct simulation, gray curve: analytical result) as a function of s.Other parameters as in Figure 2A.

FIG. 2 :
FIG.2:A Integral response of the ring rate (21) of an integrate-and-re neuron as a function of synaptic background noise σ.The upper trace is the response to a positive perturbation of magnitude s = 0.5mV, the lower to a perturbation of s = −0.5mV.