Moments of the interspike interval density

-0.5 0 0.5 1 0

0.5 1 P 0(v)

-0.5 0 0.5 1

v -0.5 0 0.5 1

τref = 0

τref = 0.1

τref = 0.5

A B C

Figure 3.5.: Influence of the refractory periodτ_refon the stationary density of an LIF neuron.

As, for the parameters used here, the threshold can only be crossed during "+" dynamics, all reset trajectories initially drift to the right forτ_ref = 0 (A). Thus, the density jumps to higher values at vR. As τ_ref is increased (B, C), the fraction of trajectories for which the noise has switched increases, so that already forτ_ref =0.1 (B), a majority of trajectories initially drifts to the left, leading to a downward jump in the density. Parameters:k+=3.3,k−=1.2,σ=3,µ= 2.5,vR=0,vT =1.

while the latter is an effect of its stationary distribution, specifically the fact that it takes discrete values. An additional difference lies in the effect that the refractory period has on the probability density: For white noise driving, it only enters through a reduced firing rate. Here, by contrast, the refractory period has a qualitative impact on the probability density aroundvR, as illustrated in Fig. 3.5.

3.4. Moments of the interspike interval density

In order to characterize neural spike train statistics, it is useful to derive expressions not only for the rate, but also for higher moments of the ISI density. A frequently used measure to quantify the regularity of spiking is the coefficient of variation (CV), which is defined in terms of the first two moments of the ISI density,

CV = q

hT²i − hTi²

hTi ^{. (3.61)}

In this section, we derive recursive relations that allow to express thenth moment of the ISI density of IF neurons driven by a DMP in terms of quadratures.

The calculation of ISI moments corresponds to solving a FPT problem for the mem-brane voltage; specifically, the ISI density corresponds to the density of first-passage times from the reset voltagev_Rto the thresholdv_T.

For the case of one-dimensional systems driven by Gaussian white noise, recursive relations for the FPT moments have been known for a long time (Siegert, 1951). An alter-native derivation was given by Lindner (2004b); here, we follow this approach and adapt it to DMP-driven IF neurons.

The starting point of our derivation is the observation that the FPT density corresponds to the time-dependent flux across the threshold of an ensemble in which all trajectories

Chapter 3. IF neurons driven by dichotomous noise

start at v_R at time t = τ_ref and where, importantly, one has made sure (through appro-priate boundary conditions) that no probability can flow back from above the threshold (which would mean that reentering trajectories are counted more than once).

We start by writing the master equation eqs. (3.18, 3.19) in terms of the total flux, J(v,t) = J+(v,t) +J−(v,t), and the flux in "-" state,J−(v,t). The latter is a useful choice to obtain simple boundary conditions, because we know that J−(v,t)needs to vanish at FPs of the "-" dynamics. After some simplification, the resulting equations read

∂_vJ(v,t) =−∂_t

We multiply both sides of eqs. (3.62, 3.63) bytⁿand integrate them overtfrom−_∞to

∞. For this, it is convenient to introduce the abbreviations J_n(v):= The term involving a time derivative can be integrated by parts. For eq. (3.63), for exam-ple, one has

Here, we have used that there is no flux before t = τ_ref and that all fluxes vanish for t → ∞, as eventually all trajectories will have crossed the threshold. The integration by parts in eq. (3.62) is carried out analogously.

One obtains a system of ODEs, d

3.4. Moments of the interspike interval density

Eq. (3.68) can be directly integrated and eq. (3.69) can be solved by variation of con-stants. Keeping in mind the necessary treatment of FPs in the "-" dynamics (see the pre-vious section), the voltage axis needs to be divided into N intervals delimited by the threshold voltagev_T, FPs in the "-" dynamics, and the lowest attainable voltage v− (as illustrated for a QIF in Fig. 3.4). Denoting byi(_v)the interval containing the voltagev, and byl_i andr_i the lower and upper boundaries of theith interval, the solutions can be given as

Here, we have already satisfied the boundary condition forJ_n(v): The total flux needs to vanish at the lowest attainable voltagev−=l₁, which impliesJ_n(l₁) =0.

The boundary conditions for Mⁱ_n(v) warrant a more detailed discussion. There are three kinds of BCs:

3. Finally, if the ith interval is delimited below (above) by an unstable FP in the

"-" dynamics, settingc_i = l_i (c_i = r_i) is needed for the proper behavior at the FP (vanishing instead of diverging fluxJ−(v_U,t)).

The three BCs thus lead to the same choice inc_i as in the calculation of the stationary density above. Again, they consistently cover all cases: For example, if f(v_T)−σ > 0, then the "no backflow" condition is automatically fulfilled and cannot be used to deter-minec_N. However, in this casec_N is still uniquely determined by one of the other BCs:

Either f(v)−σ does not have zero crossings betweenv_T andv_R, in which case the left interval boundary isv^→_R , at which f(v_R)−σ > 0 and BC 2 demandsc_N = v^→_R = l_N. Or there is a zero crossing, in which case the lower interval boundary is necessarily an un-stable FP and BC 3 demandsc_N =v_U =l_N. Note thatc_iis always one of the two interval boundaries.

Our aim was to obtain an expression for the nth moment of the ISI density. As we

Chapter 3. IF neurons driven by dichotomous noise

Figure 3.6.: Firing rate and CV of LIF and QIF neurons over k− for different values of k+. Theory (lines) is compared to simulation results (symbols). Parameters:µ = 0.5,σ = 1,v_R = 0,v_T =1,τ_ref=0 (LIF),µ=0.5,σ=1,v_R=−10,v_T =10,τ_ref=0 (QIF).

have argued, the ISI density corresponds to the flux across the threshold with boundary conditions as described above. To calculate the nth moment, we may thus recursively calculateJ_n(v)and evaluate it at the threshold,

hTⁿi=

Z∞ 0

dT TⁿJ(v_T,T) =J_n(v_T). (3.72)

From eq. (3.70), one sees directly that J₀(v_T) = 1, consistent with the requirement that the FPT density be normalized.

Finally, eqs. (3.70, 3.71) provide a quick way to calculateα. Recall thatαis the fraction of trajectories that cross the threshold in "+" state. Thus,

α=

Z∞ 0

dt J+(_vT,t) =J₀(_vT)− M₀(_vT)_{. (3.73)}

This can be solved forα(in general, the r.h.s. also depends onα!), yielding the expression given in eq. (3.59).