Abstract

Spike trains with negative interspike interval (ISI) correlations, in which long/short ISIs are more likely followed by short/long ISIs, are common in many neurons. They can be described by stochastic models with a spike-triggered adaptation variable. We analyze a phenomenon in these models where such statistically dependent ISI sequences arise in tandem with quasi-statistically independent and identically distributed (quasi-IID) adaptation variable sequences. The sequences of adaptation states and resulting ISIs are linked by a nonlinear decorrelating transformation. We establish general conditions on a family of stochastic spiking models that guarantee this quasi-IID property and establish bounds on the resulting baseline ISI correlations. Inputs that elicit weak firing rate changes in samples with many spikes are known to be more detectible when negative ISI correlations are present because they reduce spike count variance; this defines a variance-reduced firing rate coding benchmark. We performed a Fisher information analysis on these adapting models exhibiting ISI correlations to show that a spike pattern code based on the quasi-IID property achieves the upper bound of detection performance, surpassing rate codes with the same mean rate—including the variance-reduced rate code benchmark—by 20% to 30%. The information loss in rate codes arises because the benefits of reduced spike count variance cannot compensate for the lower firing rate gain due to adaptation. Since adaptation states have similar dynamics to synaptic responses, the quasi-IID decorrelation transformation of the spike train is plausibly implemented by downstream neurons through matched postsynaptic kinetics. This provides an explanation for observed coding performance in sensory systems that cannot be accounted for by rate coding, for example, at the detection threshold where rate changes can be insignificant.

1  Introduction

How a neuron's spiking pattern is produced by its intrinsic dynamics and how the resulting spike patterns encode information are distinct yet connected central questions in neuroscience (Stevens & Zador, 1995; London, Roth, Beeren, Haüsser, & Latham, 2010; Jacobs et al., 2009; Panzeri, Harvey, Piasini, Latham, & Fellin, 2017). Spike patterns arise from the confluent effects of synaptic input and intracellular firing dynamics that preserve some memory of past spike times. Point process models provide a means to assess these combined effects on both the descriptive statistics and coding properties of spiking patterns (Kass & Ventura, 2001; Keat, Reinagel, Reid, & Meister, 2001; Truccolo, Eden, Fellows, Donoghue, & Brown, 2005; Paninski, 2004; Citi, Ba, Brown, & Barbieri, 2014). Including a spike-triggered intracellular adaptation variable into a point process model is one such memory mechanism by which past spike events induce negative feedback that reduce the probability of current spike emissions (Kass & Ventura, 2001; Truccolo et al., 2005; Lüdtke & Nelson, 2006; Muller, Buesing, Schemmel, & Meier, 2007; Nesse, Maler, & Longtin, 2010; Farkhooi, Muller, & Nawrot, 2011). This adaptation can produce interspike interval (ISI) temporal dependencies in which spikes emitted in short succession (i.e., short ISIs) are commonly followed by longer ISIs and vice versa. Such negative serial correlations are commonly observed in many cell types and organisms (Perkel, Gerstein, & Moore, 1967; Nawrot et al., 2007; Farkhooi, Strube-Bloss, & Nawrot, 2009; Berry & Meister, 1998). They reduce spike count variability and enhance the detection of rate changes relative to comparable renewal interval sequences (Chacron, Longtin, & Maler, 2001; Chacron, Lindner, & Longtin, 2004; Chacron, Maler, & Bastian, 2005; Lüdtke & Nelson, 2006; Nawrot et al., 2007; Maimon & Assad, 2009; Farkhooi et al., 2009, 2011; Avila-Akerberg & Chacron, 2011; Nikitin, Stocks, & Bulsara, 2012).

Here, novel properties of stochastic models with spike-triggered adaptation that produce temporal ISI correlations are of interest. A previous study by our group (Nesse et al., 2010) demonstrated that a specific hybrid model, involving both a hazard function and adaptation dynamics, can exhibit statistically dependent successive ISIs; but paradoxically, the spike-triggered adaptation states exhibit probabilistic quasi-independence with identical distribution (quasi-IID). That work established that quasi-IID adaptation states could be utilized in lieu of correlated ISI sequences to represent the spike pattern of a cell. Adaptation dynamics thus amount to a nonlinear decorrelation transformation of the ISI sequence. The results in Nesse et al. (2010) established this quasi-IID property for a specific stochastic spiking model and confirmed its existence for a conductance-based model as well, but a more general set of sufficient conditions guaranteeing the property was not made.

We specify sufficient conditions on a general class of stochastic spiking models that guarantee the quasi-IID property and give several examples. We also establish that such quasi-IID adaptation processes necessarily induce negative serial ISI correlations and establish bounds on that correlation. Furthermore, because spike-triggered adaptation processes are dynamically similar to that of synaptic transmission, such representations of the spike pattern may be plausibly communicated postsynaptically if there are matched pre- and postsynaptic kinetics (Pfister, Dayan, & Lengyel, 2010; Khanbabaie, Nesse, Longtin, & Maler, 2010; Avila-Akerberg & Chacron, 2011; Marcoux, Clarke, Nesse, Longtin, & Maler, 2016), and our letter addresses this possibility.

An important related question is whether this quasi-IID representation of the spike pattern could provide additional information about an input that is not available in firing rate statistics (Gussin, Benda, & Maler, 2007). As stated above, numerous studies have established that negative serial ISI correlations reduce the variance of firing rate, making firing rate changes more statistically discriminable relative to comparable renewal ISI spike trains (Chacron et al., 2001, 2004; Chacron, Maler, & Bastian, 2005; Lüdtke & Nelson, 2006; Nawrot et al., 2007; Maimon & Assad, 2009; Farkhooi et al., 2009, 2011; Avila-Akerberg & Chacron, 2011; Nikitin et al., 2012). This variance-reduction phenomenon defines the rate coding performance benchmark for spike trains with negative serial interval correlations. The quasi-IID representation, because it contains a complete representation of the spike train in an IID format, provides the means to rigorously assess the statistical discrimination performance of serially correlated ISI sequences in response to weak signal inputs.

Here we show that spike pattern coding that utilizes the quasi-IID adaptation states is statistically efficient for weak signal discrimination. Further, we calculate that it exceeds the previous coding benchmark by proving an inequality between the Fisher information of each code when large numbers of spikes are utilized for this weak signal discrimination. We then show that the spike pattern Fisher information exceeds the Fisher information for all rate codes, including variance-reduced rate codes, by 20% to 30%, in this general class of adapting models. We find that this reduction in rate code Fisher information arises because, while temporal ISI correlations reduce firing rate variance, the negative feedback of the adaptation also reduces firing rate gain in response to input signal changes, offsetting some of the gain in information from the variance reduction. The additional information contained in the spike pattern, beyond that found in the firing rate, may account for the behavioral performance observed in some sensory systems where firing rate coding is insufficient to account for observed weak signal discrimination performance (Knudsen, 1974; Gussin et al., 2007; Jacobs et al., 2009; Jung, Longtin, & Maler, 2016; Fettiplace & Hackney, 2006; Heil & Peterson, 2015). Indeed this question is one of the main motivations of our study, along with the related question of what determines the extreme detection capabilities seen in various senses across species.

The letter is organized as follows. We begin by defining the general category of adapting point process models we study. We then establish model conditions that guarantee the quasi-IID property theorem and provide a proof. Three illustrative example models are then discussed. We next explore consequences that the theorem has on spike statistics and make a case for the biophysical plausibility of such a coding scheme. Finally, we prove that the pattern coding that arises from the quasi-IID property exceeds that of rate codes.

Let $λ(H):R→[0,∞)$ be a nonincreasing and Lipschitz continuous spike emission rate function. The adaptation function $H(η)$: $R→R$ controls the emission rate, where $η$ is a special “adaptation time” variable. Adaptation time $η$ increases at the same rate in lockstep with “clock” time $t$ for intervals in which no spiking occurs (see below for more details). The probability of spiking in the time interval $[η,η+dη]$ is $λ(H(η))dη+O(dη2)$. The adaptation function $H$ is differentiable and monotonic decreasing in the $η$ variable, making $H(η)$ invertible. The nonincreasing assumption of $λ(H)$ dictates that the decline of $H$ elicits a nondecreasing probability of spiking as adaptation time increases. That $λ(H(η))$ is nondecreasing in $η$ is a reasonably nonrestrictive assumption that simply expresses that spiking is increasingly likely (or same) as adaptation $H$ declines.

Adaptation time $η$ increases as clock time $t$ increases for intervals in which no spiking occurs. What distinguishes $η$ from $t$ is that both $η$ and $H$, which are in one-to-one correspondence, are instantaneously assigned to new values when a spike occurs. The reassignment of $H$ is through a spike-triggered activation function $g(H):R→R$ (see Figure 1A). That is, when $H$ is viewed as a function of clock time $t$$H(η(t))=H(t)$—the graph has discontinuous jumps at the spike times. We require $g$ to be monotonic increasing, differentiable, and thus invertible. Here we use an asterisk to denote the adaptation state and adaptation time value just before a spike, while variables without an asterisk denote values immediately after a spike. Adaptation declines monotonically after an activation. Let $H(η*)=h*$ be an instantaneous adaptation state. By inverting the function $H$, one can obtain an adaptation time $η*$ associated with this adaptation state. This specially constructed adaptation time is in some sense internal to the adaptation process: it jumps “backward” in a time frame associated with adaptation whenever the stochastic model spikes in the real external clock time. Thus, just prior to the $i$th real-time spike, the adaptation time equals $ηi*$. The spike-triggered adaptation by definition increments the adaptation state from $hi*$ just prior to the spike to the new state $hi$, to which is associated a new adaptation time $ηi$, because $H(η)$ is invertible. The function $g(h*)=h$ thus governs the jump to produce the postspike adaptation state of adaptation.
Figure 1:

Schematic of the stochastic model with spike-triggered adaptation and its output. (A) The $H(η)$ function, the invertible function $f$ that maps $ηi*→ηi$, and the invertible function $g$ that maps $hi*→hi$ ($i=0,1$). Here we assume the model produced a spike at the moment when the adaptation state was $h0*$, which corresponds to adaptation time $η0*$. The adaptation state is updated to $h0=g(h0*)$ (red activation pathway), which corresponds to an adaptation time $η0$. The function $f$ thus maps $η0*→η0$, and $g$ maps $h0*→h0$. The adaptation then decays, and the next firing of the stochastic model occurs when the adaptation state equals $h1*$ with corresponding adaptation time $η1*$. The spike-triggered adaptation then iterates the maps $f$ and $g$ again to obtain, respectively, $η1$ and $h1$ (blue activation pathway). (B) An example $λ(·)$ function (black) on the right ordinate, with area-under-curve left of the $ηε$ threshold less than $ε$ (shaded); and the survival functions $P(η|η0)$ (red) and $P(η|η1)$ (blue) on the left ordinate. Note that the starting points of the survival functions differ and determine the length of the effective refractory period. But once $η$ surpasses the $ε$-refractory threshold $ηε$, the survival functions will be very similar. (C) Independent random draws from the $q∞$-distribution (in green on the vertical axis) of the model in example 3 produce $H(t)$ traces with postspike $hn$-peaks at each spike time. (D) The resulting ISIs have statistical dependence, illustrated by the negative serial autocorrelation $ρΔt(k)$ of successive ISIs (i.e., at lag $k=1$); nonsuccessive ISIs have negligible correlation.

Figure 1:

Schematic of the stochastic model with spike-triggered adaptation and its output. (A) The $H(η)$ function, the invertible function $f$ that maps $ηi*→ηi$, and the invertible function $g$ that maps $hi*→hi$ ($i=0,1$). Here we assume the model produced a spike at the moment when the adaptation state was $h0*$, which corresponds to adaptation time $η0*$. The adaptation state is updated to $h0=g(h0*)$ (red activation pathway), which corresponds to an adaptation time $η0$. The function $f$ thus maps $η0*→η0$, and $g$ maps $h0*→h0$. The adaptation then decays, and the next firing of the stochastic model occurs when the adaptation state equals $h1*$ with corresponding adaptation time $η1*$. The spike-triggered adaptation then iterates the maps $f$ and $g$ again to obtain, respectively, $η1$ and $h1$ (blue activation pathway). (B) An example $λ(·)$ function (black) on the right ordinate, with area-under-curve left of the $ηε$ threshold less than $ε$ (shaded); and the survival functions $P(η|η0)$ (red) and $P(η|η1)$ (blue) on the left ordinate. Note that the starting points of the survival functions differ and determine the length of the effective refractory period. But once $η$ surpasses the $ε$-refractory threshold $ηε$, the survival functions will be very similar. (C) Independent random draws from the $q∞$-distribution (in green on the vertical axis) of the model in example 3 produce $H(t)$ traces with postspike $hn$-peaks at each spike time. (D) The resulting ISIs have statistical dependence, illustrated by the negative serial autocorrelation $ρΔt(k)$ of successive ISIs (i.e., at lag $k=1$); nonsuccessive ISIs have negligible correlation.

Upon spiking, the postspike adaptation time $η$ is found to depend on the prespike adaptation time through the composition function $f=H-1∘g∘H:R→R$, which is monotonic increasing, differentiable, and invertible; it represents the mapping of the adaptation time variable: $f(η*)=η$. As we will show, this formulation of the time variable $η$ with discontinuous jumps at each spike is necessary for the statistical analysis of the point process generated by the stochastic spiking model. These $η$-discontinuities at spike events allow $h$ and $η$ to be in a one-to-one correspondence over every interspike interval, such that when $h$ is activated via the $g$-function, the $η$-variable is correspondingly reset via the $f$-function.

The emission function $λ$, adaptation $H$, and activation $g$ (or equivalently $f$) together constitute the adapting point process model: ${λ,H,g}$. Spikes are emitted stochastically, followed by a postspike-triggered activation of adaptation (with a corresponding postspike time resetting $η$), and then on to the next spike, forming a cycle of spike, activate, spike again, and so on. This cycle is illustrated in Figure 1A. Let $η0*$ be a spike time associated with the level of adaptation $H(η0*)=h0*$ just prior to the clock time at which the model generates a first “reference” spike. Then $η0$ is the postspike reset adaptation time. The spike triggers an activation $η0=f(η0*)$, or, alternatively, $h0=g(h0*)$, which defines the red circuit. Eventually the next spike event occurs at $η1*$ (blue circuit) and then resets again to $η1$, and so on, producing tandem sequences ${η0,η1,η2,…}$ and ${h0,h1,h2,…}$. In what follows we focus our presentation on the postspike-activated $h$-values and reset $η$-values (nonasterisked) in order to simplify the approach, bearing in mind that the prespike values are obtained through the inverse transformations: $η*=f-1(η)$ and $h*=g-1(h)$.

Let $P(η|η0)$ be the proportion of an infinite ensemble that has not yet produced a subsequent spike event in $η≥η0$ time thereafter. The function $P(η|η0)$ is termed the survival function, where $P(η0|η0)=1$ and $P(η|η0)≤1$ is nonincreasing in $η$. The function $P$ is defined uniquely as the solution to the differential equation
$dPdη=-λ(H(η))P⟹P(η|η0)=e-K(η|η0),$
(2.1)
where $K$ in equation 2.1 is termed the potential function:
$K(η|η0)≡∫η0ηλ(H(ξ))dξ.$
(2.2)
The survival function, equation 2.1, is shown in Figure 1B, and a similar $P(η|η1)$ is also depicted.
We wish to use equation 2.1 to characterize the statistics of the resulting spike trains. The $n$th interspike interval (ISI) is defined as $Δtn=ηn+1*-ηn=f-1(ηn+1)-ηn≥0$. The ISI-based survival function is then $P(Δtn+ηn|ηn)$. The stochastic process $ηn$ and its analog $hn$ can both reconstruct the ISI sequence through the following expressions:
$Δtn=f-1(ηn+1)-ηn=f-1(H-1(hn+1))-H-1(hn)≡A(hn+1)-B(hn).$
(2.3)
We have introduced functions $A≡f-1∘H-1$ and $B≡H-1$ in equation 2.3 for notational simplicity. The probability densities for $η$ and $h$ are derived below and are utilized to study the statistics of $Δt$ using equation 2.3.
From the survival function $P(η|η0)$, the probability density $p(η|η0)$ for $η$ conditioned on $η0$ is found through differentiation,
$p(η|η0)=∂∂η(1-P)=λ(H(η))e-K(η|η0),$
(2.4)
which we will use to construct a stochastic process ${ηn}$ of postspike start times. Upon spiking at $η1*$, adaptation is activated to $g∘H(η1*)=h1$, with reset time-variable $η1=f(η1*)$ (see Figure 1B). The sequence ${ηn}$ is determined by the Markov transition operator $R$ defined as
$R(ηn+1|ηn)=p(f-1(ηn+1)|ηn)ddηf-1(ηn+1),$
(2.5)
where $ddηf-1$ is the Jacobian. If $rn(η)$ is a density describing an ensemble of activated $η$-variables on the $n$th spike (we assume the densities are sufficiently smooth), then the density on the $(n+1)$st is given by the standard Chapman-Kolmogorov equation: $rn+1(η)=∫R(η|ξ)rn(ξ)dξ.$ At each stage of the sequence, an alternative representation can be made in terms of the postspike adaptation states of the system by introducing the change of variables to the $hn$-space, where $h=g∘H(η)≡w(η)$ represents the peak activation of the adaptation state variable immediately after a spike. A Markov transition operator $Q$ for the $hn$ following this change of variables is given by
$Q(hn+1|hn)=p(w-1(hn+1)|H-1(hn))ddhw-1(hn+1),$
(2.6)
yielding the sequence of distributions $qn+1(h)=∫Q(h|ξ)qn(ξ)dξ.$ In the case of both $R$ and $Q$ operators, equations 2.5 and 2.6, Ruelle-Perron-Frobenius operator theory guarantees the existence of invariant asymptotic densities $limn→∞rn=r∞$ and $limn→∞qn=q∞$, where $r∞(η)$ and $q∞(h)$ are fixed points of the Chapman-Kolmogorov equation: $r∞(η)=∫R(η|ξ)r∞(ξ)dξ$, and $q∞(h)=∫Q(h|ξ)q∞(ξ)dξ$. Our previous work (Nesse et al., 2010) established that for a specific model, there are cases in which the $hn$-sequence exhibits near statistical independence. This near-independence occurs when the Markov operators, equations 2.5 and 2.6, exhibit very weak conditional dependence. The next section establishes a more general result.

The objective of this section is to establish the conditions for which $R(ηn+1|ηn)∼R(ηn+1)$ and $Q(hn+1|hn)∼Q(hn+1)$, that is, that both $η$ and $h$ are independent, or quasi-independent, of the previous value, amounting to an effectively trivial Markov process. This implies that the asymptotic distributions are obtained from a single application of the Markov operator: $r∞(η)∼R(η)$, and $q∞(h)∼Q(h)$.

Stated simply, an effectively trivial Markov process is guaranteed if the activation $g∘H(ηn*)=hn$ is always large enough, and in turn $ηn$ goes sufficiently below a threshold (see the vertical dashed line in Figures 1A and 1B) that the emission function $λ$ is minuscule for some time after a spike. In other words, the adaptation causes a significant relative refractory period. Consider the following theorem and proof:

Theorem 1

(the Quasi-IID Adaptation Property). Suppose that for a given tolerance $0≤ε≪1$ there exists a particular $ηε$ such that $K(η|ηε)=∫ηηελ(H(ξ))dξ<ε$ for all $η≤ηε$. Further suppose that $f(η)≤ηε$ for all $η$, and $η0<ηε$.

Then the survival function can be expressed as
$P(η|η0)=e-K(η|η0),η≤ηεP(η)[1+(e-K(ηε|η0)-1)],η>ηε$
(3.1)
where in the first case, $η≤ηε$ of equation 3.1, the survival function is $ε$-close to unity, that is, $(1-e-K(η|η0))<ε$ since the emission function is so small. In the second case, $η>ηε$ of equation 3.1, the term $(e-K(ηε|η0)-1)$ is $ε$-bounded as well: $|e-K(ηε|η0)-1|<ε$, and there is a function $P(η)$ that is independent of $η0$ and is $ε$-close to the true survival function,
$|P(η)-P(η|η0)|<ε,$
(3.2)
for all $η$ and $η0$ such that $η0<ηε<η$.

The proof of theorem 1 is found in the appendix.

Theorem 1 implies that for sufficiently small $ε$, $P(η)≈P(η|η0)$. We will use equality ($=$) in lieu of approximation ($≈$) henceforth and define $P(η)≡P(η|η0)=P(η)$, for notational simplicity. We also regard the quasi-IID variables as IID in subsequent analyses. Also note that $η≤ηε$ and $h≥hε$ define the upper and lower domain boundaries for the respective variables. We use $K$ as the potential function in both $η$ and $h$ variables and define distributions $p∞(η)=K'e-K$ where the prime means differentiation with respect to the relevant variable; the analogous postspike activated form $r∞$ is similarly defined through a change of variables. On the boundaries of the respective domains of $h$ or $η$, the potential function has the limits $K→0$ and $K→∞$ so that the probability function $P=e-K∈[0,1]$, as required.

Now we establish the consequences of quasi-independent adaptation states on spike statistics. We have shown that a spiking model with adaptation that satisfies the hypotheses of theorem 1 yields both the adaptation states $hn$ and adaptation times $ηn$ as quasi-independent random variable sequences, with distributions $q(h)$ or $r(η)$, respectively. Note that by using the activation functions, the ISIs are expressible in terms of the adaptation states or adaptation times:
$Δtn=A(hn+1)-B(hn)=f-1(ηn+1)-ηn=ηn+1*-ηn.$
(3.3)

Equation 3.3 defines a nonlinear transformation that expresses ISI sequences in terms of quasi-IID adaptation states $hn$ or quasi-IID adaptation times $ηn$.

We now turn to the correlation properties of these ISIs. The ISI density can be obtained from the spike time density $p∞(η*)=ddη(1-P(η*))$ and related postspike adaptation time density $r∞(η)$, obtained from $p∞$ through the change of variables $η=f(η*)$. The difference between two independent random draws, one from $p∞$ and the other from $r∞$ defined by equation 3.3, yields an ISI. The distribution of ISIs is then computed via the convolution:
$pISI(Δt)=∫p∞(Δt+η)r∞(η)dη.$
(3.4)

Only subsequent ISIs are statistically dependent, and they are negatively correlated, as we now show. Recall from equation 3.3 or 2.3 that ISIs are expressed as a difference of functions of two IID variables $hn$ and $hn+1$. This means that only subsequent ISIs are negatively correlated, as is common in many neuron types (Farkhooi et al., 2009; Nawrot et al., 2007). We state this fact as a corollary:

Corollary 1.

If a point process model satisfies the IID adaptation property, then the ISI sequence is serially correlated, meaning that subsequent ISIs will be statistically dependent; however, nonsubsequent ISIs $Δtn$ and $Δtk$, where $k≠n±1$, are quasi-IID and therefore uncorrelated.

The verification is simple: serial correlation results because any subsequent pair of ISIs $Δtn$ and $Δtn+1$ depends on a common adaptation variable, namely, $hn+1$; nonsubsequent ISI pairs are functions of distinct, and thus independent, adaptation variables.

Another important question is the allowable correlation range for subsequent ISIs. Let $Δtn$ and $Δtn+1$ be subsequent ISIs. The correlation range is as follows.

Corollary 2.
If a point process model satisfies the IID adaptation property, then the correlation between subsequent ISIs $Corr(Δtn,Δtn+1)$ is:
$-12≤Corr(Δtn,Δtn+1)<0.$
(3.5)

A short proof is given in the appendix.

Finally, the mean ISI, variance, and serial covariance $Cov(Δtn,Δtn+1)$ are readily computable by taking expectations with respect to the quasi-IID distribution $q∞(h)$ and using the ISI formula, equation 3.3. For notational simplicity, we set $A(hn+1)=An+1$ and $B(hn)=Bn$. By assuming the $hn$-sequence is IID, we take expectations with respect to $q∞(h)dh$ to obtain
$〈Δt〉=〈A〉-〈B〉,$
(3.6)
$Var(Δt)=Var(A)+Var(B),$
(3.7)
$Cov(Δtn,Δtn+1)=〈A〉〈B〉-〈AB〉.$
(3.8)
Note that we are able to drop the subscripts $n$ in the above equations after simplification. Of course, $Corr(Δtn,Δtn+1)$ is the ratio of equations 3.8 and 3.7.

In sum, this section has established that equation 3.3 constitutes a nonlinear decorrelation transformation between ISI sequences and adaptation states and times. As we shall show in section 6, we will use this decomposition of ISIs into IID variables to perform an analysis of spike train coding performance. But first, we illustrate some examples of adapting point process models in the following section.

4  Examples of Quasi-IID Adapting Point Process Models and Stochastic Simulations

The following three examples start from simple and less tunable to more physiologically realistic and tunable. The first example was chosen because the distributions can all be computed in terms of elementary functions, which facilitates the deduction of its behavior directly from the equations. The other two examples involve nonelementary functions but are easily depicted graphically.

For each example, we evaluate analytically derived distributions and compare them to distributions from Monte Carlo simulations. Those simulations involved incrementing by very small time steps $Δη$ and computing the probability of spiking in each time step to first order $λ(H(η))dη$, and comparing it to a random number generator draw. If a spike was determined to occur, $η$ and $H$ were activated accordingly. Many thousands of spikes were recorded to generate empirical histograms for ISIs and $h$-data. Distributions of $η$ and $η*$ are not shown because they can be obtained from $h$ through variable changes.

Importantly, Monte Carlo simulations make no assumptions on the statistics of $hn$-data, particularly the quasi-IID property. Therefore, a good agreement between the empirical distribution of $h$-data and the presumptive IID asymptotic distribution $q∞(h)$ means that the quasi-IID property is a valid description of the true $hn$-process for the $ε$-values chosen.

4.1  Example 1: Threshold Spike Emission Function, Exponential Adaptation, and Additive Activation

The emission function
$λ(H)=(1-H),H≤10,H>1$
(4.1)
features a threshold at $H=1$ above which there is zero spike emission. The adaptation decays exponentially: $H(η)=exp(-η)$. The threshold activation is additive: $g(H)=H+1$. Adaptation time activation is then $f(η*)=-ln(exp(-η*)+1)=η$ and is paired with its inverse, $f-1(η)=-ln(exp(-η)-1)=η*$.
The adaptation increases by 1 when a spike is emitted and consequently silences spiking for a time interval until $H(η)$ decays sufficiently. The starting time point from which spiking probability becomes nonzero is $ηε=0$. We set $ε=0$ and $P(η|0)=P(η)$. The potential function then is
$K(η|0)=∫0ηλ(H(ξ))dξ=∫0η(1-e-ξ)dξ=η+e-η-1.$
(4.2)
The survival function and density are, respectively,
$P(η)=e-(η+e-η-1),η≥0,$
(4.3)
$p(η)=(1-e-η)e-(η+e-η-1),η≥0.$
(4.4)
The limiting $h$-distribution is
$q∞(h)=(2-h)e-(h-2),h∈(1,2],$
(4.5)
which is obtained by using the activation transformation $h=e-η+1$ as a change from $η$ to $h$ in equation 4.4. Figure 2A shows the relevant densities for $η$, $η*$, (blue $r∞$, and red $p∞$, respectively) and the ISI distribution ($pISI$ black), which is defined in equation 3.4. The ISI autocorrelation $ρΔt(k)=Corr(Δtn,Δtn+k)$, where $k$ is the lag, is shown in Figure 2D.
Figure 2:

Examples 1 (A), 2 (B), and 3 (C) of adapting point process models that exhibit quasi-independent adaptation sequences and statistically dependent ISI sequences. Top panels show the steady-state densities $p∞(η*)$, $r∞(η)$ of the adaptation time variables and the ISIs, that is, $pISI(Δt)$. Bottom panels show the densities $q∞(h)$ of the postspike adaptation state $h$ for the three models. (D) ISI autocorrelations function $ρΔt(k)$ shown for a few lags for all examples. (C) For Example 3, the parameters are $α=2$, $τ=10$, $Λ=0.2$, and $β=50$.

Figure 2:

Examples 1 (A), 2 (B), and 3 (C) of adapting point process models that exhibit quasi-independent adaptation sequences and statistically dependent ISI sequences. Top panels show the steady-state densities $p∞(η*)$, $r∞(η)$ of the adaptation time variables and the ISIs, that is, $pISI(Δt)$. Bottom panels show the densities $q∞(h)$ of the postspike adaptation state $h$ for the three models. (D) ISI autocorrelations function $ρΔt(k)$ shown for a few lags for all examples. (C) For Example 3, the parameters are $α=2$, $τ=10$, $Λ=0.2$, and $β=50$.

Let $λ(H)=e-H$, $H(η)=e-η$, and $g(H)=H+A$. The additive term $A$ in the activation must be large enough to guarantee theorem 1's strong activation hypothesis. In this case, $A=-ln(ε)$ makes it so. This $A$ is the minimum sufficient activation, which occurs at time $ηε$, so $A=H=e-ηε$ implies $ηε=-ln(-ln(ε))$ and $hε=-ln(ε)$. The potential function is
$K(η|ηε)=∫ηεηe-e-ξdξ.$
(4.6)
Let $u=e-ξ$ be a change of variables to the activation ($h$) space, so $-du/u=dξ$ makes the potential function $K(h*|hε)=∫h*hεe-uudu$. Note that $hε=-ln(ε)$ is sufficiently large so that the integrand $e-uu$ is vanishingly small there, and the integration limit $hε$ can be replaced with $∞$, making the potential function equal to the exponential integral $Ei(·)$ function:
$K(h*|hε)=K(h*)=Ei(h*)≡∫h*∞e-uudu.$
(4.7)
Changing to the activation postspike variable $h=h*+A$, we obtain the density function $q∞(h)=Ei'(h-A)e-Ei(h-A)$. The associated $η$, $η*$ densities that are obtainable by changes of variables are depicted in Figure 2B.

4.3  Example 3: Exponential Emission, Exponential Adaptation, with Saturating Activation

Let $λ(H)=αe-βH$ and $H(η)=e-ητ$ be similar to example 2 but incorporating tunable parameters $α$, $β$, and $τ$. The following model was explored previously to a more limited extent in a previous work by our group (Nesse et al., 2010). The parameter $β$ sets the strength of $H$-dependency on emission rate, $α$ sets the overall spike emission rate, and $τ$ sets the characteristic timescale of adaptation decay. The activation function imposes an upper bound on the adaptation state that it saturates to for very high firing rates: $g(H)=Λ(1-H)+H=Λ+(1-Λ)H$, where $Λ∈(0,1)$. The saturation dictates that $g(H)∈(Λ,1)$, so for all $η$, and $H(η)=h∈(0,1)$. The parameter $Λ$ is the minimum additive activation.

This saturating activation model comprising $g(H)$ and $H(t)$ can be derived from the aggregate current from many stochastically activating and inactivating binary-state ion channels. Let the $i$th channel open/closed state be a binary stochastic process $Xi=0,1$ defining an ion channel, for $i=1,2,…,M$, for large $M$. Let $H=1M∑iXi∈[0,1]$ be the ensemble average. Each state $Xi=0,1$ undergoes fast state transitions. During the spike event, a fraction $Λ∈(0,1)$ of the zero-states $Xi=0$ will transition to $Xi=1$. Then $H$ is the fraction already activated and $1-H$ is the fraction not activated prior to a spike. Upon spiking, the transition $H→H+Λ(1-H)$ occurs, which is the activation function we have defined already. In between the spike-triggered state transitions, there is a constant rate of transition back to the zero-state with timescale $τ$ ($Xi→0$) according to $τdHdt=-H$, resulting in a continuous exponential decay of $H$ in the mean field (for very large $M$).

For a given $ε$, all four parameters $α$, $β$, $τ$, and $Λ$ must be set appropriately to satisfy the $ε$-hypotheses of theorem 1. The parameter $Λ$ anchors our analysis. Let $Λ=hε=e-ηετ$, which is inverted to get $ηε=-τln(Λ)$. For $η0<ηε$, we have
$K(ηε|η0)=α∫η0ηεe-βH(ξ)dξ.$
(4.8)
Performing the change of variables $u=βH(ξ)=βe-ξτ$ and $-τduu=dξ$ yields $u0=βH(η0)<β$ and $uε=βH(ηε)=βΛ$, so the potential function can be bounded by
$K(ηε|η0)=ατ∫βΛu0e-uudu<ατ∫βΛβe-uudu<ατe-βΛβΛ(β(1-Λ))$
(4.9)
$=ατ(1-Λ)Λe-βΛ<ε.$
(4.10)
If the inequality in equation 4.10 is satisfied, so is theorem 1. Recall that $β$ is the coupling parameter between adaptation and emission rate—a larger $β$ makes the adaptation $H$ more effective at reducing the emission rate. Of course, one can fix any three parameters, and it would determine the fourth parameter in equation 4.10. For example, given a fixed $Λ$, $τ$, and $α$, the inequality in equation 4.10 requires that $β$ be sufficiently large. Alternatively, $Λ$ must be sufficiently large and the product $ατ$ must be sufficiently small, all other things being constant. In the appendix, we perform an extensive parameter scan of the example 3 model.
To derive the distribution $q∞(h)$, we proceed as we did for example 2: let $u=βe-ξ/τ$ so $-τdu/u=dξ$ and $h*=e-η*/τ$, so
$K(h*|hε)=ατ∫βh*βhεe-uudu.$
(4.11)
And we know $h=Λ+(1-Λ)h*$, so $h*=h-Λ1-Λ$. Moreover, $βhε≥βΛ$ is sufficiently large as per equation 4.10, so it can be replaced with $∞$ in the integral without significant change, so, $K(h|hε)→K(h)$:
$K(h)=ατ∫βh-Λ1-Λ∞e-uudu=ατEiβh-Λ1-Λ$
(4.12)
so that $q∞(h)=K'(h)exp(-K(h))$, where $K'$ is the derivative with respect to $h$. The associated $η$, $η*$ densities, each derivable by a change of variable, are depicted in Figure 2C. We used parameters $α=2$, $β=50$, $τ=10$, and $Λ=0.2$.
Embedded in example 3 are convenient expressions to relate the $hn$-sequence to the ISIs. We wish to use them in the next section, so we now make them explicit. The $hn$-sequence map can be written as a function of the ISIs $Δtn$,
$hn+1=Λ+hn(1-Λ)e-Δtnτ,$
(4.13)
and it can also be inverted as
$Δtn=τlnhn(1-Λ)hn+1-Λ$
(4.14)
$=-τln(hn+1-Λ)+τln(hn(1-Λ))]=A(hn+1)-B(hn),$
(4.15)
where the $A$ and $B$ functions are clearly observed in equation 4.15.

In this letter, we focus entirely on parameter regimes in which $ε$ is suitably small, in which point process models are in an effective IID parameter regime. Consequently, in what follows, we drop the “$∞$” on the asymptotic distributions $q∞→q$ and assume the quasi-IID process is IID. To understand more about the statistical properties of the point process models under a broader range of parameters in which there is a loss of the quasi-IID property, see the prior study by our group (Nesse et al., 2010).

5  Biological Plausibility of Encoding ISI Sequences by Adaptation States: Postsynaptic Proxies

The nonlinear transformation from the IID sequence of adaptation states, or equivalently adaptation times, to correlated ISI's is an interesting operation that simplifies data analysis. But is such a representation of the data biophysically realized and plausibly transmitted synaptically? Postsynaptic coding schemes have been proposed in which presynaptic dynamics are transmitted through matched kinetics with postsynaptic dynamics (Pfister et al., 2010; Avila-Akerberg & Chacron, 2011; Lüdtke & Nelson, 2006; Marcoux et al., 2016). In our example models presented in section 4, we have utilized an instantaneously activating and exponentially decaying adaptation variable in which the dynamics of spike-activated adaptation are similar to synaptic dynamics.

In example 3, we derived a model for adaptation based on the average of independent stochastic processes $Xi∈{0,1}$ for $i=1,2,…,M$, for large $M$ that undergo spike-activating and ISI-inactivating state transitions. Note that we can similarly define a postsynaptic proxy $S$, consisting of an average of a distinct set of independent set state variables $Yi∈{0,1}$ and for $i=1,2,…,M$ (large $M$), so that $S=1M∑iYi$. Here, $Yi$ could define the conductance states of ion channels or bound-unbound states of postsynaptic receptors.

If the independent $Xi$ and $Yi$ state variables have the same kinetics, then upon spiking, the state transition of each process yields an activation function
$g(S)=Λ+(1-Λ)S,$
(5.1)
where $Λ$ is the activation fraction as in example 3. In the interspike time interval, there is a inactivation transition rate $1/τ$ that yields exponential decay in the $S$-variable, $τdSdt=-S$, yielding exponential decline from peak $sn$ values, $S(t)=sne-(t-tn)/τ$, for $t∈[tn,tn+1)$. The $S$-dynamics are identical in form to the $H$ variable defined in example 3, where $sn$ are the postsynaptic proxy peaks on each spike, similar to the $hn$ sequence peaks of the adaptation process (see Figure 3).
Figure 3:

Both adaptation and synaptic transmission can form a decorrelating transformation of the correlated ISIs. The presynaptic intracellular adaptation code $hn$ is communicated through synaptic transmission to yield a postsynaptic proxy sequence, provided the pre- and postsynaptic kinetics are matched. Here, the $hn$ are represented postsynaptically by a quickly activating neurotransmitter binding event, followed by a slower unbinding/inactivating decay, yielding the postsynaptic proxy $S(t)$ with peak values $sn$. When the proxy is initialized in the left panel so that $sn-2=hn-2$, then $sn=hn$ always with matched kinetics (sequence 1), but if initialized differently $sn-2≠hn-2$, proxy sequences 2 and 3 illustrate the quick convergence to $hn$ in a few spikes.

Figure 3:

Both adaptation and synaptic transmission can form a decorrelating transformation of the correlated ISIs. The presynaptic intracellular adaptation code $hn$ is communicated through synaptic transmission to yield a postsynaptic proxy sequence, provided the pre- and postsynaptic kinetics are matched. Here, the $hn$ are represented postsynaptically by a quickly activating neurotransmitter binding event, followed by a slower unbinding/inactivating decay, yielding the postsynaptic proxy $S(t)$ with peak values $sn$. When the proxy is initialized in the left panel so that $sn-2=hn-2$, then $sn=hn$ always with matched kinetics (sequence 1), but if initialized differently $sn-2≠hn-2$, proxy sequences 2 and 3 illustrate the quick convergence to $hn$ in a few spikes.

Note that while $H$ and $S$ have the same kinetics and are activated at the same spike times $tn$, it is not automatically true that $hn=sn$. Particularly if the $H$ and $S$ are initialized differently, then $hn≠sn$ in general. Recall, however, that with exponential decay, from a spike time we can express the proxy as $S(t)=sne-(t-tn)/τ$. Also note that the ISI is $tn+1-tn=Δtn$, so that using equation 5.1, we obtain a discrete stochastic linear map for the proxy peak sequence $sn$,
$sn+1=Λ+(1-Λ)e-Δtnτsn,$
(5.2)
where the $Δtn$ stochastic sequence is predetermined from the point process model and the factor $(1-Λ)eΔtnτ$ is a stochastic coefficient for the discrete map, equation 5.3. This map is identical to the presynaptic adaptation $hn$-sequence map, equation 4.13. Using equation 4.15, we express the ISI in terms of $hn+1$ and $hn$, and it can easily verified that if $sn=hn$, then $sn+1=hn+1$. If $s0$ is initialized differently from $h0$, the discrete map equation 5.2, can be verified to satisfy certain regularity and stationarity properties that guarantee that $sn→hn$ almost surely as $n$ increases for any initialization (Brandt, 1986), as illustrated for a few example $sn$ sequences in Figure 3 (right).
Thus far, we have discussed only the example 3 adaptation dynamics. Examples 1 and 2 involve a purely additive activation that can be expressed as a limiting case of the example 3 model that is far from saturation. Let $c≫1$ be a large, rescaling parameter defining a change of variables $cS→S$ and $csn→sn$, and let $0<Λ≪1$ be very small so that $A=cΛ$ is fixed, and so $1-Λ∼1$. Then equation 5.3 becomes
$sn+1=A+e-Δtnτsn,$
(5.3)

The above arguments show that standard models of synaptic dynamics are equivalent to how we have defined presynaptic adaptation dynamics. Provided kinetics are matched between adaptation and the postsynaptic processes, any information in the adaptation states is potentially representable and transmissible postsynaptically, as depicted schematically in Figure 3.

6  Weak Signal Detection Performance of Spike Pattern and Rate Codes Using Fisher Information

In this section, we show that the spike pattern code exhibits the best possible detection performance for weak signal perturbations, exceeding rate coding for the adapting point process model.

Consider a $λ$-rate model that contains a multiplicative excitability parameter $α>0$. Example 3 possesses this parameter ($λ=αe-βH$), but any model can be modified to include the parameter (e.g., by writing $λα(H)=αλ(H))$. The parameter $α$ is similar to parameters used in other studies: $λ1$ in Kass and Ventura (2001) and the $q$ function in Berry and Meister (1998). The parameter $α$ sets the baseline probability of spiking and can be interpreted as an input rate parameter.

Note that it is sometimes more convenient to work with the inverse rate or timescale parameter $θ≡1α$. We will use these parameters interchangeably depending on which is more convenient. We define a weak signal by a perturbation $α+Δα$ or $θ+Δθ$, where $α$ and $θ$ are considered baseline levels.

Statistical discrimination between baseline and perturbed excitability parameters is performed by analyzing the Fisher information per spike for the spike pattern and rate models, based on a sample of $N$ spikes. We define the spike pattern code by the IID sequence ${hn}n=1N$ and its associated probability distribution $q(h)$, whereas spike rate coding is based on the distribution of the sample mean ISI: $TN=1N∑n=1NΔtn$. Note that $TN$ is in units of time, so it is not itself a rate but rather an inverse rate; however, we call $TN$ a rate statistic for convenience.

Note that the baseline parameter $α$ must not be excessively large in order to guarantee the $ε$-criterion for refractoriness that permits the quasi-IID property. Hence, there is a restricted range in which our analysis is valid that we term the physiological range that is given by an upper bound $αUB$. Conversely, too small an $α$-value, while it satisfies theorem 1, is also not very interesting because the resulting $h$-sequence is near constant so quasi-IID property is trivially satisfied while the ISI sequence is an effective renewal process and uncorrelated. Therefore we define a lower bound $αLB$ that excludes these trivial cases as well. The appendix rigorously defines $αUB$ and $αLB$ in terms of parameters in the example 3 model.

We compare the pattern code with two versions of a rate-based code. In each code, we assess discrimination performance on a per spike basis, which is reasonable because each code originates from the same model and so the firing rate is equal for each code. The first rate code, which we simply refer to as the rate code is based on the $TN$ statistic that takes into account the serial ISI correlations that yield short-long-short ISI sequences that result in a reduced variance of $TN$ relative to the case where the ISIs are a renewal process (Chacron et al., 2001; Nawrot et al., 2007). The second rate code assumes, wrongly, that the ISI sequence originates from a renewal process, and so does not benefit from the variance reduction.

6.1  The Pattern Code

Statistical estimation of $θ$ is straightforward using the $hn$-sequence. Let
$q(h|θ)=1θK'exp-1θK$
(6.1)
be the density for $h$ conditioned on $θ$, where $θ=〈K〉$, where the parentheses represent the expectation taken with respect to either $h$ or $η$. The unbiased estimator of the spike pattern code $θ^=1N∑n=1NK(hn)$ is statistically efficient, meaning it meets the Cramér-Rao bound: $Var(θ^)=Var(K)=Ipattern-1$. The term $Ipattern$ is the Fisher information, which is defined below.
Let $q(h|θ)$ and $q(h|θ+Δθ)$ be densities for $h$ conditioned on the baseline and perturbed inputs, respectively. The a posteriori information gained about the perturbation $Δθ$ from a single datum $h$ depends on the Kullback-Leibler (KL) divergence, or relative entropy, between the two densities: $D(qθ+Δθ||qθ)$ (Cover & Thomas, 1991). The Kullback-Leibler (KL) divergence is Taylor expanded about $Δθ=0$ to be
$D(qθ+Δθ||qθ)=12Δθθ2+…=Ipattern(θ)2Δθ2+…$
(6.2)
whose leading order coefficient, by definition, is half the Fisher information $Ipattern$. The Fisher information of the exponential family model, equation 6.1, is
$Ipattern(θ)≡1θ2,Ipattern(α)≡1α2,$
(6.3)
which holds for both the timescale parameter $θ$ and rate parameter $α$. Because each $h$-variable is associated with a single spike, this Fisher information is a per spike quantity.

6.2  The Rate Code

The Fisher information of the spike rate statistic $TN=1N∑i=1NΔti$ is computed by reexpressing the ISIs in terms of IID adaptation states $hn$. Recall that $Δtn=A(hn+1)-B(hn)$ from equation 2.3. The statistic $TN$ is expressed as
$TN=1N∑n=1NΔtn=1N[A(hN+1)-B(h1)]+1N∑n=2N[A(hn)-B(hn)].$
(6.4)
For large $N$, the single value $A(hN+1)$ in equation 6.4 is replaced by $A(h1)$ with only negligible change. Then, by combining all terms in the sum (see equation 6.4) that are dependent on the same $hn$-variable, we get
$TN({hn}n=1N)=1N∑n=1NZ(hn),$
(6.5)
where
$Z(hn)=Zn≡A(hn)-B(hn)$
(6.6)
is a random time variable with mean identical to the mean ISI. However, unlike ISIs, the $Zn$ quantities, equation 6.6, are IID because they are functions of single IID variables $hn$.

The distribution of $Z(h)$ can be derived through a change of variables from $h$, and so it has the same Fisher information as the pattern code. However, rate code inferences are not made from isolated $Zn$-variables or isolated $Δtn$, but on the large-$N$ sample mean $TN$ equation 6.5, of the $Zn$-variables. The distribution of $TN$ is asymptotically normal for large-$N$ samples by the central limit theorem.

The mean and variance of $Z$, given as a function $θ$ (but could also be given in terms of $α$), can be related to the mean and standard deviation of ISI statistics and allows a useful means of analysis of the spike train. We present the statistics of $Z$ as a critical lemma that enables us to analyze large-$N$ ISI averages (the inverse rate statistic):

Lemma 1
(Relation of $Z$ Statistics to ISI Statistics). The first- and second-order statistics of $Z$, equation 6.6, are expressible as
$μZ(θ)=μΔt(θ)=〈Δt〉=〈Z〉,$
(6.7)
$σZ2(θ)=Var(Z)=NVar(TN)=Var(Δt)+2(N-1)NCov(Δtn,Δtn+1)$
(6.8)
$≈Var(Δt)+2Cov(Δtn,Δtn+1),$
(6.9)
where for large $N$, the factor $(N-1)/N→1$ in equation 6.8 to produce equation 6.9.

Note that for large $N$, the mean ISI random variable, equation 6.5, approaches normality, where equations 6.7 and 6.9 define the mean and variance. Importantly, the variance of $Z$, equation 6.9, represents the variance reduction phenomenon of negative serially correlated ISI sequences.

Figure 4A shows, for the example 3 model, the ISI mean ($μZ$) and standard deviation ($σZ$) as a function of $α$, with $αLB$ and $αUB$ marked on the graph representing the physiological range. There are significant negative serial correlations $ρΔt(1)<0$ within the physiological range, concomitant to a coefficient of ISI variation $σZ/μZ$ being substantially less than unity (see Figure 4B) that defines non-Poisson and negative serially correlated spiking.

The rate code Fisher information per spike (i.e., divided by $N$) with mean and standard deviation that have $θ$ (or $α$) parameter dependencies is given by a standard identity for normally distributed statistics of the variable $TN$ by using the identities from lemma 4:
$Irate(θ)=1N〈[∂θln(φ(TN;μTN(θ),σTN2(θ))]2〉$
(6.10)
$=ddθμZ(θ)σZ(θ)2+12NddθσZ2(θ)σZ2(θ)2$
(6.11)
$≈ddθμZ(θ)σZ(θ)2,$
(6.12)
where $φ$ is the normal probability density and the second term in equation 6.11 is insignificant for large $N$ and so is dropped.
Figure 4:

Information loss in the rate code is due to insufficient mean ISI change relative to ISI standard deviation for the stochastic adapting model. Here, example 3 is used with parameters $β=40$, $τ=20$, and $Λ=0.2$. (A) Mean and standard deviation, along with Monte Carlo estimates, of the rate code variable $Z$ in equation 6.6. (B) ISI serial correlation with Monte Carlo estimate and coefficient of variation. (C) The three coding efficiencies. Panels A to C are plotted as a function of the excitability parameter $α$ in the spike emission rate function $λ$. (D) The trajectory of ($μZ$, $σZ$) pairs for the rate code (blue) compared to a theoretical curve (green) that would need to be satisfied by the rate code to achieve the higher statistical efficiency of the pattern code. This theoretical code possesses the same ISI standard deviation as the stochastic adaptive model, but has an increased change in mean ISI as a function of $α$. Insets show gaussian ISI densities for two baseline excitation values $α=10$ (dark blue and dark green squares) and $α=100$ (light blue and light green squares) for both the rate code and the theoretical curve.

Figure 4:

Information loss in the rate code is due to insufficient mean ISI change relative to ISI standard deviation for the stochastic adapting model. Here, example 3 is used with parameters $β=40$, $τ=20$, and $Λ=0.2$. (A) Mean and standard deviation, along with Monte Carlo estimates, of the rate code variable $Z$ in equation 6.6. (B) ISI serial correlation with Monte Carlo estimate and coefficient of variation. (C) The three coding efficiencies. Panels A to C are plotted as a function of the excitability parameter $α$ in the spike emission rate function $λ$. (D) The trajectory of ($μZ$, $σZ$) pairs for the rate code (blue) compared to a theoretical curve (green) that would need to be satisfied by the rate code to achieve the higher statistical efficiency of the pattern code. This theoretical code possesses the same ISI standard deviation as the stochastic adaptive model, but has an increased change in mean ISI as a function of $α$. Insets show gaussian ISI densities for two baseline excitation values $α=10$ (dark blue and dark green squares) and $α=100$ (light blue and light green squares) for both the rate code and the theoretical curve.

The ratio of the rate code Fisher information, equation 6.12 over pattern code Fisher information, equation 6.3, is termed the statistical efficiency:
$Eratevs.pattern≡Irate(θ)Ipattern(θ).$
(6.13)
The efficiency measures the fraction of information contained in the rate code based on the ${Tn}n=1N$ relative to the spike pattern code represented by the ${hn}n=1N$.

We note that the rate code Fisher information in equation 6.12 captures the variance-reduction effect of serially correlated ISI sequences. That is, $σZ2(θ)$ in equation 6.12 is smaller than what it would be if the ISI sequence were a renewal process, making the ratio in equation 6.12 larger than it would be otherwise (see below for a detailed analysis of this point). Despite this variance-reduction effect, the Fisher information, equation 6.12, is still bounded above by the pattern code Fisher information, equation 6.3, due to the gain change reflected in the numerator, as we rigorously establish in the following theorem.

Theorem 2.

The statistical efficiency $Eratevs.pattern$, equation 6.13, is bounded by unity, provided the adaptation process $hn$ is IID and $N$ is large.

A proof of theorem 5 is given in the appendix.

6.3  The Renewal Rate Code

The renewal rate code differs from the rate code in that it does not benefit from reduced variance, so it necessarily has a smaller Fisher information. Note that $σZ$ captures the reduced standard deviation of the serially correlated ISI sequence relative to the comparative renewal sequence. That is, $Var(Z)=NVar(TN)$. Note that by solving for $Var(Δt)$ in equation 6.8, we find
$Var(Δtrenewal)=Var(Z)-2Cov(Δtn,Δtn+1),$
(6.14)
where $Var(Δtrenewal)=Var(Δt)$—the “renewal” subscript is added for clarity, and $Cov(Δtn,Δtn+1)$ is the negative-valued serial ISI covariance that accounts for the increase in renewal-sequence ISI variance relative in equation 6.14 to the nonrenewal sequence—that is,
$Var(Z)≤Var(Δtrenewal),$
(6.15)
which expresses a similar finding to that in previous studies (Chacron et al., 2001, 2004; Farkhooi et al., 2011). The inequality 6.15 yields a corollary. Replacing $Var(Z)$ with $Var(Δtrenewal)$ in equation 6.12 implies that the Fisher information of the renewal spike train is bounded above by that for the nonrenewal rate code:
Corollary 3.
The Fisher information is ordered
$Irenewal≤Irate≤Ipattern;$
(6.16)
consequently, the efficiencies are ordered
$Erenewalvs.pattern≤Eratevs.pattern≤1,$
(6.17)
$Erenewalvs.rate≤1.$
(6.18)

We performed a numerical investigation to determine the magnitude of the aforementioned inequalities in theorem 5 and corollary 2. Figure 4C shows that the ordering of efficiencies over the physiological range is $Erenewalvs.pattern using representative parameters $β=40$, $τ=20$, and $Λ=0.2$. The rate code outperforms the renewal rate code, consistent with numerous prior studies (Ratnam & Nelson, 2000; Chacron et al., 2001, 2005; Goense & Ratnam, 2003; Nikitin et al., 2012), and the pattern code outperforms both of the rate codes. The pattern code performance premium above any rate code is roughly 20% to 30% over the physiological range. Note also that this 20% to 30% factor holds over a wide range of parameters $α$, $β$, $τ$, and $Λ$, as detailed in the appendix.

What accounts for the diminished Fisher information of the rate code? In terms of statistical estimation, the simple answer is that the pattern code, interpreted through the IID $hn$-sequence data and distribution $q$, can form a sufficient statistic for $θ$; the rate code utilizing $TN$ cannot. However, the diminished efficiency can also be understood in terms of physiology. The introduction of the variance-reducing effect of negative serial ISI correlations on the spike rate is a direct result of the adaptation $hn$ that defines the spike pattern code. It would seem that on first inspection, the coding performance of the rate code would be commensurate with the spike pattern code. However, the introduction of adaptation also reduces spike emission. So the adaptation that produces variance reduction in the denominator of equation 6.12 also reduces the firing rate change $dμZdθ$ in the numerator in equation 6.12, thereby offsetting at least some of the benefit of variance-reduction.

Figure 4D depicts this interaction between firing rate and firing rate variability with the curve $(μZ,σZ)$ parameterized by $α$ (the blue line with black squares marking the physiological range). The curve goes from upper right to lower left for increasing $α$-values. We used $β=40$, $τ=20$, and $Λ=0.2$ for this example. Any curve that curves further to the lower left corner, where a larger mean change relative to a fixed variance change, will necessarily yield a more discriminable rate code.

Consider example points $α=10$ and $α=100$ in the physiological range depicted in Figure 4D. The means of the ISI distributions obtained for these two excitation values differs by 5.25 msec (see the lower right inset in Figure 4D, blue lines). In order for a rate code to reach information parity with a spike pattern code, it would require a significantly greater mean ISI change than what was observed. The theoretical mean ISI change for Fisher information parity can be computed from the relation between a small input change $dα$ and the resulting change in mean $dμZ$. By equating equation 6.12 with $Ipattern$ (ensuring $E=1$) and solving for $dμZ$, we get
$dμZ(α)=Ipattern(α)σZ(α)dα=1ασZ(α)dα.$
(6.19)
Integrating equation 6.19 yields a maximally efficient theoretical mean and standard deviation curve $(μZ,σZ)$ shown in Figure 4D (green line). This curve represents a putative rate code that could achieve Fisher information parity with the pattern code. The theoretical mean ISI change was 6.02 msec over the example points (see the central inset) which is 14.67% larger than the mean ISI change achievable through the actual rate model. This establishes that while adaptation can reduce the variance in the firing rate statistic (see equation 6.15), it also reduces firing rate gain in response to input, and the resulting statistical discriminability is substantially less than what is available by the more statistically efficient pattern code.

7  Discussion

Point process models incorporating a spike-triggered adaptation produce the kind of serial ISI correlations present in many neuron types (Nawrot et al., 2007; Farkhooi et al., 2009) and can be modeled using point process models incorporating adaptation (Muller et al., 2007; Nesse et al., 2010; Farkhooi et al., 2011). In broad classes of these models, the IID sequence of random variables $hn$ representing adaptation levels serves as a decorrelating transformation of the correlated ISI sequence. This IID representation requires that adaptation be strong enough to produce a refractory period, a natural constraint for neuron models.

The quasi-IID property also requires that the spike emission function $λ(·)$ be a monotonically decreasing function of the postspike adaptation $H(t)$. Furthermore, the adaptation must monotonically increase following a spike. Such requirements on adaptation are entirely plausible, and postspike exponential decay of adaptation is utilized in some point process models (Keat et al., 2001). However, this is not strictly observed in practice. Studies by Kass and Ventura (2001) and Truccolo et al. (2005) fitted point process models to observe empirically estimated nonmonotonic spike probabilities after an initial spike. That is, their postspike probability functions increased from zero to reach a peak but then declined modestly to a steady state probability thereafter. However, in both cases, the typical ISI fell mostly within the monotonic zone. Similarly, using a linear response model approximating the Hodgkin-Huxley model, fitted kernels exhibit near-zero influence of inputs postspike followed by a postspike rebound peak (Kistler, Gerstner, & Hemmen, 1997).

Numerous studies have established that serial ISI correlations reduce spike count variance, thereby making a given firing rate change more discriminable (Chacron et al., 2001, 2004, 2005; Lüdtke & Nelson, 2006; Nawrot et al., 2007; Maimon & Assad, 2009; Farkhooi et al., 2009, 2011; Avila-Akerberg & Chacron, 2011; Nikitin et al., 2012). This long-standing finding represented a performance benchmark for negative serially correlated spike trains. Our results here establish that there is a level of coding performance beyond this previous benchmark. Firing rate coding, whether from correlated ISIs or renewal ISIs, is less informative than spike pattern coding as represented by adaptation states $hn$ for the detection of weak input changes. This is because the presence of adaptation, while it can reduce spike count variance, also reduces firing rate gain, thereby offsetting some performance. The stronger performance of the spike pattern code relative to rate codes arises because the firing rate is not a sufficient statistic for adapting cell excitability, while the pattern code based on the $hn$ adaptation states is equipped with a sufficient statistic.

Our work has focused exclusively on the negative ISI correlations because it has been deeply studied both experimentally and theoretically. However, it is known that ISI correlations come in various forms, such as a mixture of positive and negative correlations at different lags, and even oscillatory ISI correlations. The precise makeup of correlations determines spike train statistics such as the Fano factor over asymptotically large counting windows. Recent theoretical work (Schwalger & Lindner, 2013; Shiau, Schwalger, & Lindner, 2015) and modeling studies (Fuwape & Neiman, 2008) have explored the dynamical origin of these correlations and their impact on information transfer, concluding that the nonrenewal property of the ISI sequence can often enhance information transfer. It remains to be seen how our results on the quasi-IID properties of large classes of adaptation processes can provide further insight into the single cell and cellular network coding abilities of those dynamical processes. This poses a serious challenge given the more complex nature of the putative dynamics responsible for the more complex ISI correlations.

It is noteworthy that the inverse-square relationship $1/α2$ of the spike pattern Fisher information is identical to a nonadapting, nonrefractory, renewal Poisson point process $λ=α$, owing to the fact that both are in the exponential family of distributions. With Poisson spiking, information available to discriminate inputs is reflected purely by varying spike rates with no upper bound on the rate or spike count and no refractory periods that characterize real neurons. In the adapting model, this same information is reflected in the adaptation states and requires smaller spike rate variation. Neuron dynamics permit only a limited dynamic range of spike rate variation. The adaptation-based scheme therefore permits maximum coding performance by neurons operating under physiological constraints on their dynamic range.

We can also speculate how such neural spike pattern representations would need to be decoded. Decoding the firing rate necessarily involves some physiological proxy of a time-windowed spike count. In contrast, the decoding of spike patterns is achieved by assessing the amplitude of the IID postspike adaptation amplitudes $hn$. Because the kinetics of spike-triggered adaptation are much like that of synaptic neurotransmission, there could be a corresponding quasi-IID postsynaptic proxy of each $hn$ amplitude, provided there are matched kinetics (see Figure 3). A similar approach of representing ISI patterns by synaptic dynamics with appropriately tuned kinetics to infer presynaptic states has been proposed (Pfister et al., 2010; Lüdtke & Nelson, 2006). Furthermore, a recent experimental study has established that postsynaptic potential peaks, akin to the $hn$-sequence of adaptation states, are decorrelated when the driving presynaptic ISI structure exhibits negative serial correlations, but only if the synaptic kinetics of the postsynaptic potentials, consisting of multiple neurotransmitter systems, were not pharmacologically perturbed (Marcoux et al., 2016). This result suggests that theoretical findings in this letter may be applied to the decoding of neuron data. Moreover, such pattern-based discrimination may be necessary to account for performance in neural systems where rate coding involving spike patterns with negative serial ISI correlations is demonstrably insufficient to explain behavioral performance, especially near the threshold of sensory detection (Knudsen, 1974; Gussin et al., 2007; Jacobs et al., 2009; Jung et al., 2016; Fettiplace & Hackney, 2006; Heil & Peterson, 2015).

Finally, an interesting avenue to pursue in future work is to understand more fully the encoding of temporal signals in the presence of the quasi-IID adaptation states and the negatively correlated ISIs. This would go beyond the simple stimuli in the form of step increases in amplitude that were considered here. Previous work (Ratnam & Nelson, 2000; Chacron et al., 2001, 2004) comparing models with and without these negative ISI correlations reveals a noise-shaping effect. This causes a decrease in the baseline of the power spectral density of the spontaneous spike trains from the model at lower frequencies. This leads to a signal detection and encoding enhancement, since the signal-to-noise ratio increases in that range. The effect is related to the regularization of spike counts by the anticorrelation discussed here (Chacron et al., 2001; Nesse et al., 2010). However, the decrease of the noise floor at lower frequencies is accompanied by an increase at the higher frequencies. Consequently, the coding benefits gained at low frequencies lead to losses at higher frequencies, which effectively limits the bandwidth of the system to the lower range. The precise manner in which the theorems described here for the quasi-IID behavior translate into signal encoding benefits by the adaptation states and the ISI sequence remains to be explored, a task left for future work.

Appendix

A.1  Proof of Theorem 1, the Quasi-IID Adaptation Property

The value $η=η0$ is the starting value after activation from some previous spike. We can mandate that $η0≤ηε$ because the time activation function $f$ never exceeds $ηε$.

Consider the first case in equation 3.1: $η≤ηε$ with survival function $P(η|η0)=e-K(η|η0)$. Note that $K(η|η0)<ε$, so
$|1-e-K(η|η0)|=|1-(1-K+…)|≤K<ε.$
(A.1)
The inequality is true by the alternating series estimation theorem, owing to the fact that $e-K$ can be represented as an alternating series expanded about $K=0$. The theorem states in our case here that the difference between the full series and the zeroth-order approximation of $e-K$ (which is 1) is bounded by the first-order term, which is $K$.

This expansion establishes the first case, equation 3.1. The statement means that for the initial interval of time $[η0,ηε)$, less than an $ε$-fraction of the ensemble spikes, so the survival function is $ε$-close to unity.

Now consider the second case in equation 3.1: $η>ηε$. The potential function is then
$K(η|η0)=∫η0ηελ(H(ξ))dξ+∫ηεηλ(H(ξ))dξ=K(ηε|η0)+K(η|ηε).$
(A.2)
Let $P(η)=e-K(η|ηε)$, so
$P(η|η0)=e-K(η|ηε)e-K(ηε|η0)=P(η)(1+[e-K(ηε|η0)-1]).$
(A.3)
Therefore,
$|P(η)-P(η|η0)|=|P(η)-P(η)(1+[e-K(ηε|η0)-1])|$
(A.4)
$=|P(η)||e-K(ηε|η0)-1|<ε$
(A.5)
because $|e-K(ηε|η0)-1|<ε$ and $0≤P(η)≤1$. $□$

A.2  Proof of Corollary 3

The covariance of subsequent ISIs is computed using the representation $Δtn=A(hn+1)-B(hn)$ as follows:
$Cov(Δtn,Δtn+1)=〈[A(hn+1)-B(hn)][A(hn+2)-B(hn+1)]〉-〈Δt〉2$
(A.6)
$=〈A〉〈B〉-〈AB〉,$
(A.7)
where the brackets represent expectation with respect to the distribution $q(h)$. The correlation is then $Corr(Δtn,Δtn+1)=Cov(Δtn,Δtn+1)/Var(Δt)$.
We establish that $Cov(Δtn,Δtn+1)<0$ first. Note that $A(h)$ and $B(h)$ are both monotonic decreasing because $H-1$ is monotonic decreasing and $f-1$ is monotonic increasing. This satisfies the hypotheses of Chebyshev's integral inequality:
$〈A〉〈B〉<〈AB〉.$
(A.8)
Therefore, $Cov(Δtn,Δtn+1)$ is negative (see equation A.7).

A serial correlation no smaller than $-12$ is an established fact of event sequences (Cox & Lewis, 1966; Chacron et al., 2004; Nawrot et al., 2007; Farkhooi et al., 2009), which concludes the proof. $□$

A.3  Proof of Theorem 5

The numerator of equation 6.12 involves the squared derivative of the mean $μZ$ with respect to $θ$. Because $Z$ is a function of $h$, we can compute the mean using $q(h)$. Note that $q(h)=K'(h)θe-K(h)/θ$, so
$μZ=〈Z〉=∫Z(h)K'(h)θe-K(h)θdh.$
(A.9)
However, in order to make progress with computing $dμZdθ$, the first two moments of $K(h)$ must be computed. The mean, or first moment, of $K$ is computed to be
$〈K〉=∫KK'θe-K/θdh$
(A.10)
$=θ∫Kθe-K/θdh=θ.$
(A.11)
The second moment is
$〈K2〉=∫K2K'θe-K/θdh$
(A.12)
$=[-K2e-K/θ|K=0K=∞+2θ∫KK'θe-K/θdh=2θ2.$
(A.13)
So the variance is then computed to meet the Cramér-Rao bound $Var(K)=2θ2-θ2=θ2=Ipattern-1$.
The squared-derivative $dμZdθ2$ in equation 6.12, which is the square of the gain of the Z-rate code, can be expressed as
$dμZdθ=ddθ∫ZK'θe-K/θdh$
(A.14)
$=-1θ∫ZK'θe-K/θdh+1θ2∫ZKK'θe-K/θdh$
(A.15)
$=-1θ〈Z〉+1θ2〈ZK〉$
(A.16)
$-1θ2〈K〉〈Z〉+1θ2〈ZK〉$
(A.17)
$=-〈K〉〈Z〉+〈ZK〉Var(K)$
(A.18)
$=Cov(K,Z)Var(K).$
(A.19)
Squaring the derivative, we get
$dμZdθ2=Cov(K,Z)2Var(K)2$
(A.20)
$≤Var(K)Var(Z)Var(K)2=Var(Z)Var(K)$
(A.21)
$=Var(Z)1θ2=σZ2Ipattern(θ).$
(A.22)
The inequality in equation A.21 is the classical covariance inequality. Also note $Var(K)=θ2=Ipattern(θ)-1$ in equation A.22 by the Cramér-Rao identity for efficient estimators. By dividing equations A.20 to A.22 by $σZ2$, the expression becomes equation 6.12, and the inequality states $Irate≤Ipattern$ or, equivalently $Eratevs.pattern≡Irate(θ)/Ipattern(θ)≤1$, which proves theorem 5. $□$

A.4  The Effect of Auxiliary Parameters on Spike Statistics and Fisher Information

The coding efficiency of the rate code was found to be 20% to 30% less than that of the pattern code for the more biophysically realistic example 3 model for a specific set of parameters $β=40$, $τ=20$, and $Λ=0.2$, with $α$ varying over what we termed the physiological range. Here we rigorously define this range and perform an extensive parameter scan to understand how all parameters $β$, $Λ$, and $τ$ interact with $α$ in the example 3 model to affect coding efficiency. We establish that the 20% to 30% efficiency loss we reported for the rate code is representative over an exhaustive scan of the parameters.

Physiologically realistic spiking is characterized by (1) exhibiting significant serial ISI correlations between subsequent ISIs, with $ρΔt(1)$ being sufficiently negative; (2) excitability is not so high that the quasi-IID property holds—a postspiking refractory period is present; and (3) the adaptation pattern sequence $hn$ exhibits nontrivial variability levels so that the ISI sequence and the $hn$ sequence are not a constant or near-constant sequence. Conditions 1 to 3 imply that the ISI coefficient of variation (CV) is sufficiently below unity, unlike a Poisson process, yet also not too close to zero like that of a constant ISI sequence. Over a very large range of parameter combinations, we show that there exist $α$ intervals where the spiking behavior is consistent with conditions 1 to 3.

Within this physiological range satisfied by conditions 1 to 3, we find the coding efficiency $Eratevs.pattern$ is consistently in the 70% to 80% range. In other more extreme parameter ranges, coding efficiencies were much lower and higher, but the associated spiking statistics did not meet criteria 1 to 3. This is explored rigorously as follows.

The adapting point process model comprises four parameters: $τ$, $β$, $Λ$, and $α$. The parameter $τ$ represents the characteristic timescale of the system and can be set to unity without loss of generality; however, in what follows, we include $τ$ for clarity. The parameter $Λ$ is the fraction of ion channels that transition to the open state (see the example 3 derivation), so by definition, $Λ∈(0,1)$. The remaining three parameters must be considered in relation to conditions 1 to 3. Condition 1 is in part satisfied when $α$ is sufficiently large to produce robust negative serial ISI correlations. This occurs when $α$-values are larger than some scaling factor $c$ of the inverse timescale $τ-1$:
$cτ-1<α.$
(A.23)
This ensures that spikes will be in short enough succession that the adaptation process $H(t)$ can visibly summate rather than decay to near-zero between spikes.
Requirement 2 guarantees the quasi-IID condition and is contingent on $β$ and $Λ$ being sufficiently large. This condition is guaranteed when spiking probability immediately after a spike is very small, as guaranteed by the inequality
$ατe-βΛ(1-Λ)Λ<ε,$
(A.24)
which was derived in equation 4.10. Equation A.24 yields an upper bound on $ατ$ for fixed $β$ and $Λ$. Combining the above two inequalities, equations A.23 and A.24, we get
$cτ-1<α<ετ-1eβΛΛ(1-Λ).$
(A.25)
However, for our purposes here, it is best to express the $α$-range on a logarithmic scale: $α=eβs$, where $s$ is a change of variable of the excitability parameter, resulting in
$ln(cτ-1)β
(A.26)
where, without loss of generality, $c=4$ and $ln(ε)=-3$ ($ε∼0.0498$) were chosen to form a specific $s$-range for every $β$ and $Λ$ (where the $s$-range exists). Figure 5 shows the resulting “islands” of $s$-$β$ parameter space that conform to the bound given by equation A.26. Each column shows the islands for a range of $Λ$-values. The island for the smallest $Λ=0.1$-value was very small, so it is enlarged in an inset. Each row shows a color-coded map of various spike statistics and coding efficiencies (see below).
Figure 5:

Effect of the other adapting point process model parameters on spike statistics and coding efficiency, depicted for increasing $Λ$-values (indicated at the top of each column) from left to right. Each panel scans a range of $s$-excitability and $β$-adaptation parameters. Top row: Serial correlation coefficient at lag 1, that is, $ρΔt(1)$. Second row: The coefficient of variation CV. Third and fourth rows: Coding efficiency $E$ comparing the pattern code to the renewal rate code, and pattern code to the rate code, respectively. Each plot depicts one of the quantities noted over a physiological range of $s$-$β$ values given by the bounds in equation A.26. The dashed black and white lines in the bottom two rows depict the parameters yielding realistic coding efficiencies; they are subsets of the areas framed by the thick black lines obtained in the top two rows that correspond to realistic baseline spike statistics.

Figure 5:

Effect of the other adapting point process model parameters on spike statistics and coding efficiency, depicted for increasing $Λ$-values (indicated at the top of each column) from left to right. Each panel scans a range of $s$-excitability and $β$-adaptation parameters. Top row: Serial correlation coefficient at lag 1, that is, $ρΔt(1)$. Second row: The coefficient of variation CV. Third and fourth rows: Coding efficiency $E$ comparing the pattern code to the renewal rate code, and pattern code to the rate code, respectively. Each plot depicts one of the quantities noted over a physiological range of $s$-$β$ values given by the bounds in equation A.26. The dashed black and white lines in the bottom two rows depict the parameters yielding realistic coding efficiencies; they are subsets of the areas framed by the thick black lines obtained in the top two rows that correspond to realistic baseline spike statistics.

The choice of $c=4$ in equation A.26 was not very restrictive. Increasing $c$ further would result in a narrower admissible parameter range; however, a broader range was chosen so that the behavior of the adapting point process model is better illustrated. On the left-hand sides of all the parameter islands, the spike patterns possess weak, negative serial ISI correlations (see the darker red region of the color map, top row of Figure 5), but increasing $s$ yields stronger negative correlations. We will only consider realistic spiking that exhibits serial correlations below $-0.2$ that is typical of reported serial correlations (Farkhooi et al., 2009; Nawrot et al., 2007), a region to the right of the thick black line in the top row of Figure 5 that is considered physiological. The right-hand boundary of the parameter island defines the upper bound of equation A.26 that satisfies condition 2.

The condition on the CV is shown in the second row of Figure 5. The admissible CV range was defined to be 0.2 to 0.3, an intermediate level of spiking variability in this model. Studies of real neurons typically report larger CV-values (Nawrot et al., 2007; Farkhooi et al., 2009), but these must be estimated from long sampling windows, typically much longer than used for sensory detection, which could include intermediate-to-longer-timescale fluctuations that are not present in a stationary model. These fluctuations will necessarily elevate experimentally computed estimates of variability (Perkel et al., 1967). This 0.2 to 0.3 CV range was framed by the lower and upper black lines in the second row of Figure 5. Increasing $β$- and $s$-levels for any $Λ$ had the effect of increasing the regularity of the ISI sequence, resulting in a lower CV.

Together, conditions 1 to 3 formed an admissible parameter range for realistic serially correlated spiking (Nawrot et al., 2007; Farkhooi et al., 2009). These combinations occurred to the right side of the ISI correlation boundary and to the left of the edge of the parameter island formed by equation A.26, and between the upper and lower CV boundaries; this defines a region of parameter space for particular levels of $Λ$-values. For some of the larger $Λ$-values, there was no admissible parameter region because the increased $Λ$-value “pushed” the physiological parameter range off the island. A very large $Λ$-value in the range $0.6<Λ<1$ results in reduced variability for both the ISIs and the $hn$-sequence because the adaptation levels become effectively saturated (not shown).

Two coding efficiencies were assessed and plotted in the third and fourth rows of Figure 5 within the admissible parameter range island. The border enclosing the physiological range (satisfying properties 1 to 3) was drawn as the white-black dashed lines. The first efficiency (third row) compared the adaptation pattern code with the renewal rate code. The renewal rate code exhibited a consistently low efficiency across a range of $Λ$-values, with an upper limit around 0.45 to 0.5, and a lower limit around 0.23 to 0.39 over the $Λ$-columns.

The second efficiency compared the pattern code to the spike rate code for the adapting model, with results shown in the fourth row of Figure 5. The negative ISI correlations reduce the variance of the ISI sequence, and thus of spike counts, and improve detection performance relative to renewal ISIs. The upper coding efficiency level was observed to be approximately 0.75 to 0.85, while the lower coding efficiency level was approximately 0.71 to 0.75 over the $Λ$-columns. We conclude that for realistic spike trains exhibiting negative serial correlations and IID adaptation, the coding efficiency over this exhaustive parameter scan for the example 3 model is roughly 70% to 80%; in other words, the rate code detection performance is 20% to 30% less than that of the pattern code.

Acknowledgments

We thank Gregory Rice (University of Waterloo), and Lajos Horvath (University of Utah) for their helpful discussions. A.L. acknowledges support from the NSERC Canada DG program, and L.M. acknowledges support from CIHR 153143.

References

Avila-Akerberg
,
O.
, &
Chacron
,
M. J.
(
2011
).
Nonrenewal spike train statistics: Causes and functional consequences on neural coding
.
Exp. Brain Res.
,
210
,
353
371
.
Berry
,
M.
, &
Meister
,
M.
(
1998
).
Refractoriness and neural precision
.
Journal of Neuroscience
,
18
,
2200
2211
.
Brandt
,
A.
(
1986
).
The stochastic equation $yn+1=anyn+bn$ with stationary coefficients
.
,
18
,
211
220
.
Chacron
,
M. J.
,
Lindner
,
B.
, &
Longtin
,
A.
(
2004
).
Noise shaping by interval correlations increases information transfer
.
Phys. Rev. Lett.
,
92
, 080601.
Chacron
,
M. J.
,
Longtin
,
A.
, &
Maler
,
L.
(
2001
).
Negative interspike interval correlations increase the neuronal capacity for encoding time-dependent stimuli
.
Journal of Neuroscience
,
21
(
14
),
5328
5343
.
Chacron
,
M. J.
,
Maler
,
L.
, &
Bastian
,
J.
(
2005
).
Electroreceptor neuron dynamics shape information transmission
.
Nat. Neurosci.
,
21
,
673
678
.
Citi
,
L.
,
Ba
,
D.
,
Brown
,
E. N.
, &
Barbieri
,
R.
(
2014
).
Likelihood methods for point processes with refractoriness
.
Neural Computation
,
26
(
2
),
237
263
.
Cover
,
T.
, &
Thomas
,
J.
(
1991
).
Elements of information theory
.
New York
:
Wiley
.
Cox
,
D. R.
, &
Lewis
,
P. A. W.
(
1966
).
The statistical analysis of series of events
.
New York
:
Wiley
.
Farkhooi
,
F.
,
Muller
,
E.
, &
Nawrot
,
M. P.
(
2011
).
Adaptation reduces variability of the neuronal population code
.
Phys. Rev. E
,
83
, 050905.
Farkhooi
,
F.
,
Strube-Bloss
,
M. F.
, &
Nawrot
,
M. P.
(
2009
).
Serial correlation in neural spike trains: Experimental evidence, stochastic modelling, and single neuron variability
.
Phys. Rev. E
,
79
, 021905.
Fettiplace
,
R.
, &
Hackney
,
C. M.
(
2006
).
The sensory and motor roles of auditory hair cells
.
Nat. Rev. Neurosci.
,
7
,
19
29
.
Fuwape
,
I.
, &
Neiman
,
A. B.
(
2008
).
Spontaneous firing statistics and information transfer in electroreceptors of paddlefish
.
Phys. Rev. E
,
78
, 051922.
Goense
,
J. B. M.
, &
Ratnam
,
R.
(
2003
).
Continuous detection of weak sensory signals in afferent spike trains: The role of anti-correlated interspike intervals in detection performance
.
J. Comp. Physiol. A
,
189
,
741
759
.
Gussin
,
G.
,
Benda
,
J.
, &
Maler
,
L.
(
2007
).
Limits of linear rate coding of dynamic stimuli by electroreceptor afferents
.
J. Neurophysiol.
,
97
,
2917
2929
.
Heil
,
P.
, &
Peterson
,
A. J.
(
2015
).
Basic response properties of auditory nerve fibers: A review
.
Cell Tissue Res.
,
361
,
129
158
.
Jacobs
,
A. L.
,
Fridman
,
G.
,
Douglas
,
R. M.
,
Alam
,
N. M.
,
Latham
,
P.
,
Prusky
,
G. T.
, &
Nirenberg
,
S.
(
2009
).
Ruling out and ruling in neural codes
.
,
106
,
5936
5941
.
Jung
,
S. N.
,
Longtin
,
A.
, &
Maler
,
L.
(
2016
).
Weak signal amplification and detection by higher-order sensory neurons
.
J. Neurophysiol.
,
115
,
2158
2175
.
Kass
,
R. E.
, &
Ventura
,
V.
(
2001
).
A spike-train probability model
.
Neural Computation
,
13
(
8
),
1713
1720
.
Keat
,
J.
,
Reinagel
,
P.
,
Reid
,
R.
, &
Meister
,
M.
(
2001
).
Predicting every spike: A model for the responses of visual neurons
.
Neuron
,
30
(
3
),
803
817
.
Khanbabaie
,
R.
,
Nesse
,
W. H.
,
Longtin
,
A.
, &
Maler
,
L.
(
2010
).
Kinetics of fast short-term depression are matched to spike train statistics to reduce noise
.
J. Neurophysiol.
,
103
,
3337
3348
.
Kistler
,
W. M.
,
Gerstner
,
W.
, &
Hemmen
,
J. L. v.
(
1997
).
Reduction of the Hodgkin-Huxley equations to a single-variable threshold model
.
Neural Computation
,
9
(
5
),
1015
1045
.
Knudsen
,
E. I.
(
1974
).
Behavioral thresholds to electric signals in high frequency electric fish
.
J. Comp. Physiol.
,
91
,
333
353
.
London
,
M.
,
Roth
,
A.
,
Beeren
,
L.
,
Haüsser
,
M.
, &
Latham
,
P. E.
(
2010
).
Sensitivity to perturbations in vivo implies high noise and suggests rate coding in cortex
.
Nature
,
466
,
123
127
.
Lüdtke
,
N.
, &
Nelson
,
M. E.
(
2006
).
Short-term synaptic plasticity can enhance weak signal detectability in nonrenewal spike trains
.
Neural Computation
,
18
(
12
),
2879
2916
.
Maimon
,
G.
, &
,
J. A.
(
2009
).
Beyond Poisson: Increased spike-time regularity across primate parietal cortex
.
Neuron
,
62
,
426
440
.
Marcoux
,
C. M.
,
Clarke
,
S. E.
,
Nesse
,
W. H.
,
Longtin
,
A.
, &
Maler
,
L.
(
2016
).
Balanced ionotropic receptor dynamics support signal estimation via voltage-dependent membrane noise
.
J. Neurophysiol.
,
115
,
530
545
.
Muller
,
E.
,
Buesing
,
L.
,
Schemmel
,
J.
, &
Meier
,
K.
(
2007
).
.
Neural Computation
,
19
(
11
),
2958
3010
.
Nawrot
,
M. P.
,
Boucsein
,
C.
,
Rodriguez-Molina
,
V.
,
Aertsen
,
A.
,
Grün
,
S.
, &
Rotter
,
S.
(
2007
).
Serial interval statistics of spontaneous activity in cortical neurons in vivo and in vitro
.
Neurocomputing
,
70
(
10
),
1717
1722
.
Nesse
,
W. H.
,
Maler
,
L.
, &
Longtin
,
A.
(
2010
).
Biophysical information representation in temporally correlated spike trains
.
,
107
,
21973
21978
.
Nikitin
,
P.
,
Stocks
,
N.
, &
Bulsara
,
A.
(
2012
).
Enhancing the resolution of a sensor via negative correlation: A biologically inspired approach
.
Phys. Rev. Lett.
,
109
, 238103.
Paninski
,
L.
(
2004
).
Maximum likelihood estimation of cascade point-process neural encoding models
.
Network: Computation in Neural Systems
,
15
(
4
),
243
262
.
Panzeri
,
S.
,
Harvey
,
C.
,
Piasini
,
E.
,
Latham
,
P.
, &
Fellin
,
T.
(
2017
).
Cracking the neural code for sensory perception by combining statistics, intervention, and behaviour
.
Neuron
,
93
,
491
507
.
Perkel
,
D. H.
,
Gerstein
,
G. L.
, &
Moore
,
G. P.
(
1967
).
Neuronal spike trains and stochastic point processes. I. The single spike train
.
Biophysical Journal
,
7
(
4
),
391
418
.
Pfister
,
J. P.
,
Dayan
,
P.
, &
Lengyel
,
M.
(
2010
).
Synapses with short-term plasticity are optimal estimators of presynaptic membrane potentials
.
Nat. Neurosci.
,
13
,
1271
1275
.
Ratnam
,
R.
, &
Nelson
,
M. E.
(
2000
).
Nonrenewal statistics of electrosensory afferent spike trains: Implications for the detection of weak sensory signals
.
J. Neurosci.
,
20
,
6672
6683
.
Schwalger
,
T.
, &
Lindner
,
B.
(
2013
).
Patterns of interval correlations in neural oscillators with adaptation
.
Frontiers in Computational Neuroscience
,
7
, 164.
Shiau
,
L.
,
Schwalger
,
T.
, &
Lindner
,
B.
(
2015
).
Interspike interval correlation in a stochastic exponential integrate-and-fire model with subthreshold and spike-triggered adaptation
.
Journal of Computational Neuroscience
,
38
(
3
),
589
600
.
Stevens
,
C. F.
, &
,
A.
(
1995
).
Neural coding: The enigma of the brain
.
Curr. Biol.
,
5
,
1370
1371
.
Truccolo
,
W.
,
Eden
,
U. T.
,
Fellows
,
M. R.
,
Donoghue
,
J. P.
, &
Brown
,
E. N.
(
2005
).
A point process framework for relating neural spiking activity to spiking history, neural ensemble, and extrinsic covariate effects
.
Journal of Neurophysiology
,
93
(
2
),
1074
1089
.