History dependence disentangles the effects of input activation, reactivation and temporal depth of a binary autoregressive process.

To show the behavior of the measures in a well controlled setup, we analyzed a simple binary autoregressive process with varying temporal depth l (Fig 3A). The process evolves in discrete time steps, and has an active (1) or inactive (0) state. Active states are evoked either by external input with probability h, or by internal reactivations that are triggered by activity within the past l steps. Each past activation increases the reactivation probability by m, which regulates the strength of history dependence in the process. In the following, we describe how the measures behave as we vary each of the different model parameters, and then summarize the key difference between the measures.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. History dependence disentangles the effects of input activation, reactivation and temporal depth of a binary autoregressive process.

(A) In the binary autoregressive process, the state of the next time step (grey box) is active (one) either because of an input activation with probability h, or because of an internal reactivation. The internal activation is triggered by activity in the past l time steps (green), where each active state increases the activation probability by m. (B) Increasing the input activation probability h increases the total mutual information I_tot, although input activations are random and therefore not predictable. Normalizing the total mutual information by the entropy yields the total history dependence R_tot, which decreases mildly with h. (C) Autocorrelation C(T), lagged mutual information L(T) and gain in history dependence ΔR(T) decay differently with the time lag T. For l = 1 and m = 0.8 (top), autocorrelation C(T) decays exponentially with autocorrelation time τ_C, whereas L(T) decays faster due to the non-linearity of the mutual information. For l = 5 (bottom), C(T) and L(T) plateau over the temporal depth, and then decay much slower than for l = 1. In contrast, ΔR(T) is non-zero only for T shorter or equal to the temporal depth of the process, with much shorter timescale τ_R. Parameters m and h were adapted to match the firing rate and total history dependence between l = 1 and l = 5. (D) When increasing the reactivation probability m for l = 1, timescales of time-lagged measures τ_C and τ_L increase. For history dependence, the information timescale τ_R remains constant, but the total history dependence R_tot increases. (E) When varying the temporal depth l, all timescales increased. Parameters h and m were adapted to hold the firing rate and R_tot constant.

https://doi.org/10.1371/journal.pcbi.1008927.g003

The input strength h increases the firing rate and thus the spiking entropy H(spiking). This leads to a strong increase in the total mutual information , whereas the total history dependence R_tot is normalized by the entropy and does slightly decrease (Fig 3B). This slight decrease is expected from a sensible measure of history dependence, because the input is random and has no temporal dependence. In addition, input activations may fall together with internal activations, which slightly reduces the total history dependence.

In contrast, the total history dependence R_tot increases with the reactivation probability m, as expected (Fig 3D). For the autocorrelation, the reactivation probability m not only influences the magnitude of the correlation coefficients, but also the decay of the coefficients. For autoregressive processes (and l = 1), autocorrelation coefficients C(T) decay exponentially [14] (Fig 3C), where the autocorrelation time τ_C = −Δt/log(m) increases with m and diverges as m → 1 (Fig 3D). The lagged mutual information L(T) is a non-linear measure of time-lagged dependence, and has a very similar behavior as the autocorrelation, with a slightly faster decay and thus smaller generalized timescale τ_L (Fig 3C and 3D). Note that we normalized L(T) by the spiking entropy H to make it directly comparable to ΔR(T). In contrast to the time-lagged measures, the gain in history dependence ΔR(T) is only non-zero for T smaller or equal to the true temporal depth l of the process (Fig 3C). As a consequence, the information timescale τ_R does not increase with m for fixed l (Fig 3D).

Finally, the temporal depth l controls how far into the past activations depend on their preceding activity. Indeed, we find that the information timescale τ_R increases with l as expected (Fig 3C and 3E). Similarly, the timescales of the time-lagged measures τ_C and τ_L increase with the temporal depth l. Note that parameters m and h were adapted for each l to keep the firing rate and total history dependence R_tot constant, hence differences in the timescale can be unambiguously attributed to the increase in l.

To conclude, history dependence disentangles the effects of input activation, reactivation and temporal depth, which provides a comprehensive characterization of past dependencies in the autoregressive model. This is different from the total mutual information, which lacks the entropy normalization and is sensitive to the firing rate. This is also different from time-lagged measures, whose timescales are sensitive to both, the reactivation probability m and the temporal depth l. The confusion of effects in the timescales is rooted in the time-lagged nature of the measures—by quantifying past dependencies out of context, C(T) and L(T) also capture indirect, redundant dependencies onto past events. Indirect, redundant dependencies arise from unique dependencies, because past states that are uniquely predictive of future activities were in turn uniquely dependent on their own past. The stronger the unique dependence, the longer the indirect dependencies reach into the past, which increases the timescale of time-lagged measures. In contrast, indirect dependencies do not contribute to the history dependence, because they add no predictive information once more-recent past is taken into account.

History dependence dismisses redundant past dependencies and captures synergistic effects.

A key property of history dependence is that it evaluates past dependencies in the light of more-recent past. This allows the measure to dismiss indirect, redundant past dependencies and to capture synergistic effects. In three common models of neural spiking activity, we demonstrate how this leads to a substantially different characterization of past dependencies compared to time-lagged measures of temporal dependence.

First, we simulated a subsampled branching process [14], which is a minimal model for activity propagation in neural networks and captures key properties of spiking dynamics in cortex [15]. Similar to the binary autoregressive process, active neurons activate neurons in the next time step with probability m, the so-called branching parameter, and are activated externally with some probability h. The process was simulated in time steps of Δt = 4 ms with a population activity of 500 Hz, which was subsampled to obtain a single spike train with a firing rate of 5 Hz (Fig 4A). Similar to the binary autoregressive process, the autocorrelation decays exponentially with autocorrelation time τ_C = −Δt/log(m) = 198 ms, and the lagged mutual information decays slightly faster (Fig 4B). In comparison, the gain in history dependence ΔR decays much faster. When increasing the branching parameter m (for fixed firing rate), the total history dependence increased, as in the autoregressive process (S11 Fig). Strikingly, the timescale τ_R remained constant or even decreased for larger m > 0.967 and thus higher autocorrelation time τ_C > 120 ms (S11 Fig), which is different from the binary autoregressive process. The reason is that the branching process evolves at the population level, whereas history dependence is quantified at the single neuron level. Thereby, history dependence also captures indirect dependencies, because the own spiking history reflects the population activity. The higher the branching parameter m, the more informative past spikes are about the population activity, and the shorter is the timescale τ_R over which all the relevant information about the population activity can be collected. Thus, for the branching process, the total history dependence R_tot captures the influence of the branching parameter, whereas the information timescale τ_R behaves very differently from the timescales of time-lagged measures.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. History dependence dismisses redundant past dependencies and captures synergistic effects.

(A,B) Analysis of a subsampled branching process. (A) The population activity was simulated as a branching process (m = 0.98) and subsampled to yield the spike train of a single neuron (Materials and methods). (B) Autocorrelation C(T) and lagged mutual information L(T) include redundant dependencies and decay much slower than the gain ΔR(T), with much longer timescales (vertical dashed lines). (C,D) Analysis of an Izhikevich neuron in chattering mode with constant input and small voltage fluctuations. The neuron fires in regular bursts of activity. (D) Time-lagged measures C(T) and L(T) measure both, intra- (T < 10 ms) and inter-burst (T > 10 ms) dependencies, which decay very slowly due to regularity of the firing. The gain ΔR(T) reflects that most spiking can already be predicted from intra-burst dependencies, whereas inter-burst dependencies are highly redundant. In this case, only ΔR(T) yields a sensible time scale (blue dashed line). (E,F) Analysis of a generalized leaky integrate-and-fire neuron with long-lasting adaptation filter ξ [3, 43] and constant input. (F) Here, ΔR(T) decays slower to zero than the autocorrelation C(T), and is higher than L(T) for long T. Therefore, the dependence on past spikes is stronger when taking more-recent past spikes into account (ΔR(T)), as when considering them independently (L(T)). Due to these synergistic past dependencies, ΔR(T) is the only measure that captures the long-range nature of the spike adaptation.

https://doi.org/10.1371/journal.pcbi.1008927.g004

Second, we demonstrate the difference of history dependence to time-lagged measures on an Izhikevich neuron, which is a flexible model that can produce different neural firing patterns similar to those observed for real neurons [44]. Here, parameters were chosen according to the “chattering mode” [44], with constant input and small voltage fluctuations (Materials and methods). The neuron fires in regular bursts of activity, with consistent timing between spikes within and between bursts (Fig 4C). While time-lagged measures capture all the regularities in spiking and oscillate with the bursts of activity, history dependence correctly captures that spiking can almost be entirely predicted from intra-burst dependencies alone (Fig 4D). History dependence dismisses the redundant inter-burst dependencies and thereby yields a sensible measure of a timescale (blue dashed line).

Finally, we analyzed a generalized leaky integrate-and-fire neuron with long-range spike adaptation (22 seconds) (Fig 4E), which reproduces spike-frequency adaptation as observed for somatosensory pyramidal neurons [3, 43]. For this model, time-lagged measures C(T) and L(T) actually decay to zero much faster than the gain in history dependence ΔR(T), which is the only measure that captures the long-range adaptation effects of the model (Fig 4F). This shows that past dependencies in this model include synergistic effects, where the dependence is stronger in the context of more-recent spikes. This is most likely due to the non-linearity of the model, where past spikes cause a different adaptation when taken together as when considered as the sum of their contributions.

Thus, due to its ability to dismiss redundant past dependencies and to capture synergistic effects, history dependence really provides a complementary characterization of past dependencies compared to time-lagged measures. Importantly, because the approach better disentangles the effects of timescale and total history dependence, the results remain interpretable for very different models, and provide a more comprehensive view on past dependencies.

Without regularization, history dependence is severely overestimated for high-dimensional embeddings.

For demonstration, we estimated the history dependence R(τ, d) for varying numbers of bins d and a constant bin size τ = 20 ms (i.e. κ = 0 and T = d ⋅ τ). We compared estimates obtained by maximum likelihood (ML) estimation [28], or Bayesian estimation using the NSB estimator [33], with the model’s true R(τ, d) (Fig 5A). Both estimators accurately estimate R(τ, d) for up to d ≈ 20 past bins. As expected, the NSB estimator starts to be biased at higher d than the ML estimator. For embedding dimensions d > 30, both estimators severely overestimate R(τ, d). Note that ± two standard deviations are plotted as shaded areas, but are too small to be visible. Therefore, any deviation of estimates from the model’s true history dependence R(τ, d) can be attributed to positive estimation bias, i.e. a systematic overestimation of the true history dependence due to limited data.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Embedding optimization captures history dependence for a neuron model with long-lasting spike adaptation.

Results are shown for a generalized leaky integrate-and-fire (GLIF) model with long-lasting spike frequency adaptation [3, 43] with a temporal depth of one second (Methods and material). (A) For illustration, history dependence R(τ, d) was estimated on a simulated 90 minute recording for different embedding dimensions d and a fixed bin width τ = 20 ms. Maximum likelihood (ML) [28] and Bayesian (NSB) [33] estimators display the insufficient embedding versus estimation bias trade-off: For small embedding dimensions d, the estimated history dependence is much smaller, but agrees well with the true history dependence R(τ, d) for the given embedding. For larger d, the estimated history dependence increases, but when d is too high (d > 20), it severely overestimates the true R(τ, d). The Bayesian bias criterion (BBC) selects NSB estimates for which the difference between ML and NSB estimate is small (red solid line). All selected estimates are unbiased and agree well with the true R(τ, d) (grey line). Embedding optimization selects the highest, yet unbiased estimate (red diamond). (B) The Shuffling estimator (blue solid line) subtracts estimation bias on surrogate data (yellow dashed line) from the ML estimator (blue dashed line). Since the surrogate bias is higher than the systematic overestimation in the ML estimator (difference between grey and blue dashed lines), the Shuffling estimator is a lower bound to R(τ, d). Embedding optimization selects the highest estimate, which is still a lower bound (blue diamond). For A and B, shaded areas indicate ± two standard deviations obtained from 50 repeated simulations, which are very small and thus hardly visible. (C) Embedding-optimized BBC estimates (red line) yield accurate estimates of the model neuron’s true history dependence R(T), total history dependence R_tot and information timescale τ_R (horizontal and vertical dashed lines). The zoom-in (right panel) shows robustness of both regularization methods: For all T the model neuron’s R(T, d*, κ*) lies within errorbars (BBC), or consistently above the Shuffling estimator that provides a lower bound. Here, the model’s R(T, d*, κ*) was computed for the optimized embedding parameters d*, κ* that were selected via BBC or Shuffling, respectively (dashed lines). Shaded areas indicate ± two standard deviations obtained by bootstrapping, and colored vertical bars indicate past ranges over which estimates were averaged to compute (Materials and methods).

https://doi.org/10.1371/journal.pcbi.1008927.g005

The aim is now to identify the largest embedding dimension d* for which the estimate of R(τ, d) is not yet biased. A biased estimate is expected as soon as the two estimates ML and NSB start to differ significantly from each other (Fig 5A, red diamond), which is formalized by the Bayesian bias criterion (BBC) (Materials and methods). According to the BBC, all NSB estimates with d lower or equal to d* are unbiased (solid red line). We find that indeed all BBC estimates agree well with the true R(τ, d) (grey line), but d* yields the largest unbiased estimate.

The problem of estimation bias has also been addressed previously by the so-called Shuffling estimator [31]. The Shuffling estimator is based on the ML estimator and applies a bias correction term (Fig 5B). In detail, one approximates the estimation bias using surrogate data, which are obtained by shuffling of the embedded spiking history. The surrogate estimation bias (yellow dashed line) is proven to be larger than the actual estimation bias (difference between grey solid and blue dashed line). Therefore, Shuffling estimates provide lower bounds to the true history dependence R(τ, d). As with the BBC, one can safely maximize Shuffling estimates over d to find the embedding dimension d* that provides the largest lower bound to the model’s total history dependence R_tot (Fig 5B, blue diamond).

Thus, using a model neuron, we illustrated that history dependence can be severely overestimated if the embedding is chosen too complex. Only when overestimation is tamed by one of the two regularization methods, BBC or Shuffling, embedding parameters can be safely optimized to yield better estimates of history dependence.

Optimized embeddings capture the model’s true history dependence.

In the previous part, we demonstrated how embedding parameters are optimized for the example of fixed κ and τ. Now, we optimize all embedding parameters for fixed past range T to obtain embedding-optimized estimates of R(T). We find that embedding-optimized BBC estimates agree well with the true R(T), hence the model’s total history dependence R_tot and information timescale τ_R are well estimated (Fig 5C, vertical and horizontal dashed lines). In contrast, the Shuffling estimator underestimates the true R(T) for past ranges T > 200 ms, hence the model’s R_tot and τ_R are underestimated (blue dashed lines). For large past ranges T > 1000 ms, estimates of both estimators decrease again, because no additional history dependence is uncovered, whereas the constraint of an unbiased estimation decreases the temporal resolution of the embedding.

Embedding-optimized estimates are robust to overestimation despite maximization over complex embeddings.

In the previous part, we investigated how much of the true history dependence for different past ranges T (grey solid line) we miss by embedding the spiking history. An additional source of error is the estimation of history dependence from limited data. In particular, estimates are prone to overestimate history dependence systematically (Fig 5A and 5B).

To test explicitly for overestimation, we computed the true history dependence R(T, d*, κ*) for exactly the same sets of embedding parameters T, d*, κ* that were found during embedding optimization with BBC (grey dash-dotted line), and the Shuffling estimator (grey dotted line, Fig 5C, zoom-in). We expect that BBC estimates are unbiased, i.e. the true history dependence should lie within errorbars of the BBC estimates (red shaded area) for a given T. In contrast, Shuffling estimates are a lower bound, i.e. estimates should lie below the true history dependence (given the same T, d*, κ*). We find that this is indeed the case for all T. Note that this is a strong result, because it requires that the regularization methods work reliably for every single set of embedding parameters used for optimization—otherwise, parameters that cause overestimation would be selected.

Thus, we can confirm that the embedding-optimized estimates do not systematically overestimate the model neuron’s history dependence, and are on average lower bounds to the true history dependence. This is important for the interpretation of the results.

Mild overfitting can occur during embedding optimization on short recordings, but can be overcome with cross-validation.

We also tested whether the recording length affects the reliability of embedding-optimized estimates, and found very mild overestimation (1–3%) of history dependence for BBC for recordings as short as 3 minutes (S1 and S4 Figs). The overestimation is a consequence of overfitting during embedding optimization: Variance in the estimates increases for shorter recordings, hence maximizing over estimates selects embedding parameters that have high history dependence by chance. Therefore, the overestimation can be overcome by cross-validation, e.g. by optimizing embedding parameters on the first half, and computing estimates on the second half of the data (S1 Fig). In contrast, we found that for the model neuron, Shuffling estimates do not overestimate the true history dependence even for recordings as short as 3 minutes (S1 Fig). This is because the effect of overfitting was small compared to the systematic underestimation of Shuffling estimates. Here, all experimental recordings where we apply BBC are long enough (≈ 90 minutes), thus no cross-validation was applied on the experimental data.

Estimates of the information timescale are sensitive to the recording length.

Finally, we also tested the impact of the recording length on estimates of the total history dependence as well as estimates of the information timescale. While on recordings of 3 minutes embedding optimization still estimated ≈ 95% of the true R_tot, estimates were only ≈ 75% of the true τ_R (S2 Fig). Thus, estimates of the information timescale τ_R are more sensitive to the recording length, because they depend on the small additional contributions to R(T) for high past ranges T, which are hard to estimate for short recordings. Therefore, we advice to analyze recordings of similar length to make results on τ_R comparable across experiments. In the following, we explicitly shorten some recordings such that all recordings have approximately the same recording length.

In conclusion, embedding optimization accurately estimated the model neuron’s true history dependence. Moreover, for all past ranges, embedding-optimized estimates were robust to systematic overestimation. Embedding optimization is thus a promising approach to quantify history dependence and the information timescale in experimental spike recordings.

Embedding optimization reveals history dependence that is not captured by a generalized linear model or a single past bin.

We use embedding optimization to estimate history dependence R(T) as a function of the past range T (see Fig 6B for an example single unit from hippocampus layer CA1, and S6, S7 and S8 Figs for all analyzed sorted units). In this example, BBC and Shuffling with a maximum of d_max = 20 past bins led to very similar estimates for all T. Notably, embedding optimization with both regularization methods estimated high total history dependence of almost R_tot ≈ 40% with a temporal depth of almost a second, and an information timescale of τ_R ≈ 100 ms (Fig 6B). This indicates that embedding-optimized estimates capture a substantial part of history dependence also in experimental spike recordings.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Embedding optimization is key to estimate long-lasting history dependence in extra-cellular spike recordings.

(A) Example of recorded spiking activity in rat dorsal hippocampus layer CA1. (B) Estimates of history dependence R(T) for various estimators, as well as estimates of the total history dependence R_tot and information timescale τ_R (dashed lines) (example single unit from CA1). Embedding optimization with BBC (red) and Shuffling (blue) for d_max = 20 yields consistent estimates. Embedding-optimized Shuffling estimates with a maximum of d_max = 5 past bins (green) are very similar to estimates obtained with d_max = 20 (blue). In contrast, using a single past bin (d_max = 1, yellow), or fitting a GLM for the temporal depth found with BBC (violet dot), estimates much lower total history dependence. Shaded areas indicate ± two standard deviations obtained by bootstrapping, and small vertical bars indicate past ranges over which estimates of R(T) were averaged to estimate R_tot (Materials and methods). (C) An exponential past embedding is crucial to capture history dependence for high past ranges T. For T > 100 ms, uniform embeddings strongly underestimate history dependence. Shown is the median of embedding-optimized estimates of R(T) with uniform embeddings, relative to estimates obtained by optimizing exponential embeddings, for BBC with d_max = 20 (red) and Shuffling with d_max = 20 (blue) and d_max = 5 (green). Shaded areas show 95% percentiles. Median and percentiles were computed over analyzed sorted units in CA1 (n = 28). (D) Comparison of total history dependence R_tot for different estimation and embedding techniques for three different experimental recordings. For each sorted unit (grey dots), estimates are plotted relative to embedding-optimized estimates for BBC and d_max = 20. Embedding optimization with Shuffling and d_max = 20 yields consistent but slightly higher estimates than BBC. Strikingly, Shuffling estimates for as little as d_max = 5 past bins (green) capture more than 95% of the estimates for d_max = 20 (BBC). In contrast, estimates obtained by optimizing a single past bin, or fitting a GLM, are considerably lower. Bars indicate the median and lines indicate 95% bootstrap confidence intervals on the median over analyzed sorted units (CA1: n = 28; retina: n = 111; culture: n = 48).

https://doi.org/10.1371/journal.pcbi.1008927.g006

Importantly, other common estimation approaches fail to capture the same amount of history dependence (Fig 6B and 6D). To compare how well the different estimation approaches could capture the total history dependence, we plotted for each sorted unit the different estimates of R_tot relative to the corresponding BBC estimate (Fig 6D). Embedding optimization with Shuffling yields estimates that agree well with BBC estimates. The Shuffling estimator even yields slightly higher values on the experimental data. Interestingly, embedding optimization with the Shuffling estimator and as little as d_max = 5 past bins captures almost the same history dependence as BBC with d_max = 20, with a median above 95% for all neural systems. In contrast, we find that a single past bin only accounts for 70% to 80% of the total history dependence. A GLM bears little additional advantage with a slightly higher median of ≈ 85%. To save computation time, GLM estimates were only computed for the temporal depth that was estimated using BBC (Fig 6B, violet square). The remaining embedding parameters d and κ of the GLM’s history kernel were separately optimized using the Bayesian information criterion (Materials and methods). Since parameters were optimized, we argue that the GLM underestimates history dependence because of its specific model assumptions, i.e. no interactions between past spikes. Moreover, we found that the GLM performs worse than embedding optimization with only five past bins. Therefore, we conclude that for typical experimental spike trains, interactions between past spikes are important, but do not require very high temporal resolution. In the remainder of this paper we use the reduced approach (Shuffling d_max = 5) to compare history dependence among different neural systems.

Increasing bin sizes exponentially is crucial to estimate long-lasting history dependence.

To demonstrate this, we plotted embedding-optimized BBC estimates of R(T) using a uniform embedding, i.e. equal bin sizes, relative to estimates obtained with exponential embedding (Fig 6C), both for BBC with d_max = 20 (red) and Shuffling with d_max = 20 (blue) or d_max = 5 (green). For past ranges T > 100 ms, estimates using a uniform embedding miss considerable history dependence, which becomes more severe the longer the past range. In the case of d_max = 5, a uniform embedding captures around 80% for T = 1 s, and only around 60% for T = 5 s (median over analyzed sorted units in CA1). Therefore, we argue that an exponential embedding is crucial for estimating long-lasting history dependence.

Definition of history dependence.

We quantify history dependence R(T) as the mutual information I(X, X^−T) between present and past spiking X and X^−T, normalized by the binary Shannon information of spiking H(X), i.e. (6) Under the assumption of stationarity and ergodicity the mutual information can be computed either as the average over the stationary distribution p(x, x^−T), or the time average [21, 58], i.e. Here, indicates the presence of a spike in a small interval [t_n, t_n + Δt) with Δt = 5 ms throughout the paper, and encodes the spiking history in a time window [t_n − T, t_n) at times t_n = nΔt that are shifted by Δt.

Definition of lagged mutual information.

The lagged mutual information L(T) [41] for a stationary neural spike trains is defined as the mutual information between present spiking X and past spiking X_−T with time lag T, i.e. (11) (12) (13) Here, indicates the presence of a spike in a time bin [t_n, t_n + Δt) and the presence of a spike in a single past bin [t_n − T, t_n − T + Δt) at times t_n = nΔt that are shifted by Δt. In analogy to R(T), one can apply the generalized timescale to the lagged mutual information to obtain a timescale τ_L with (14)

Definition of autocorrelation.

The autocorrelation C(T) for a stationary neural spike train is defined as (15) with time lag T and and as above. For an exponentially decaying autocorrelation , τ_C is called autocorrelation time.

Past embedding.

Here, we encode the spiking history in a finite time window [t − T, t) as a binary sequence of binary spike counts in d past bins (Fig 2). When more than one spike can occur in a single bin, is chosen for spike counts larger than the median activity in the ith bin. This type of temporal binning is more generally referred to as past embedding. It is formally defined as a mapping (16) from the set of all possible spiking histories , i.e. the sigma algebra generated by the point process (neural spiking) in the time interval [t − T, t), to the set of d-dimensional binary sequences S^d. We can drop the dependence on the time t because we assume stationarity of the point process. Here, T is the embedded past range, d the embedding dimension, and θ denotes all the embedding parameters that govern the mapping, i.e. θ = (d, …). The resulting binary sequence at time t for given embedding θ and past range T will be denoted by In this paper, we consider the following two embeddings for the estimation of history dependence.

Uniform embedding.

If all bins have the same bin width τ = T/d, the embedding is called uniform. The main drawback of the uniform embedding is that higher past ranges T enforce a uniform decrease in resolution when d is fixed.

Exponential embedding.

One can generalize the uniform embedding by letting bin widths increase exponentially with bin index j = 1, …, d according to τ_j = τ₁10^(j−1)κ. Here, τ₁ gives the bin size of the first past bin, and is uniquely determined when T, d and κ are specified. Note that κ = 0 yields a uniform embedding, whereas κ > 0 decreases resolution on distant past spikes. For fixed embedding dimension d and past range T, this allows to retain a higher resolution on spikes in the more-recent past.

Sufficient embedding.

Ideally, the past embedding preserves all the information that the spiking history in the past range T has about the present spiking dynamics. In that case, no additional past information has an influence on the probability for x_t once the embedded spiking history is given, i.e. (17) for any other past embedding . If Eq (17) holds for all times t, the embedding Γ_T(θ) is called a sufficient embedding. For the remainder of this paper, the sequences of sufficient embeddings are denoted by .

Insufficient embeddings cause underestimation of history dependence.

The past embedding is essential when inferring history dependence from recordings, because an insufficient embedding causes underestimation of history dependence. To show this, we note that for any embedding parameters θ and past range T the Kullback-Leibler divergence between the spiking probability for the sufficient embedding and cannot be negative [60], i.e. (18) with equality iff . By taking the average over all times t_n, we arrive at (19) (20) (21) (22) where the last step follows from stationarity and ergodicity and marginalizing out in the first term. From here, it follows that one always underestimates the history dependence in neural spiking, as long as the embedding is not sufficient, i.e. (23)

Embedding optimization.

For given T, find the optimal embedding θ* that maximizes the estimated history dependence (24) This yields an embedding-optimized estimate of the true history dependence R(T).

Requirements.

Embedding optimization can only give sensible results if the optimized estimates are guaranteed to be unbiased or a lower bound to the true R(T, θ). Otherwise, embeddings will be chosen that strongly overestimate history dependence. In this paper, we therefore use two estimators, BBC and Shuffling, the former of which is designed to be unbiased, and the latter a lower bound to the true R(T, θ) (see below). In addition, embedding optimization works only if the estimation variance is sufficiently small. Otherwise, maximizing over variable estimates can lead to a mild overestimation. We found for a benchmark model that this overestimation was negligibly small for a recording length of 90 minutes for a model neuron with a 4 Hz average firing rate (S1 Fig). For smaller recording lengths, potential overfitting can be avoided by cross-validation, i.e. optimizing embeddings on one half of the recording and computing embedding-optimized estimates on the other half.

Implementation.

For the optimization, we compute estimates for a range of embedding dimensions d ∈ [1, 2, …, d_max] and scaling parameters κ = [0, …, κ_max]. For each T, we then choose the optimal parameter combination d*, κ* for each T that maximizes the estimated history dependence , and use as the best estimate of R(T).

Estimation of total history dependence and the information timescale.

When estimating history dependence R(T) from data, there are some adjustments required to estimate the total history dependence R_tot and the information timescale τ_R.

First, estimates are not guaranteed to converge for large past ranges T, but might decrease due to a reduced resolution of embeddings for higher T (Fig 2D). Thus, we estimated an interval [T_D, T_max] for which estimates have converged. Here, the temporal depth T_D and the upper bound T_max are the first and the last past ranges T for which estimates are within one standard deviation of the highest estimate , i.e. (Fig 2D, vertical blue bars). The standard deviation was estimated by bootstrapping (see Bootstrap confidence intervals). From this interval, an estimate of the total history dependence is obtained by averaging over past ranges T ∈ [T_D, T_max] (Fig 2D, horizontal dashed blue line).

Second, noisy estimates are not guaranteed to be monotonously increasing, hence increments can be negative. Moreover, noisy estimates can lead to positive even though the true R(T) has already converged to R_tot. This can have a huge effect on the estimated information timescale if one simply uses these estimates in Eq (5). To avoid this, we use knowledge about the behavior of the true R(T) when estimating ΔR(T). In particular, we set estimates equal to the largest previous estimate for T′ < T if they fall below it, and equal to if they are larger than . This enforces that the estimated gain is non-negative, and excludes spurious gain for high T due to noisy estimates.

Finally, the information timescale τ_R can crucially depend on the choice of the minimum past range T₀ in the sum in Eq (5). A T₀ > 0 larger than zero allows to ignore short term effects on the history dependence such as the refractory period or different firing modes, which we found beneficial for resolving differences in the timescale among different neural systems (S15 Fig). In contrast, if the decay is truly exponential, then τ_R is independent of T₀. In this paper, we chose T₀ = 10ms to exclude short term effects, while also not excluding too much past information.

Workflow.

The estimation workflow using embedding optimization is summarized in Fig 10.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Workflow of past-embedding optimization to estimate history dependence and the information timescale.

1) Define a set of embedding parameters d, κ for fixed past range T. 2) For each embedding d, κ, record sequences of current and past spiking for all time steps t_n in the recording. 3) Use the frequencies of the recorded sequences to estimate history dependence for each embedding, either using maximum likelihood (ML), or fully Bayesian estimation (NSB). 4) Apply regularization, i.e. the Bayesian bias criterion (BBC) or Shuffling bias correction, to ensure that all estimates are unbiased or lower bounds to the true history dependence. 5) Select the optimal embedding to obtain an embedding-optimized estimate of R(T). 6) Repeat the estimation for a set of past ranges T to compute estimates of the information timescale τ_R and the total history dependence R_tot.

https://doi.org/10.1371/journal.pcbi.1008927.g010

Maximum likelihood estimation.

Most commonly, probabilities of spike sequences a_k are then estimated as the relative frequencies of their occurrence in the observed data. It is the maximum likelihood (ML) estimator of π for the multinomial likelihood (28) Plugging the estimates into the definition of the entropy in Eq (27) results in the ML estimator of the entropy (29) or history dependence (30) The ML estimator has the right asymptotic properties [28, 61], but is known to underestimate the entropy severely when data is limited [28, 62]. This is because all probability mass is assumed to be concentrated on the observed outcomes. A more concentrated probability distribution results in a smaller entropy, in particular if many outcomes have not been observed. This results in a systematic underestimation or negative bias (31) The negative bias in the entropy, which is largest for the highest-dimensional joint entropy , then typically leads to severe overestimation of the mutual information and history dependence [27, 63]. Because of this severe overestimation, we cannot use the ML estimator for embedding optimization.

Bayesian Nemenman-Shafee-Bialek (NSB) estimator.

In a Bayesian framework, the entropy is estimated as the posterior mean or minimum mean square error (MMSE) (32) The posterior mean is the mean of the entropy with respect to the posterior distribution on the probability vector π given the observed frequencies of spike sequences n (33) The probability p(n|π) for i.i.d. observations n from an underlying distribution π is given by the multinomial distribution in Eq (28).

If the prior p(π) is a conjugate prior to the multinomial likelihood, then the high-dimensional integral of Eq (32) can be evaluated analytically [32]. This is true for a class of priors called Dirichlet priors, and in particular for symmetric Dirichlet priors (34) The prior p(π|β) gives every outcome the same a priori weight, but controls the weight β > 0 of uniform prior pseudo-counts. A β = 1 corresponds to a flat prior on all probability distributions π, whereas β → 0 gives maximum likelihood estimation (no prior pseudo-count).

It has been shown that the choice of β is highly informative with respect to the entropy, in particular when the number of outcomes K becomes large [64]. This is because the a priori variance of the entropy vanishes for K → ∞, thus for any π ∼ p(π|β) the entropy H(π) is very close to the a priori expected entropy (35) where are the polygamma functions. In addition, a lot of data is required to counter-balance this a priori expectation. This is because the prior adds pseudo-counts on every outcome, i.e. it assumes that every outcome has been observed β times prior to inference. In order to influence a prior that constitutes K pseudo-counts, one needs at least N > K samples, with more data required the sparser the true underlying distribution. Therefore, an estimator of the entropy for little data and fixed concentration parameter β is highly biased towards the a priori expected entropy ξ(β).

Nemenman et al. [33] exploited the tight link between concentration parameter β and the a priori expected entropy to derive a mixture prior (36) (37) that weights Dirichlet priors to be flat with respect to the expected entropy ξ(β). Since the variance of this expectation vanishes for K ≫ 1 [64], for high K the prior is also approximately flat with respect to the entropy, i.e. for π ∼ p_NSB(π). The resulting MMSE estimator for the entropy is referred to as the NSB estimator (38) (39) Here, ρ(β, n) is proportional to the evidence for given concentration parameter (40) (41) where Γ(x) is the gamma function. The posterior mean of the entropy for given concentration parameter is (42)

From the Bayesian entropy estimate, we obtain an NSB estimator for history dependence (43) where the marginal and joint entropies are estimated individually using the NSB method.

To compute the NSB entropy estimator, one has to perform a one-dimensional integral over all possible concentration parameters β. This is crucial to be unbiased with respect to the entropy. An implementation of the NSB estimator for Python3 is published alongside the paper with our toolbox [37]. To compute the integral, we use a Gaussian approximation around the maximum a posteriori β* to define sensible integration bounds when the likelihood is highly peaked, as proposed in [34].

Bayesian bias criterion.

The goal of the Bayesian bias criterion (BBC) is to indicate when estimates of history dependence are potentially biased. It might indicate bias even when estimates are unbiased, but the opposite should never be true.

To indicate a potential estimation bias, the BBC compares ML and NSB estimates of the history dependence. ML estimates are biased when too few joint sequences have been observed, such that the probability for unobserved or undersampled joint outcomes is underestimated. To counterbalance this effect, the NSB estimate adds β pseudo-counts to every outcome, and then infers β with an uninformative prior. For the BBC, we turn the idea around: when the assumption of no pseudo-counts (ML) versus a posterior belief on non-zero pseudo-counts (NSB) yield different estimates of history dependence, then too few sequences have been observed and estimates are potentially biased. This motivates the following definition of the BBC.

The NSB estimator R_NSB(T, θ) is biased with tolerance p > 0, if (44) Similarly, we define the BBC estimator (45) This estimator is designed to be unbiased, and thus can be used for embedding optimization in Eq (24). We use the NSB estimator for R(T, θ) instead of the ML estimator, because it is generally less biased. A tolerance p > 0 accounts for this, and accepts NSB estimates when there is only a small difference between the estimates. The bound for the difference is multiplied by , because this provides the scale on which one should be sensitive to estimation bias. We found that a tolerance of p = 0.05 was small enough to avoid overestimation by BBC estimates on the benchmark model (Fig 5 and S1 Fig).

Shuffling estimator.

The Shuffling estimator was originally proposed in [31] to reduce the sampling bias of the ML mutual information estimator. It has the desirable property that it is negatively biased in leading order of the inverse number of samples. Because of this property, Shuffling estimates can safely be maximized during embedding optimization without the risk of overestimation. Here, we therefore propose to use the Shuffling estimator for embedding-optimized estimation of history dependence.

The idea behind the Shuffling estimator is to rewrite the ML estimator of history dependence as (46) and to correct for bias in the entropy estimate . Since X is well sampled and thus is unbiased, and the bias of the ML entropy estimator is always negative [28, 61], we know that (47) (48) Therefore, if we find a correction term of the magnitude of , we can turn the bias in the estimate of the history dependence from positive to negative, thus obtaining an estimator that is a lower bound of the true history dependence. This can be achieved by subtracting a lower bound of the estimation bias from .

In the following, we describe how [31] obtain a lower bound of the bias in the conditional entropy by computing the estimation bias for shuffled surrogate data.

Surrogate data are created by shuffling recorded spike sequences such that statistical dependencies between past bins are eliminated. This is achieved by taking all past sequences that were followed by a spike, and permuting past observations of the same bin index j. The same is repeated for all past sequences that were followed by no spike. The underlying probability distribution can then be computed as (49) and the corresponding entropy is (50) The pairwise probabilities are well sampled, and thus each conditional entropy in the sum can be estimated with high precision. This way, the true conditional entropy for the shuffled surrogate data can be computed and compared to the ML estimate on the shuffled data. The difference between the two (51) yields a correction term that is on average equal to the bias of the ML estimator on the shuffled data.

Importantly, the bias of the ML estimator on the shuffled data is in leading order more negative than on the original data. To see this, we consider an expansion of the bias on the conditional entropy in inverse powers of the sample size N [27, 63] (52) Here, denotes the number of past sequences with nonzero probability of being observed when followed by a spike (x = 1) or no spike (x = 0), respectively. Notably, the bias is negative in leading order, and depends only on the number of possible sequences . For the shuffled surrogate data, we know that implies , but Shuffling may lead to novel sequences that have zero probability otherwise. Hence the number of possible sequences under Shuffling can only increase, i.e. , and thus the bias of the ML estimator under Shuffling to first order is always more negative than for the original data (53) Terms that could render it higher are of order and are assumed to have no practical relevance.

This motivates the following definition of the Shuffling estimator: Compute the difference between the ML estimator on the shuffled and original data to yield a bias-corrected Shuffling estimate (54) and use this to estimate history dependence (55) Because of Eqs (48) and (53), we know that this estimator is negatively biased in leading order (56) and can safely be used for embedding optimization.

Estimation of history dependence by fitting a generalized linear model (GLM).

Another approach to the estimation of history dependence is to model the dependence of neural spiking onto past spikes explicitly, and to fit model parameters to maximize the likelihood of the observed spiking activity [21]. For a given probability distribution of the model with parameters ν, the conditional entropy can be estimated as (57) which one can plug into Eq (6) to obtain an estimate of the history dependence. The strong law of large numbers [59] ensures that if the model is correct, i.e. for all t, this estimator converges to the entropy H(X|X^−T) for N → ∞. However, any deviations from the true distribution due to an incorrect model will lead to an underestimation of history dependence, similar to choosing an insufficient embedding. Therefore, model parameters should be chosen to maximize the history dependence, or to maximize the likelihood (58) We here consider a generalized linear model (GLM) with exponential link function that has successfully been applied to make predictions in neural spiking data [20] and can be used for the estimation of directed, causal information [21]. In a GLM with past dependencies, the spiking probability at time t is described by the instantaneous rate or conditional intensity function (59) Since we discretize spiking activity in time as spiking or non-spiking in a small time window Δt, the spiking probability is given by the binomial probability (60) The idea of the GLM is that past events contribute independently to the probability of spiking, such that the conditional intensity function factorizes over their contributions. Hence, it can be written as (61) where h_j gives the contribution of past activity in past time bin j to the firing rate, and μ is an offset that is adapted to match the average firing rate.

Although fitting GLM parameters is more data-efficient than computing non-parametric estimates, overfitting may occur for limited data and high embedding dimensions d, hence d cannot be chosen arbitrarily high. In order to estimate a maximum of history dependence for limited d, we apply the same type of binary past embedding as we use for the other estimators, and optimize the embedding parameters by minimizing the Bayesian information criterion [65]. In particular, for given past range T, we choose embedding parameters d*, κ* that minimize (62) where N is the number of samples and (63) is the maximized log-likelihood of the recorded spike sequences for optimal model parameters μ*, h*. We then use the optimized embedding parameters to estimate the conditional entropy according to (64) which results in the GLM estimator of history dependence (65)

Bootstrap confidence intervals.

In order to estimate confidence intervals of estimates for given past embeddings, we apply the blocks of blocks bootstrap method [66]. To obtain bootstrap samples, we first compute all the binary sequences for n = 1, …, N that result from discretizing the spike recording in N time steps Δt and applying the past embedding. We then randomly draw N/l blocks of length l of the recorded binary sequences such that the total number of redrawn sequences is the same as the in the original data. We choose l to be the average interspike interval (ISI) in units of time steps Δt, i.e. l = 1/(rΔt) with average firing rate r. Sampling successive sequences over the typical ISI ensures that bootstrap samples are representative of the original data, while also providing a high number of distinct blocks that can be drawn.

The different estimators (but not the bias criterion) are then applied to each bootstrap sample to obtain confidence intervals of the estimates. Instead of computing the 95% confidence interval via the 2.5 and 97.5 percentiles of the bootstrapped estimates, we assumed a Gaussian distribution and approximated the interval via , where is the standard deviation over the bootstrapped estimates.

We found that the true standard deviation of estimates for the model neuron was well estimated by the bootstrapping procedure, irrespective of the recording length (S10 Fig). Furthermore, we simulated 100 recordings of the same recording length, and for each computed confidence interval for the past range T with the highest estimated history dependence R(T). By measuring how often the model’s true value for the same embedding was included in these intervals, we found that the Gaussian confidence intervals are indeed close to the claimed confidence level (S10 Fig). This indicates that the bootstrap confidence intervals approximate well the uncertainty associated with estimates of history dependence.

Cross-validation.

For small recording lengths, embedding optimization may cause overfitting through the maximization of variable estimates (S1 Fig). To avoid this type of overestimation, we apply one round of cross-validation, i.e. we optimize embeddings over the first half of the recording, and evaluate estimates for the optimal past embedding on the second half. We chose this separation of training and evaluation data sets, because it allows the fastest computation of binary sequences for the different embeddings during optimization. We found that none of the cross-validated embedding-optimized estimates were systematically overestimating the true history dependence for the benchmark model for recordings as short as three minutes (S1 Fig). Therefore, cross-validation allows to apply embedding optimization to estimate history dependence even for short recordings.

Generalized leaky integrate-and-fire neuron with spike-frequency adaptation.

As a benchmark model, we chose a generalized leaky integrate-and-fire model (GLIF) with an additional adaptation filter ξ (GLIF-ξ) that captures spike-frequency adaptation over 20 seconds [43].

For a standard leaky integrate-and-fire neuron, the neuron’s membrane is formalized as an RC circuit, where the cell’s lipid membrane is modeled as a capacitance C, and the ion channels as a resistance that admits a leak current with effective conductance g_L. Hence, the temporal evolution of the membrane’s voltage V is governed by (66) Here, V_R denotes the resting potential and I_ext(t) external currents that are induced by some external drive. The neuron emits an action potential (spike) once the neuron crosses a voltage threshold V_T, where a spike is described as a delta pulse at the time of emission . After spike emission, the neuron returns to a reset potential V₀. Here, we do not incorporate an explicit refractory period, because interspike intervals in the simulation were all larger than 10 ms. For constant input current I_ext, integrating Eq (66) yields the membrane potential between two spiking events (67) where is the time of the most recent spike, γ = g_L/C the inverse membrane timescale and V_∞ = V_R + I_ext/γ the equilibrium potential.

In contrast to the LIF, the GLIF models the spike emission with a soft spiking threshold. To do that, spiking is described by an inhomogeneous Poisson process, where the spiking probability in a time window of width δt ≪ 1 is given by (68) Here, the spiking probability is governed by the time dependent firing rate (69) The idea is that once the membrane potential V(t) approaches the firing threshold V_T(t), the firing probability increases exponentially, where the exponential increase is modulated by 1/ΔV. For ΔV → 0, we recover the deterministic LIF, while for larger ΔV the emission becomes increasingly random.

In the GLIF-ξ, the otherwise constant threshold is modulated by the neuron’s own past activity according to (70) Thus, depending on their spike times , emitted action potentials increase or decrease the threshold additively and independently according to an adaptation filter ξ(t). In the experiments conducted in [43], the following functional form for the adaptation filter was extracted: (71) The filter is an effective model not only for the measured increase in firing threshold, but also for spike-triggered currents that reduce the membrane potential. When mapped to the effective adaptation filter ξ, it turned out that past spikes lead to a decrease in firing probability that is approximately constant over a period T_ξ = 8.3 ms, after which it decays like a power-law with exponent β_ξ = 0.93, until the contributions are set to zero after 22 s.

Model variant with 1s past kernel.

For demonstration, we also simulated a variant of the above model with a 1s past kernel (72) All parameters are identical apart from the strength of the kernel , which was adapted to maintain a firing rate of 4 Hz despite the shorter kernel.

Simulation details.

In order to ensure stationarity, we simulated the model neuron exposed to a constant external current I_ext = const. over a total duration of T_rec = 900 min. Thereby, the current I_ext was chosen such that the neuron fired with a realistic average firing rate of 4 Hz. During the simulation, Eq (66) was integrated using simple Runge-Kutta integration with an integration time step of δt = 0.5 ms. At every time step, random spiking was modeled as a binary variable with probability as in Eq (68). After a burning-in time of 100 s, spike times were recorded and used for the estimation of history dependence. The detailed simulation parameters can be found in Table 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Simulation parameters of the GLIF-ξ model.

https://doi.org/10.1371/journal.pcbi.1008927.t001

Computation of the total history dependence.

In order to determine the total history dependence in the simulated spiking activity, we computed the conditional entropy H(X|X^−∞) from the conditional spiking probability in Eq (68) that was used for the simulation. Note that this is only possible because of the constant input current, otherwise the conditional spiking probability would also capture information about the external input.

Since the conditional probability of spiking used in the simulation computes the probability in a simulation step δt = 0.5 ms, we first have to transform this to a probability of spiking in the analysis time step Δt = 5 ms. To do so, we compute the probability of no spike in a time step [t, t + Δt) according to (73) and then compute the probability of at least one spike by . Here, the rate is computed as λ(t) in Eq (69), but only with respect to past spikes that are emitted at times . This is because no spike that occurs within [t, t + Δt) must be considered when computing .

For sufficiently long simulations, one can make use of the SLLN to compute the conditional entropy (74) and thus the total history dependence (75) which gives an upper bound to the history dependence for any past embedding.

Computation of history dependence for given past embedding.

To compute history dependence for given past embedding, we use that the model neuron can be well approximated by a generalized linear model (GLM) within the parameter regime of our simulation. We can thus fit a GLM to the simulated data for the given past embedding T, d, κ to obtain a good approximation of the corresponding true history dependence R(T, d, κ). Note that this is a specific property if this model and does not hold in general. For example in experiments, we found that the GLM accounted for less history dependence than model-free estimates (Fig 6).

To map the model neuron to a GLM, we plug the membrane and threshold dynamics of Eqs (67) and (70) into the equation for the firing rate Eq (69), i.e. (76) For the parameters used in the simulation, the decay time of the reset term V₀ − V_∞ is 1/γ = 15.3 ms. When compared to the minimum and mean inter-spike intervals of ISI_min = 25 ms and , it is apparent that the probability for two spikes to occur within the decay time window is negligibly small. Therefore, one can safely approximate (77) i.e. describing the potential reset after a spike as independent of other past spikes, because contributions beyond the last spike (j > 0) are effectively zero. Using the above approximation, one can formulate the rate as in a generalized linear model with (78) where (79) (80) and indicates whether the neuron spiked in [t − jδt, t − (j − 1)δt). Therefore, the true spiking probability of the model is well described by a GLM.

We use this relation to approximate the history dependence R(T, d, κ) for any past embedding T, d, κ with a GLM with the same past embedding. Since in that case the parameters μ and h are not known, we fitted them to the simulated 900 minute recording via maximum likelihood (see above) and computed the history dependence according to (81)

Computation of history dependence as a function of the past range.

To approximate the model’s true history dependence R(T), for each T we computed GLM estimates (Eq 81) for a varying number of past bins d ∈ [25, 50, 75, 100, 125, 150]. For each d, the scaling κ was chosen such that the size of the first past bin was equal or less than 0.5 ms. To save computation time, and to reduce the effect of overfitting, the GLM parameters where fitted on 300 minutes of the simulation, whereas estimates were computed on the full 900 minutes of the simulated recording. For each T, we then chose the highest estimate among the estimates for different d as the best estimate of the true R(T).

Rat cortical culture.

Neurons were extracted from rat cortex (1st day postpartum) and recorded in vitro on an electrode array 2–3 weeks after plating day. We took data from five consecutive sessions (L_Prg035_txt_nounstim.txt, L_Prg036_txt_nounstim.txt, …, L_Prg039_txt_nounstim.txt) with a total duration of about T_rec ≈ 203 min. However, we only analyzed the first 90 minutes to make the results comparable to the other neural systems. We analyzed in total n = 48 sorted units that satisfied our requirement on the firing rate. More details on the recording procedure can be found in [67], and details on the data set proper can be found in [49].

Salamander retina.

Spikes from larval tiger salamander retinal ganglion cells were recorded in vitro by extracting the entire retina on an electrode array [68], while a non-repeated natural movie (leaves moving in the wind) was projected onto the retina. The recording had a total length of about T_rec ≈ 82 min, and we analyzed in total n = 111 sorted units that satisfied our requirement on the firing rate. More details on the recording procedure and the data set can be found in [47, 48]. The spike recording was obtained from the Dryad database [47].

Rat dorsal hippocampus (layer CA1).

We evaluated spike trains from a multichannel simultaneous recording made from layer CA1 of the right dorsal hippocampus of a Long-Evans rat during an open field task (data set ec014.277). The data set provided sorted spikes from 8 shanks with 64 channels. The recording had a total length of about T_rec ≈ 90 min. We analyzed in total n = 28 sorted units that were indicated as single units and satisfied our requirement on the firing rate. More details on the experimental procedure and the data set can be found in [45, 46]. The spike recording was obtained from the NSF-founded CRCNS data sharing website.

Mouse primary visual cortex.

Neurons were recorded in vivo during spontaneous behavior, while face expressions were monitored. Recordings were obtained by 8 simultaneously implanted Neuropixel probes, and sorted units were located using the location of the electrode contacts provided in [50], and the Allen Mouse Common Coordinate Framework [69]. We analyzed in total n = 142 sorted units from the mouse “Waksman” that belonged to primary visual cortex (irrespective of their layer) and satisfied our requirement on the firing rate. Second, we only selected units that were recorded for more than T_rec ≈ 40 min (difference between the last and first recorded spike time). Details on the recording procedure and the data set can be found in [58] and [50].

Details to Fig 3.

For Fig 3B, the process was considered for l = 1 and an reactivation probability of m = 0.8. For l = 1, all probabilities can easily be calculated, with marginal probability to be active p(x_t = 1) = h/(1 − m + mh), and conditional probabilities p(x_t = 1|x_t−1 = 1) = h + (1 − h)m and p(x_t = 1|x_t−1 = 0) = h. From these probabilities, the total mutual information I_tot and total history dependence R_tot could be directly computed. We then plotted these quantities as a function of h, where values of h were chosen to vary the firing rate between 0.5 and 10 Hz, with a bin size of Δt = 5 ms. For Fig 3C, the binary autoregressive process was simulated for n = 10⁷ time steps with m = 0.8 (l = 1), whereas for l = 5, m was adapted to yield approximately the same R_tot as for l = 1. The input activation probability h was chosen to lead to a fixed probability p(x = 1) ≈ 0.025, corresponding to 5 Hz firing rate with Δt = 5 ms. Autocorrelation C(T) was computed using the MR.estimator toolbox [52], and ΔR(T) and L(T) were estimated using plugin estimation. For Fig 3D, the same procedures were applied as in Fig 3C, but now m was varied between 0.5 and 0.95, and h was adapted for each m to hold the firing rate fixed at 5 Hz. For Fig 3E, the same procedures were applied as in Fig 3C, but now l was varied between 1 and 10, and h and m were adapted for each l to hold the firing rate fixed at 5 Hz and R_tot fixed at the value for l = 1 and m = 0.8.

Details to Fig 4A and 4B.

The branching process was simulated using the MR.estimator toolbox, with a time step of Δt = 4ms, population rate of 500 Hz and subsampling probability of 0.01. Thus, the subsampled spike train had a firing rate of ≈ 5 Hz. The branching parameter was set to m = 0.98 with analytic autocorrelation time τ_C(m) = 198ms. For a long simulation, autocorrelation C(T) was computed using the MR.estimator toolbox, L(T) using plugin estimation, and R(T) using the embedding-optimized Shuffling estimator with d_max = 20. The generalized timescales τ_R and τ_L were computed with T₀ = 10ms.

Details to Fig 4C and 4D.

The Izhikevich model was simulated with the PyNN toolbox [70], with parameters set to the chattering mode (a = 0.02, b = 0.2, c = −50, d = 2), simulation time bin dt = 0.01 ms, and noisy input with mean 0.011 and standard deviation 0.001. For the analysis, a time step of Δt = 1 ms was chosen. Apart from that, C(T) and L(T) were computed as for Fig 4B. Here, R(T) was computed with BBC and d_max = 20, which revealed higher R_tot than Shuffling. To compute τ_R, we set T₀ = 0.

Details to Fig 4E and 4F.

The GLIF model was simulated as described in Benchmark neuron model (model with 22s past kernel). The analysis time step was Δt = 5 ms. Apart from that, C(T) and L(T) were computed as for Fig 4B. History dependence R(T) was estimated using a GLM as described in Benchmark neuron model. To compute τ_R, we set T₀ = 10 ms.

Details to Fig 5A and 5B.

In Fig 5A and 5B, we applied the ML, NSB, BBC and Shuffling estimators of R(τ, d) to a simulated recording of 90 minutes. Embedding parameters were T = d ⋅ τ and κ = 0, with τ = 20 ms and d ∈ [1, 60]. Since the goal was to show the properties of the estimators, confidence intervals were estimated from 50 repeated 90 minute simulations instead of bootstrap samples from the same recording. Each simulation had a burning in period of 100 seconds. To estimate the true R(τ, d), a GLM was fitted on a 300 minute recording and evaluated on the full 900 minute recording for the estimation of R.

Details to Fig 5C.

In Fig 5C, history dependence R(T) was estimated on a 90 minute recording for 57 different values of T in a range T ∈ [10 ms, 3 s]. Embedding-optimized estimates were computed with up to d_max = 25 past bins, and 95% confidence intervals were computed using the standard deviation over n = 100 bootstrap samples (see Bootstrap confidence intervals). To estimate the true R(T, d*, κ*) for the optimized embedding parameters d*, κ* with either BBC or Shuffling, a GLM was fitted for the same embedding parameters on a 300 minute recording and evaluated on 900 minutes recording for the estimation of R.

Details to Fig 6.

For Fig 6, history dependence R(T) was estimated for 61 different values of T in a range T ∈ [10 ms, 5 s]. For each recording, we only analyzed the first 90 minutes to have a comparable recording length. For embedding optimization, we used d_max = 20 as a default for BBC and Shuffling, and compared the estimates with the Shuffling estimator optimized for d_max = 5 (max five bins) and d_max = 1 (one bin). For the GLM, we only estimated R(T_D) for the temporal depth T_D that was estimated with BBC. To optimize the estimate, we computed GLM estimates of R(T_D) with the optimal embedding found by BBC, and for varying embedding dimension d ∈ [1, 2, 3, ‥, 20, 25, 30, 35, 40, 45, 50], where for each d we chose κ such that τ₁ = Δt. We then chose the embedding that minimized the BIC, and took the corresponding estimate as a best estimate for R_tot. For Fig 6A, we plotted only spike trains of channels that were identified as single units. For Fig 6B, 95% confidence intervals were computed using the standard deviation over n = 100 bootstrap samples. For Fig 6C, embedding-optimized estimates with uniform embedding (κ = 0) were computed with d_max = 20 (BBC and Shuffling) or d_max = 5 (Shuffling). Medians were computed over the n = 28 sorted units in CA1.

Details to Figs 7 and 8.

For Figs 7 and 8, history dependence was R(T) was estimated for 61 different values of T in a range T ∈ [10 ms, 5 s] using the Shuffling estimator with d_max = 5. The autocorrelation coefficients C(T) were computed with the MR.Estimator toolbox [52], and the autocorrelation time τ_C was obtained using the exponential_offset fitting function. For each recording, we only analyzed the first 40 minutes to have a comparable recording length. For Fig 7, medians of τ_R, τ_C and R_tot were computed over all sorted units that were analyzed, and 95% confidence intervals on the medians were obtained by bootstrapping with n = 10000 resamples of the median. For Fig 8, 95% confidence intervals were computed using the standard deviation over n = 100 bootstrap samples.

1) The embedding of past spiking activity should be individually optimized to account for very different spiking statistics.

It is crucial to optimize the embedding for each neuron individually, because history dependence can strongly differ for neurons from different areas or neural systems (Fig 7), or even among neurons within a single area (see examples in Fig 8). Individual optimization enables a meaningful comparison of information timescale and history dependency R between neurons.

2) The estimation has to capture any non-linear or higher-order statistical dependencies.

Embedding optimization using both, the BBC or Shuffling estimators, is based on non-parametric estimation, in which the joint probabilities of current and past spiking are directly estimated from data. Thereby, it can account for any higher-order or non-linear dependency among all bins. In contrast, the classical generalized linear model (GLM) that is commonly used to model statistical dependencies in neural spiking activity [20, 21] does not account for higher-order dependencies. We found that the GLM recovered consistently less total history dependence R_tot (Fig 6D). Hence, to capture single-neuron history dependence, higher-order and non-linear dependencies are important, and thus a non-parametric approach is advantageous.

3) Estimation has to be computationally feasible even for a high number of recorded neurons.

Strikingly, while higher-order and non-linear dependencies are important, the estimation of history dependence does not require high temporal resolution. Optimizing up to d_max = 5 past bins with variable exponential scaling κ could account for most of the total history dependence that was estimated with up to d_max = 20 bins (Fig 6D). With this reduced setup, embedding optimization is feasible within reasonable computation time. Computing embedding-optimized estimates of the history dependence R(T) for 61 different values of T (for 40 minute recordings, the approach used for Figs 7 and 8) took around 10 minutes for the Shuffling estimator, and about 8.5 minutes for the BBC per neuron on a single computing node. Therefore, we recommend using d_max = 5 past bins when computation time is a constraint. Ideally, however, one should check for a few recordings if higher choices of d_max lead to different results, in order to cross-validate the choice of d_max = 5 for the given data set.

4) Estimates have to be reliable lower bounds, otherwise one cannot interpret the results.

It is required that embedding-optimized estimates do not systematically overestimate history dependence for any given embedding. Otherwise, one cannot guarantee that on average estimates are lower bounds to the total history dependence, and that an increase in history dependence for higher past ranges is not simply caused by overestimation. This guarantee is an important aspect for the interpretation of the results.

For BBC, we found that embedding-optimized estimates are unbiased if the variance of estimators is sufficiently small (S1 Fig). The variance was sufficiently small for recordings of 90 minutes duration. When the variance was too high (short recordings with 3–45 minutes recording length), maximizing estimates for different embedding parameters introduced very mild overestimation due to overfitting (1–3%) (S1 Fig). The overfitting can, however, be avoided by cross-validation, i.e. optimizing the embedding on one half of the recording and computing estimates on the other half. Using cross-validation, we found that embedding-optimized BBC estimates were unbiased even for recordings as short as 3 minutes (S1 Fig).

For Shuffling, we also observed overfitting, but the overestimation was small compared to the inherent systematic underestimation of Shuffling estimates. Therefore, we observed no systematic overestimation by embedding-optimized Shuffling estimates on the model neuron, even for shorter recordings (3 minutes and more). Thus, for the Shuffling estimator, we advice to apply the estimator without cross-validation as long as recordings are sufficiently long (10 minutes and more, see next point).

5) Spike recordings must be sufficiently long (at least 10 minutes), and of similar length, in order to allow for a meaningful comparison of total history dependence and information timescale across experiments.

The recording length affects estimates of the total history dependence R_tot, and especially of the information timescale τ_R. This is because more data allow more-complex embeddings, hence higher history dependence can be estimated. Moreover, complex embeddings are particular relevant for long past ranges T. Therefore, if recordings are shorter, smaller R(T) will be estimated for long past ranges T, leading to smaller estimates of τ_R. We found that for shorter recordings, estimates of R_tot were roughly the same as for 90 minutes, but estimates of τ_R were considerably smaller (S2 and S3 Figs).

To enable a meaningful comparison of the information timescale between neurons, one thus has to ensure that recordings are sufficiently long (in our experience at least 10 minutes), otherwise differences in τ_R may not be well resolved. Below 10 minutes, we found that estimates of τ_R could be less than half of the value that was estimated for 90 minutes, and also estimates of R_tot showed a notable decrease. In addition, all recordings should have comparable length to prevent that differences in history dependence or timescale are due to different recording lengths.