A Stochastic Dynamic Operator Framework That Improves the Precision of Analysis and Prediction Relative to the Classical Spike-Triggered Average Method, Extending the Toolkit

Data preprocessing

Off-line cluster cutting of both interneuron and motor unit spikes was performed using Plexon Neurotechnology Research System's Offline Sorter software. The first three principal components of 1.73-ms-long waveforms, centered on each recorded spike, were clustered using an automated, expectation–maximization sorting algorithm (Figueiredo and Jain, 2002; Shoham et al., 2003) followed by manual curation (Fig. 2B). Sorted spike trains and EMG signals were imported into MATLAB for further analysis (Fig. 2C). EMG signals were further processed off-line using a series of zero-phase filters: (1) a 60 Hz notch filter, (2) a 10 Hz high-pass filter (fourth order Butterworth), and (3) a 20-point root mean square (RMS) filter, which smoothed and rectified the EMG signal. We assigned EMG signal “state” based on signal amplitude. Since the smoothed and rectified EMG signal amplitude distribution fell into a roughly exponential distribution, we quantized EMG amplitude into signal states using logarithmic binning to allocate probabilities of state closer to a normal distribution: We defined 20 signal states for each EMG channel, representing equal intervals of the log-transformed EMG signal amplitude (Fig. 2D).

Theoretical consideration of when the spike-triggered average may fail

Biological systems are dynamic and incompletely observed. Consequently, signals measured within the system are often stochastic and could be spuriously correlated from an outside perspective. When determining a neurons’ contribution or sensitivity to a measured covariate signal, various assumptions are made regarding the structure of the relationship. The STA supports different causal models. It provides one method to reverse-correlate a neuron's spikes to preceding changes in signal amplitude. The STA also provides the means to predict future signal state based on a spike occurrence. The latter is the focus of the methods compared here. When spiking events (or at least the result of neural spikes) are considered independent and the stochastic component of the signal is uncorrelated with a spike, averaging signal amplitude around spiking events provides a straightforward method of isolating a single “response” of the signal corresponding to spiking events.

The time series data used in the STA comprises signal amplitude, x, indexed at time t, i.e., x(t) (generally, t is collected at the limit of digital signal sampling in the recording system; Fig. 3A). To mechanistically compute the STA, the pre-spike and post-spike signal levels are sampled over some finite interval, Δt, relative to spike time, s. For every time point relative to a spike, ti∈(s−Δt,s+Δt) , the concurrent signal amplitude, x(t_i), is averaged over all spiking events to find the mean, x̄(t_i) (Fig. 3C). These signal waveforms may be variously preprocessed prior to this averaging operation (e.g., mean-leveling) to further isolate spike-associated correlations. Performing this operation on all relative time points thus generates a mean waveform, x̄(t), which is effectively the signal impulse response to the spiking impulse.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

A, Electrophysiological recordings presented as time series data: Time-varying amplitude values collected over a duration. The relationship of a point process, such as a spike train, and a time series is often evaluated as the spike-triggered average (STA). Here, intervals of the continuous signal are sampled around spiking events, with samples averaged together. C, Under the STA framework, if a spike–signal relationship is significant, an “impulse response” may be present, showing a difference in the average signal behavior in the pre-spike versus post-spike averaged interval. B, A complementary approach is to treat the electrophysiological signal as a discrete random variable. By applying a quantization scheme, continuous signal amplitudes may be converted into a discrete time series of finite signal “states.” This permits signal behavior over an interval to then be described as a probability distribution of signal state. D, As with the classical spike-triggered average (C), the relationship between a spike and signal may be discerned from intervals of the discrete signal before and after spike. The analogous spike-triggered average probability distributions describe signal state prior to, and after, a spike. The correlated effect of a spike event on signal behavior can thus be described as the shift in the spike-triggered pre-spike and post-spike signal distributions.

In this framework, the neuron's contribution to the correlated signal (or the signal's effective stimuli, if testing how signal facilitates a spike in reverse correlation) is the convolution of the spiking impulses with the expected impulse response (i.e., the STA). Here, spikes are implicitly considered as independent and equivalent events with the spike rate reflecting neural excitation (i.e., the neuron behaves within the assumptions of the linear–nonlinear Poisson model; Schwartz et al., 2006). If these assumptions are upheld, then the STA should provide the maximum likelihood estimation of the spike–signal correlation (Paninski, 2004).

These assumptions bear consideration. Is a neural spike equally potent in predicting future signal state across all conditions? An initial exploration of this question may be directed at the covariate signal itself: Does the spike–signal correlation depend on signal amplitude? It is important to recognize that the operation for computing the STA enforces a type of dimensionality reduction: For each position in time relative to spike, t, the ensemble of observed signal responses forms a distribution of amplitudes, X(t). By definition, the classical STA merely utilizes the mean of this distribution, a single value, when reporting the “average” response at t. Further, this distribution is often assumed to be Gaussian (Normal) for statistical significance. If the STA is calculated from the change in signal amplitudes (i.e., Δx instead of x), then any consideration of amplitude-dependent spike effect correlations is also lost. In estimating a single response, the STA thus contains less complete information about signal behavior at time t compared with the full distribution of signal values around the spike.

The spike-triggered average distribution of signal states

We here propose a method of utilizing the complete spike-triggered state probability distribution to describe spike-triggered effects. Our process is as follows: We first apply an amplitude binning scheme (i.e., quantization) to our recorded time series data (e.g., EMG), such that each observation of the recorded signal amplitude may be assigned to a single discrete amplitude bin, x_i. [The number of such bins can be chosen for the signal based on signal range and precision required (Fig. 3A,B).] Here we define signal “state” as a binned range of signal amplitudes, resulting from the quantization scheme applied to the original, continuous, time series data. (Hence, “higher states” correspond to higher levels of original signal amplitude.) Here, the frequency of observing state over the entire time series is described by the distribution P(x) (Fig. 3B). Next, we infer the spike-triggered effect, classically inferred from the change in the amplitude (waveform) before and after a spike as the change in the distributions of signal amplitude around each spike (Fig. 3D). As with the classical STA, signal is sampled over a short interval (s − Δt, s + Δt) relative to the spike event, s (Fig. 4A). By quantizing signal into discrete states, potentially non-normally distributed signal amplitudes in the STA (Fig. 4B) can be described explicitly for any time point, or interval of time, as a distribution of state (Fig. 4C). We define the pre-spike distribution as the distribution of discrete signal states over the entire sampled interval leading to a spike p[x,(s−Δt,s)] , and the post-spike distribution as the distribution of states in the entire chosen sample interval after a spike, p[x,(s,s+Δt)] (Fig. 4D). For brevity, we will represent the discrete pre-spike and post-spike distributions as p(x0) and p(x1) , respectively. Thus, the effect of a spike, classically captured in the STA as the change of mean signal amplitude in time, may be equivalently described as the change of the probability distributions relative to a spike, Δp(x) (Eq. 1).p(x1)−p(x0)=Δp(x).

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

The spike-triggered average (STA) of a signal may vary based on background signal behavior. In this example, spiking events from a single motor unit recorded in the vastus externus (VE) muscle were compared against the aggregate EMG signal recorded from the biceps femoris (BI). A window 10 ms before spike to 10 ms after spike was used for identifying spike-triggered effects. A, The simple spike-triggered average of BI EMG at spike (red line) within this window has wide variability. The mean STA signal is given as a thick purple line, while ±1 standard deviation of the windowed signal is described by the breadth of the colored envelope. B, Separating the aggregate STA by whether the signal is above or below the mean aggregate STA (A) at time of spike provides two STAs (STA-A, STA-B) which exhibit different behavior. C, Signal amplitudes were converted to discrete signal states using a logarithmic binning scheme. By taking a histogram of binned signal amplitudes (“states”) at an arbitrary single time point (here, 5 ms after spike; position indicated by dashed blue vertical lines in A and B), this bimodality corresponding to STA-A and STA-B is also seen in state space. D, If signal amplitude over the entire perispike interval is quantized into amplitude states, the signal amplitude in the pre-spike and post-spike intervals may be represented as a pre-spike and post-spike distribution of state. E, When pre-spike and post-spike distributions are described using a joint probability distribution, the theoretical peak value corresponds to the maximum likelihood estimate (MLE) of the post-spike state given pre-spike state. Here, in theoretical data, the STA effect is the migration of the peak from the matrix diagonal (i.e., the average change in state). F, If, however, the distribution of pre-spike or post-spike signal distributions are multi-peaked, as in this physiological example (A–D), transitions to the global peak (X_B; STA-B) may be suboptimal to describe transitions around local peaks (X_A; STA-A). In this example, the best prediction of post-spike state is contingent on whether pre-spike state is closer to 8 or 12.

Practically, because the signal has been quantized into a finite number of states, the pre-spike, post-spike, and difference probability distributions are each captured in finite length vectors. The average pre-spike, post-spike, and difference distributions vectors could then be averaged across spiking events, as with the STA.

In this formulation, the maximum likelihood estimation of signal state around a spike in quantized state space is given as the global maximum of the joint distribution of the pre-spike and post-spike state distributions (Fig. 4E). By this description, the STA corresponds to the movement of the global mean off the joint state diagonal (i.e., the difference in the average pre-spike and post-spike state), which assumes a stationary and unimodal outcome. If, however, p(x0) or p(x1) are not unimodal (as in this example taken from physiological data), local maxima may exist in the joint distribution, independent of the global maximum (Fig. 4F). In this case, clearly, consideration of the pre-spike state can better inform predictions of post-spike state distributions than the average of Δp(x) alone. In these instances, the classical STA may perform as a suboptimal description for spike–signal prediction.

Stochastic dynamic operator framework

To capture state-dependent relationships in spike–signal correlations more precisely, measures of the changes from pre-spike to post-spike signal distributions should also be state sensitive. Specifically, we desire to know p(x1)=∑xp(x1|x0)p(x0) . That is, we seek the post-spike distribution, given the conditional probability of signal state after a spike and the prior pre-spike signal state distribution.

To map the state-dependent shift in pre-spike signal state distribution, we define a linear operator, L. L is implemented as a matrix which uses the pre-spike state distribution p(x0) to generate a column vector providing the change in the state distribution, Δp(x) , given an observed pre-spike probability distribution (Eq. 2).Δp(x)=Lp(x0). The relationship between the pre-spike and post-spike distributions can be described via the simple update p(x1)=p(x0)+Δp(x) . By combining and rearranging Equations 1 and 2, the post-spike state distribution may be predicted from a given pre-spike state distribution by operator L (Eq. 3).p(x1)=Lp(x0)+p(x0). Next, we need to derive L. Theoretically, L may be recovered as a standard linear estimate across all observed pre-spike to post-spike distributions. Here, we further elaborate on this process. Ultimately, Equation 14 is the main result of this derivation. For a set of k spikes, we may define the sampled ensemble of pre-spike distributions as the vector concatenation P(x0)=[p(x0)1…p(x0)k] , the post-spike distributions as P(x1)=[p(x1)1…p(x1)k] , and the change in state distributions asΔP(x)=[p(Δx)1…p(Δx)k] . For n total states and k spiking events, this produces n-by-k matrices. L is then estimated as follows:ΔP(x)≈LP(x0). Postmultiplying by the transpose of the pre-spike distribution establishes the following:ΔP(x)P(x0)T=LP(x0)P(x0)T. This could, in principle, be simplified to a standard linear regression problem as follows:Rxy=LRxx, where Rxy=ΔP(x)P(x0)T and Rxx=P(x0)P(x0)T and are square matrices of order n. L could thus be derived as follows:L=RxyRxx−1. However, the linear estimation of the SDO derived in Equation 7 is often ill-posed because the covariance matrix of pre-spike state distribution (R_xx in Eq. 7) may be noninvertible. Because this matrix often could not be defined in our datasets via Equation 7, we developed an alternative linear estimation method, which we will describe below. To introduce that method, we must consider how L behaves. In practice, L is a square matrix whose order equals the number of signal state levels used for analysis. It provides the averaged spike-triggered change in probability of state, given pre-spike state, as estimated from the observed data (Eq. 4). Thus, L is effectively a kind of state-dependent spike-triggered average, albeit for state distributions rather than waveforms. Because L operates on probability distributions to describe stochastic signal dynamic changes, L is termed as SDO (see Sanger, 2011 for development and definitions). SDOs were originally developed to describe and manipulate continuous-time stochastic dynamics governed by differential equations. Here, we use the properties of the discrete time SDO to create a difference equation which describes how the spike-triggered average change in probability of state depends on the prior state.

The matrix operator, L, is not uniquely defined. The SDO matrix represents a “probability flow” from the pre-spike to the post-spike distribution. Probability flows out from each state j, to other states i, i≠j . The total outward and inward flow component from each state must sum to 0 to conserve the total probability, and states cannot increase their own probability by flowing into themselves. Four conditions are necessary and sufficient for an SDO matrix to behave as a linear operator (Sanger, 2011): (1) The matrix diagonal must contain nonpositive elements, (2) off-diagonal elements must be non-negative, (3) each column must sum to 0 (if the SDO is modeled as a left-stochastic matrix, as used here), and (4) the sum of positive elements within a column must be no greater than 1. From the perspective of a single initial state, x₀ = j, these constraints mean that the change of probability of maintaining state at x₁ must be nonpositive (Δp(x1=j|x0=j)≤0) , while the change of probability of transitioning to any other subsequent state must be non-negative; Δp(x1≠j|x0=j)≥0 (i.e., one may only increase, or sustain, the probability of transitioning to another state and only decrease, or sustain, the probability of maintaining state). A valid estimation of the SDO must thus uphold these constraints.

SDO matrix linear estimation method

Having identified how the SDO matrix should behave, here we describe an alternative equation which efficiently generates compliant SDO matrices regardless of input state probability distributions and bypasses the invertibility problem faced by the classical least-squares method in Equations 6 and 7. In Equation 2, the SDO is a matrix that provides the change in the state probability, Δp(x)=Lp(x0) while upholding the listed constraints. If R_xx cannot be directly inverted to derive L via Equation 7, let us instead consider when R_xx is not necessary. When the pre-spike state distribution is binary, a compliant SDO matrix (L) can always be estimated from the outer product of the change in state distribution and the binary pre-spike state distribution (Eq. 8). This corresponds to R_xy given in Equation 6.L=Δp(x)p(x0)T. In Equation 1, this can alternatively be represented as the difference of the initial distributions:L=[p(x1)−p(x0)]p(x0)T. However, when p(x₀) is not binary and |Δp(x)|>0 (as is usually the case), Equations 8 and 9 produce a matrix with positive elements occurring on the diagonal and negative elements elsewhere than the diagonal. Hence while the vector Δp(x) can be predicted from the resultant transformation matrix (i.e., upholds Eq. 2), it is not an SDO. However, a compliant SDO can be generated by the sum of SDO matrices generated for each of n discrete states in the pre-spike distribution (Eq. 10). (This is equivalent to iteratively and independently populating each column of the SDO in Eq. 9.)L=∑j=1n[p(x1)−p(x0=j)]p(x0=j)T. Here, the probability of state j in the pre-spike state distribution is given as p(x0=j) , and this distribution contains zeros everywhere except at element j. (Alternatively, this can be thought of as the summation of n SDO matrices generated by Eq. 8 via binary pre-spike distributions, with the jth SDO weighted by the probability of observing state j in the pre-spike interval.) Because vector p(x0=j) contains a single nonzero element, the jth term only populates column j of the SDO matrix in the net summation, thus each term is effectively independent. This technically permits a special case where the outer product can be distributed to the minuend and subtrahend. Thus, Equation 9 may be rewritten as follows (Eq. 11):L=∑j=1np(x1)p(x0=j)T−p(x0=j)p(x0=j)T. The sum of the differences may be rewritten as the difference of two sums (Eq. 12):L=∑j=1np(x1)p(x0=j)T−∑j=1np(x0=j)p(x0=j)T. Here, because the jth term of the left sum of outer products again independently contributes to the jth column of the matrix, the summation is equivalent to the outer product of the pre-spike and post-spike distributions directly, p(x1)p(x0)T . Similarly, the jth term generated via the outer product operation of binary vectors in the right sum only operates on a single element on the array diagonal at L_jj. Hence, the right summation produces a diagonal matrix, whose diagonal contains the pre-spike distribution. By employing these two simplifications, we can then derive L as follows (Eq. 13):L=p(x1)p(x0)T−diag(p(x0)). This form utilizes the original pre-spike and post-spike state probability distributions directly without any further requirements. Finally, we average the SDO matrix across all spikes. For k spiking events, if the ensemble of pre-spike distributions is provided via the vector concatenation P(x0)=[p(x0)1…p(x0)k] , and the ensemble of post-spike distributions is the similar concatenation P(x1)=[p(x1)1…p(x1)k] , the SDO is computed as follows:L=1k(P(x1)P(x0)T−diag∑P(x0)). Equation 14 is the important final result that serves as the basis of the SDO analysis used below and produces compliant SDOs which uphold both Equations 2 and 3. It is direct and algorithmically straightforward once formulated this way. The prior development equations are important to understand how Equation 14 was derived, but Equation 14 is the main one for the reader less interested in the derivation methods.

A further technical point can be made. Due to the structural assumptions made in Equation 10, and consistent with Sanger, by this formulation of Equation 14, for each spike observed, every state in pre-spike distribution probabilistically flows to all states in the observed post-spike distribution. Consequently, when calculated for a single spike, the columns of the SDO matrix are substantially identical, save for the matrix diagonal, with the jth column scaled according to the probability of observing state j in the pre-spike distribution, p(xo=j) . State-dependent behavior of the final SDO matrix, averaged over these per-spike effects, is thus not derived from the sum of state-dependent transitions of pre-spike and post-spike distributions by individual spiking events, but rather emerges from differences in sampled pre-spike and post-spike distributions of data across the entire spike train.

Next, we elaborate SDO analysis as a flexible tool to describe stochastic behavior correlations within simulated and physiological data. The SDO methods and analyses described below, including the implementation of Equation 14 for SDO matrix estimation, were generated using the SDO Analysis Toolkit, a package we wrote and assembled. Analyses presented were generated in MATLAB 2023a, running on a Dell Precision Tower 5810 using an RTX3060 GPU with a Windows 10 operating system.

Validation of SDO methods in simulated data

We first tested our inference that the classical spike-triggered average method may be insufficient to consistently detect state-dependent or paradoxical interactions within the nervous system using computational models of different stochastic processes in silico. We generated data via well-characterized systems and assessed the statistical specificity and sensitivity of the STA and SDO methods to detect these relationships. Eight methods of generating stochastic time series data (Generators Y1–Y7), corresponding roughly to the seven hypotheses (above) and an additional ARIMA process (Y8), were used. To correlate spike-triggered effects, one putative spike train was drawn from random indices for each simulation, in common for all the time series methods (Generators Y1–Y8). If a simulated stochastic time series had a spike-triggered effect (Generators Y2, Y3, Y6, Y7), the associated spike effects were applied at these time indices. We describe these stochastic generator processes as follows:

1. [Y1] Low-pass filtered white noise: Pre-spike and post-spike distributions of state should not systematically vary. Spike has no effect (Fig. 5A).

2. [Y2] Low-pass filtered white noise with spike-triggered random effects: This control is similar to Y1, but at time of spike, there is an added higher amplitude random (white noise) response added during the 10 ms duration after spiking event. The effect of spike increases signal variance (Fig. 5B).

3. [Y3] Low-pass filtered white noise with consistent spike-triggered control effects: This is the same as Y1, but with an added consistent impulse response during the 10 ms following time of spike. The spike has a consistent effect (Fig. 5C).

4. [Y4] Low-pass filtered white noise with stabilizing dynamics independent of spikes: Here, the signal was iteratively synthesized using the update equation Y4[t]=Y4[t−1]+ΔY4[t] . To maintain consistency between generator models, for every time t, we utilized the differential of Y1 to derive the update to Y4 (ΔY4[t]) . This generator model was stabilized toward a fixed point using the dynamic equation: ΔY4[t]=k(x0−Y4[t])+ΔY1[t] . Here X₀ is a fixed point (0), and k is the stabilizing coefficient (1 > k > 0). Essentially, this equation biases future signal toward zero, with an effect magnitude proportional to its current distance from zero. The differential of the Y1 signal (ΔY1[t]) is thus used here as an added random noise term. Spike has no effect (Fig. 5D).

5. [Y5] Markov process independent of spikes: Here, a sequence was generated via a random walk of a Markov matrix. The Markov transition matrix was composed of a 100 × 100 identity matrix convoluted with a Gaussian kernel, biasing transitions to nearby states. The signal was mean-leveled around zero. Spike has no effect (Fig. 5E).

6. [Y6] Stabilizing dynamics with consistent spike-triggered effects: This is the same as Y4, with the spike impulse responses of Y3 added. Spike has a consistent effect (Fig. 5F).

7. [Y7] Stabilizing dynamics with spike-triggered dynamics: This is the same as Y4; however, in the 10 ms duration after the spike, a second stabilizing equation is additionally active (of the same form in Y4). For the spike-triggered dynamics, the equilibrium of the control is positive, X₀ > 0, such that background dynamics and spike-triggered dynamics stabilize toward different points. A spike activates a consistent change in dynamics (Fig. 5G).

8. [Y8] Autoregressive Integrated Moving Average (ARIMA) stationary stochastic process: ARIMA models are commonly used to parsimoniously describe and forecast stochastic time series data. Here, a stationary ARIMA(3,2) model was used. The signal was mean leveled around zero. Spike has no effect (Fig. 5H).

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

Comparison of the STA impulse response and SDO matrices obtained from simulated stochastic control data. Time series and point process data were simulated with different stochastic processes and spike–signal control relationships (or lack thereof). For each simulated spiking event, a sample of time series data was extracted in the window 10 ms before spike to 10 ms after a spike. The STA impulse response was estimated as the average signal amplitude across all spike-triggered samples. The pre-spike window was defined as 10 ms before spike until the time of spike, and the post-spike interval as the 10 ms following the spiking event. In each case, the SDO was generated using the distribution of states extracted across each interval. Twenty states were linearly defined from the minimum and maximum of the observed signal values. Time series were stationary about a mean value (usually zero), with similar amplitude ranges on either side of the mean. Hence, the average signal level corresponds to an intermediate signal state (i.e., 10–11), and SDOs had values around the middle of the matrix. Representative examples of STAs and SDOs from simulated data are plotted here. The 500 spike-triggered window signals are overlaid as gray traces in time. A, [Y1] Time series was generated using low-pass filtered white noise. Spikes have no effect. STA is flat and SDO is symmetrical about the diagonal. B, [Y2] Spikes evoke increased noise level impulses in otherwise white noise. STA is flat, while the SDO is elongated in the horizontal dimension (corresponding to the increased signal variance after spike). C, [Y3] Spikes have consistent effects within a white noise signal. STA extracts the impulse response. The SDO is elongated above the matrix diagonal, showing an increase in transition toward higher states (amplitudes). D, Signal is generated by a stabilizing dynamic equation. Spike has no effect and STA is flat. The SDO demonstrates more stabilizing effects than in white noise (Y1); above-average input states (11–15) are pushed toward lower states (5–10) while below-average input states (5–10) are pushed toward higher states (10–15), as indicated by SDO elements. E, [Y5] Markov process. Subsequent signal amplitude transitions to similar amplitude as prior amplitude. Spike has no effect. STA is flat. The SDO indicates that changes in signal distribution at spike time are minor and bidirectional (indicated by the limited band on the matrix diagonals). F, [Y6] Spike has consistent effect within a process described by a stabilizing dynamical equation. As with Y3, STA extracts the impulse response and SDO is elongated. G, [Y7] Spike induces dynamic effects, while signal is also generated by a stabilizing dynamic equation. STA is flat. SDO is elongated in the horizontal dimension, indicating pre-spike signal variance is less than post-spike variance. H, [Y8] Signal is generated from an ARIMA(3,2) process. STA is flat. As with the white noise (Y1), the SDO is relatively symmetrical about the matrix diagonal.

We compared the STA and SDO methods of analyzing data from these various generators’ simulations of model systems in which spikes were independent of dynamics or participated in the control of dynamics.

Computational implementation of the STA methods

We used three methods of calculating the spike-triggered average between a spike and (simulated) rectified EMG and one measure of spike tuning. A spike-triggered effect was considered “significant” if any test of effect returned a positive result using α = 0.05. To avoid spurious correlations, only spike trains with a minimum of 500 spikes were included for analysis.

Computational implementation of SDO methods

To generate, visualize, and analyze SDOs, we developed the SDO Analysis Toolkit, written in the MATLAB programming language. For each time series signal (Generators Y1–Y8 in simulated data, or rectified aggregate EMG in the physiological data), signal states were defined as 20 intervals between the channel-wise minimum and maximum values of signal amplitude observed over the entire dataset. (Note that simulated data used a linear scale while in the collected EMG data we used a logarithmic scale, based on observed signal amplitude distributions.) Under the SDO model, each spike was treated as an independent event. For each spike, we extracted signal states in a short interval (10 ms) before (“pre-spike”) or after (“post-spike”) the spike event. The frequency histogram of states within these observed intervals relative to spike comprised the pre-spike and post-spike probability distributions. The pre-spike distribution included the time bin containing the timestamp as the final element. SDO matrices were then generated for each combination of spiking unit (i.e., neuron) and signal. Note that as a first demonstration of the SDO methods explored here, we consistently used 20 states for quantizing signal amplitude, and pre-spike and post-spike durations of 10 ms each, although these parameters are among those which may be readily modified within the SDO Analysis Toolkit.

Linearly estimated SDOs were used for analysis

In addition to the linear estimation method described above (Eq. 14), we also included and tested a constrained optimization method utilizing the MATLAB Optimization Toolbox applied within the SDO Analysis Toolkit which currently uses the fmincon solver. Here, the SDO is sought as a matrix which minimizes the least mean-squared prediction error of the observed change of the pre-spike distribution, Δp(x), given the pre-spike distribution as in Equation 4, while also upholding SDO matrix constraints. However, we found that the SDOs generated by linear estimation had significantly more stable matrix element values than SDOs generated via constrained optimization, measured by element-wise distances between matrices, across different numbers of spikes (250, 500, 1,000, 2,000) within the same stochastic realizations (Table 1). We found this result across all types of simulated data (Table 2). This suggested that the optimization algorithm was potentially overfitting to the data or was following different exploitations of matrix redundancies within a dataset rather than converging onto a single and more generalizable solution. Indeed, prediction errors of the constraint-optimized SDOs were often higher than those derived from the linear-estimated SDOs generated from the same simulated data, indicating that the solver is not identifying the globally optimal SDO.

We infer that although the four SDO matrix constraints are sufficient to define an SDO, for constraint-optimized estimation there may be additional uncharacterized constraints, preferable optimization techniques, or else much longer datasets than tested here (with concomitantly slower computation) may be necessary to derive a stable and generalizable SDO. We mention the optimization method for completeness in describing the SDO Analysis Toolkit and include it as a method of deriving the SDO but caution that significant care will be needed for its use with physiological data. Since the direct linear equation (Eq. 14) seems to be the faster and perhaps better approach in our estimation, we used Equation 14 for deriving and characterizing the SDO throughout results, as described below.

Generating spike-independent SDOs for significance testing

Spike-triggered SDOs were tested for significant spike–signal correlations against null hypotheses generated using Monte Carlo sampling. We used the following process: First, we generated 1,000 shuffled-spike trains, for each observed spiking unit (e.g., neuron). We included two methods of shuffling spike times, both included within the SDO Analysis Toolkit. The first is to directly shuffle the observed interspike intervals on a trial-wise basis, thus preserving the durations between impulses. The second method is to estimate the rate process underlying the observed spike train, assuming spikes were generated by a renewal process, then to resample this process 1,000 times. (The SDO Analysis Toolkit allows for various filtering kernels for this purpose.) In both cases, the 1,000 “shuffled” spike trains from the tested spiking unit were then used to produce 1,000 “shuffled” spike train SDOs, using the same signal time series as used for the tested experimental SDO. Hence, in the null hypotheses we are testing if the changes in signal dynamics measured at a spike are a result of sampling (interspike intervals) or dynamics near a spike, but not directly associated with it (renewal process). The number of shuffles used for significance testing is a tunable parameter within the SDO Analysis Toolkit.

SDO matrix normalization

The joint probability of the pre-spike and post-spike state distributions describes the overall probability of observing state transitions in the pre-spike and post-spike intervals (Fig. 6A). The column sum of the joint distribution provides P(x₀), the average pre-spike state distribution from the ensemble, and the row sum of the joint distribution provides P(x₁), the average post-spike state distribution from the ensemble. As estimated in Equation 4, the SDO matrix (Fig. 6B) describes the change in state probability distributions across the ensemble and will reflect the probability distributions used to construct it. Specifically, Δp(x1=i,x0=j) depends on both the conditional change in probability given an initial state Δp(x1=i|x0=j) and the probability of observing the initial state, p(x0=j) , in the pre-spike distribution. Both the value of Δp(x) and confidence in this value are affected by the probability of observing it in the data. This latter quantity may differ between the spike-triggered (method H7), spike-shuffled (null), and background (all time points, method H4) SDOs.

• Download figure
• Open in new tab
• Download powerpoint

Figure 6.

SDOs alter background state transition matrices to compose spike-triggered transition matrices. In this example, pre-spike and post-spike distributions of signal state were generated from the spikes of a single motor unit in the vastus externus against the analog EMG signal recorded in the biceps femoris, using a time interval of 10 ms for both the pre-spike and post-spike distributions. A, The spike-triggered average joint distribution of the pre-spike and post-spike state distributions behaves as the description of state transitions around spike. B, A spike-triggered SDO captures the change of state probability across the 10 ms pre-spike and post-spike distribution of states, given an occurrence of spike. This SDO shows strong directional effects, evidenced by the asymmetrical banding relative to the diagonal. C, Normalization of each column of the joint distribution to 1 creates a left transition matrix, representing the expected probability of post-spike state distribution given a pre-spike state distribution. D, Normalizing the columns of the SDO by the same factor used in the joint distribution (C) creates the normalized SDO. This form is used to predict the conditional change of state distributions.

To account for differences in the sampling of pre-spike state (column scaling) between the spike-triggered, background, and “spike-shuffled” SDOs, the columns of the SDO matrices are first normalized to the conditional form, (Δp(x1|x0)) . When normalizing an SDO matrix, each column j of the SDO is scaled by the reciprocal of the probability of state j in the pre-spike distribution, 1p(x0=j|s) . [When applied to the joint distribution of state p(x1|x0) , this is equivalent to normalizing each column to 1 (Fig. 6C), producing a transition matrix.] In this normalized form, the SDO represents the conditional change of the probability of state, given an initial state distribution, Δp(x0|x0) , and hence emphasizes magnitude of correlated signal change, by state (Fig. 6D). The normalized (conditional) SDO matrix is the form that is used for generating predictions of Δp(x) as in Equation 2 for all H1–H7. However, it is not ideal to directly compare between normalized SDOs or use these to estimate significance here because this form may distort the contributions of very rarely sampled states to the measures of significance.

Instead of normalization, for comparison and testing significance, the conditional SDOs of the spike-triggered (method H7), “shuffled”, and background (method H4) SDOs are rescaled by the average pre-spike state distribution of the observed spike-triggered ensemble, P(x₀), simulating what the “shuffled” and background SDOs would resemble if they had been taken from the same pre-spike state distributions as the spike-triggered SDO (method H7). [Note that this process produces the same method H7 spike-triggered SDO as originally sampled (Fig. 6B).] These rescaled SDOs can then be used for evaluating the matrix similarity-based measures of significance (as mentioned above). Both the normalized and rescaled SDOs can be readily plotted and examined within the SDO Analysis Toolkit.

Identifying significant SDOs

We first identified significant SDOs by looking for differences in the spike-triggered (method H7) and spike-shuffled SDO matrices. [Differences between spike-triggered SDOs (method H7) and background SDOs (method H4) were also measured but were not used to directly determine significance of spike effects]. In the rescaled spike-triggered and shuffled-spike SDOs, we measured four features of SDO effect and one feature of state tuning:

1. Matrix element values: The sum of squares error (SSE) calculated between each element of the spike-triggered SDO matrix versus the distribution of values for that element in the shuffled distribution.

2. Matrix similarity: Overall cumulative SSE of the spike-triggered SDO matrix and the shuffled SDO matrices or between the respective joint distributions.

3. State-wise post-spike distribution directional shift/bias: For each input state (SDO column), the sum of the elements above the matrix diagonal versus the elements below the matrix diagonal. (This corresponds to the expected change in mean imposed by the SDO, conditional on input state.)

4. Array-wise total post-spike distribution directional shift/bias: The sum of all SDO matrix elements above the diagonal versus the elements below the diagonal (i.e., does the spike coarsely facilitate or inhibit signal?).

5. Pre-spike state tuning: The Kullback–Leibler divergence (KLD) between the observed state at spike, p(x|s) , versus shuffled-spike times. (Note that this is not a measure of the spike effect but can be used to identify significant relationships between signal level and a spike occurrence.)

The significance of the spike-triggered SDO (method H7) was determined by comparing the spike-triggered test statistic versus the null distribution of test statistics from the 1,000 shuffled-spike SDOs. Here, the difference between the spike-triggered SDO statistic and the mean of the shuffled-spike SDO statistics was compared with the internal variance of the population of 1,000 “shuffled” SDOs (i.e., is the test result from the spike-triggered SDO significantly different than would be predicted from the shuffled distribution?) A p value of <0.05, with Bonferroni’s correction for the number of states (here, 20), was used to determine significance. Within the SDO Analysis Toolkit, the number of shuffles and the alpha value for significance testing can be readily modified.

Visualization of SDO and SDO effect predictions

Visualization of the correlation between spikes and signal is important for intuition. The spike impulse response of the STA may be readily visualized and interpreted from the average time-varying waveform of sampled signal amplitude across all spiking events. In contrast, the SDO describes the change in state probability distributions without prescribing a specific temporal waveform structure. Thus, any time effects would not be immediately apparent from the pre-spike and post-spike distributions for the SDO. Nonetheless, considering the underlying state transition probability and signal behavior over time is important.

We observed that the spike-triggered response may vary by signal state at time of a spike (Fig. 7A). To visualize and identify such variations of signal behavior relative to spiking events in the SDO framework, we developed the spike-triggered impulse response probability distribution (STIRPD) description (Fig. 7B). The STIRPD is effectively a transformation of the classical STA into the probability state space but has the advantage of more robustly describing the relationship between the spiking source and signal behavior in time. To calculate the STIRPD, as with the STA, for each spiking event, s, the state-quantized signal is sampled on an interval around spike from (s−Δt,s+Δt) . Subsequently, the distribution of signal states at every time point t, p(x,t), is calculated from the same relative time bin across all sampled signals. The value of each element represents the probability of observing a given state (as opposed to another state) at a given time after a spike, STIRPDij=p(x=i|t=j) .

• Download figure
• Open in new tab
• Download powerpoint

Figure 7.

Expressing STAs and probability distributions. The spike-triggered impulse response probability distribution (STIRPD) is a combination of the classical spike-triggered average and state probability distributions. A, The relationship between the spike train of an interneuron is compared against EMG (filtered acausally). In the state-dependent STA, a set interval [s − Δt, s + Δt] of signal behavior is extracted for every spiking event (s), sampled in discrete time. A state, x, is assigned to every observed time point, t, (t∈[s−Δt,sΔt]) within this interval to provide a sequence of states, x(t) . When segregating the STA by state at time of spike, x(s) , state-dependent behavior may be observed. B, The spike-triggered impulse response probability distribution extends the STA by calculating the distribution of states at each measured time point p(x,t) across all spiking events. Because both states and time are discrete, the STIRPD is a finite matrix and may be displayed as an image. Here, the “arching” signal behavior present at the higher states in the STA is present as a dark band, but the STIRPD nonetheless captures the relatively invariant signal behavior in states 6–10 in the perispike interval. Horizontal separation of state mass in the STIRPD (i.e., banding) will result in multimodal pre-spike or post-spike distributions elements of the joint probability matrix and may be indicative of cases where the STA is insufficient. Here, the STA is overlaid on the STIRPD in blue.

Because both time and state are defined discretely, the STIRPD may be displayed as a raster image and interpreted qualitatively. The familiar spike-triggered average (in state space: mean state) is laid over the probability distribution. The probability distribution of state at time of a spike, p(x|s), is captured as a single column vector for the time bin containing the spike on the STIRPD and thus can reveal if the distribution is expected to be unimodal or multimodal or if there is a coarse signal behavior change around the spike. The summation of state probability over time in the STIRPD during the pre-spike interval, [s−Δt,s] provides p(x0) , while summation of the STIRPD over a given post-spike interval provides p(x1) , used in the calculation of a specific SDO. The STIRPD description thus links the SDO and STA visualization methods.

As with the STIRPD, we could also qualitatively infer the relationship between a target signal and the spike source neuron for the pre-state (x₀) to the post-spike (x₁) state change using visualizations of the SDO matrix directly (Fig. 8A). When inspecting the SDO, we also found it is intuitively helpful to plot a shear of the SDO matrix such that the prior matrix diagonal, corresponding to (i=j) , is aligned horizontally (Fig. 8B). This transformation emphasizes the direction of the SDO effects, Δp(Δx,t) : Elements above the shear SDO horizontal correspond to transitions toward higher states (x1>x0) , while below-diagonal corresponds to transition toward lower states (x1<x0) . The sign and magnitude of these elements may be used to interpret SDO effect for a given input state. [For example, for an input state which has positive elements (increased probability) above the shear-horizontal “toward higher states” and negative elements (decreased probability) below the shear-horizontal “toward lower states,” the net effect for that input state is a shift toward the higher states].

• Download figure
• Open in new tab
• Download powerpoint

Figure 8.

Alternative visualizations of the SDO matrix. The effect of a spike-triggered SDO on the probability of signal state between two observations x₀ and x₁ may be qualitatively inferred from the SDO directly. A, The SDO matrix may be divided into three regions, the primary diagonal (red line), representing the probability of maintaining state; the “above-diagonal” (orange) region, representing the probability of increasing state; and the “below-diagonal” (purple) region, representing the probability of decreasing state. B, The square SDO matrix may be sheared to orient the prior matrix diagonal to the horizontal. The shear SDO visualization emphasizes the change in probability of state, Δp(x), relative to the directional change of state, Δx. As SDO operators are sparse matrices and usually concentrated around the diagonal when time intervals are short, the shear SDO may be trimmed to a few rows to compactly display the effects of an SDO. C, The magnitude and directional effects of an SDO can be coarsely described by summing Δp(x₁|x₀) separately for the regions above (orange) and below (purple) the SDO matrix diagonal. In the quiver SDO visualization, these column sums are represented as vertical vectors oriented above or below the abscissa, respectively. The length of each vector for state x is the magnitude of the summed change in probability. The vector heading projecting from the abscissa corresponds to whether the change of state, Δx, is positive (orange) or negative (purple). The elements of the diagonal are plotted as a line within the same axes. The orientation and magnitude of vectors for each state, x, provide a coarse measure of influence of Δp(x₁|x₀).

A third visualization of the SDO we explored is the “Quiver” representation of an SDO (Fig. 8C), which emphasizes the magnitude and coarse directional biases of Δp(x) for a given input state. Here, in the quiver representation, each column of the shear SDO matrix is summed over (1) the above-horizontal elements and (2) the below-horizontal elements, and each represented by a vector whose magnitude is proportional to this sum, originating at the abscissa, and projecting up and down, respectively. The elements of the reference horizontal row (the original SDO matrix diagonal), corresponding to Δp(x1=x0) , are coplotted by a simple line (because the SDO matrix diagonal is defined to be nonpositive, this trace is negative). These vectors indicate the coarse directional effect of the SDO for a given pre-spike state: If both vectors above and below have an equivalent magnitude for a state, then the change in post-spike state overall is not directional (i.e., there is no change in the mean of the post-spike state distribution contributed by this input state), while unbalanced vectors indicate a coarse direction of change in the post-spike predicted state distribution.

Assessing the prediction utility of significant SDOs

Finally, in comparing methods, we tested the capability of the SDO (method H7) to predict post-spike signal behavior within our datasets. We assessed two types of prediction: (1) predictions of the post-spike state distributions, p^(x1) , and (2) predictions of most likely single state within the same post-spike interval, x^1=argmaxxi:p^(x1=xi) . Predictions of the post-spike state distribution (1) capture overall signal dynamics, while predictions of single post-spike states (2) reflect utilization of such a prediction toward directing an exclusive action or decision (e.g., “go”/”no-go”, tracking-target, etc.). (Note that even in instances where the spike-triggered SDOs did not significantly differ from shuffled-spike or background, SDOs could still be used to predict signal behavior, but we have limited our scope here to significant cases to compare the STA vs the SDO directly.)

When evaluating predictions of an SDO (i.e., the correlation between a particular spiking unit and signal), the observed pre-spike and post-spike distributions of signal state were extracted for each spiking event. The predicted post-spike state distribution, p^(x1|H) , were generated on an event-wise basis for each of the seven matrix model hypotheses (methods H1–H7) described above (including the spike-driven SDO as method H7), in accordance with Equation 14. The “single best state” of each predicted post-spike distribution was defined as the most likely state (i.e., peak). When evaluating model-predicted post-spike distributions, the experimentally observed post-spike distribution, p(x1) , served as ground truth. Hypotheses which minimized the reconstruction error of the post-spike distributions with the ground truth model were considered the model of best fit (argminH:p(x1)−p^(x1|H) for H1–H7). (Note that here, the H3 matrix is used to predict the effects of the classical STA represented in the probability space, corresponding to the average change in signal distributions which would otherwise be predicted from the time-varying amplitude.) Spiking events, and associated errors, were treated as independent observations. The difference between predicted and observed distributions was quantified using the KLD and log-likelihood. The difference between predicted and observed single states was quantified using the error frequency (e0), absolute error (e1), and squared (e2) error. Because the spike-wise error rate for these predictions to single states are integer values, and the population of errors is not expected to resemble a continuous distribution (and hence may be ill-suited for tests of means or medians), we instead calculate the cumulative error and bootstrap this statistic to evaluate the confidence intervals. If the distribution of two model hypothesis errors did not overlap at the 95% confidence interval, they were considered significantly different. Observed error rates were considered significantly different if there was no overlap of the estimated 95% confidence intervals, corresponding to a p value of <0.05. We used Cohen's D statistic (Cohen, 1988) to measure the effect size between bootstrapped cumulative error distributions. We hypothesized and expected that, by capturing information of the spike effect more effectively, the spike-triggered SDO (H7) predictions would be the best most of the time.