Variance Based Measure for Optimization of Parametric Realignment Algorithms.

Neuronal responses to sensory stimuli or neuronal responses related to behaviour are often extracted by averaging neuronal activity over large number of experimental trials. Such trial-averaging is carried out to reduce noise and to diminish the influence of other signals unrelated to the corresponding stimulus or behaviour. However, if the recorded neuronal responses are jittered in time with respect to the corresponding stimulus or behaviour, averaging over trials may distort the estimation of the underlying neuronal response. Temporal jitter between single trial neural responses can be partially or completely removed using realignment algorithms. Here, we present a measure, named difference of time-averaged variance (dTAV), which can be used to evaluate the performance of a realignment algorithm without knowing the internal triggers of neural responses. Using simulated data, we show that using dTAV to optimize the parameter values for an established parametric realignment algorithm improved its efficacy and, therefore, reduced the jitter of neuronal responses. By removing the jitter more effectively and, therefore, enabling more accurate estimation of neuronal responses, dTAV can improve analysis and interpretation of the neural responses.

Introduction

Many neurophysiological studies are investigating neuronal responses to external events. These studies range from simple stimulus evoked neuronal responses in the corresponding primary sensory areas, e.g. neuronal responses to light flashes in the primary visual cortex [1], to neuronal activity correlated to complex behaviours, e.g. neuronal correlates of abstract problem solving [2, 3]. In such studies, neuronal responses are usually extracted by averaging the neuronal signal in order to reduce the “noise”, i.e. parts of the neuronal signal that are not correlated to the stimulus or behaviour that is being investigated. This procedure relies on the assumption that neuronal responses are time locked to the corresponding stimulus or behaviour. This assumption can be challenged, however, as neuronal responses show temporal variability in relation to the corresponding stimulus or behaviour [4-8]. Depending on the amount of the temporal jitter, the underlying neuronal response estimated by averaging may be distorted (Fig 1), possibly leading to mistakes in subsequent analyses and incorrect conclusions.

Fig 1

Effect of single trial jitter on the estimation of the underlying neuronal response.

A: Neuronal responses are related to the external event (in this example stimulus; S), but are triggered (R) by an internal process, which is not precisely time-locked to the onset of the event. B: A certain amount of noise is recorded together with the relevant neuronal responses. C: During the experiment, the external event occurs multiple times, while the neuronal activity is recorded. D: When neuronal responses are aligned on the response start, the trial average response (F: blue line) is a good approximation of the real neuronal response. However, the response onset is unknown. The trial-averaged response aligned on the event onset triggers (F: red line) does not correctly reproduce the real neuronal response. In addition, the standard deviation across trials calculated using the event onset triggers (F: blue and red shaded tubes) is an incorrect estimate of the variability of neuronal responses. This example was generated using Gaussian white noise with a standard deviation equal to 5% of the maximum response amplitude (SNR = 20). Differences between the response starts and stimulus onsets were modelled using a Gaussian distribution with a standard deviation equal to 1/7 of the response standard deviation (σ = 700ms, σ = 100ms).

Several attempts have been made to realign the jittered neural responses using realignment algorithms [9, 10]. By aligning the single trial neural responses, the realignment algorithms can be used to accurately estimate the shape and the variability of the neural response [11], estimate response latencies [12, 13] and enable processing techniques that require single trial responses to be aligned [14].

Typically, parametric realignment algorithms offer a selection of different parameter values, potentially providing greater flexibility in term of neural responses for which the algorithm can be applied to reduce the jitter. However, selecting parameter values suitable for a particular neural response without having accurate knowledge of the response itself can be difficult. Specifically, in situations where the neural recordings are noisy and jittered, accurately estimating the shape of the neural response may not be possible.

Here, we demonstrate that the reduction of variability across trials (difference of time-averaged variance; dTAV) can be used as a measure of jitter reduction. This property is crucial for the operation of a parametric realignment algorithm as it can be used to identify a set of parameter values for which the algorithm performed an accurate realignment. We demonstrate this property by comparing the jitter reduction obtained by a parametric algorithm for a range of parameter values versus the jitter reduction obtained using parameter values selected by dTAV from the same range. This analysis was performed using an established parametric realignment algorithm [9, 10], referred to as MaxCorr in the rest of the text.

Methods

The Methods are presented in the following order. First, we describe our simple model of neuronal responses to an external event (stimulus or behaviour). Second, we present the motivation behind the design of the dTAV measure using analytic tools. Since the exact triggers of neuronal responses are not known in a real-world application of the algorithm, it is necessary to design such a measure of jitter reduction in order to optimize the parameters of parametric realignment algorithms. Third, we give a brief description of the previously published MaxCorr realignment algorithm [9, 10]. Finally, we describe the details of simulated data used to assess the performance improvements obtained using dTAV to optimize the parameters of the MaxCorr realignment algorithm.

A measure of jitter reduction

We assume that the neuronal signal is a superposition of neuronal responses r(t) evoked at response onset times t plus the Gaussian white noise signal η: where σ is the standard deviation of the Gaussian white noise process and N(μ,σ) is a normal distribution with a mean of μ and a standard deviation of σ. After recording the neuronal signal and if internal neuronal response triggers are known, one can estimate the neuronal response by calculating the response-triggered average : where M is the number of responses used to calculate the average; η(t) is the noise in i-th trial; and is the random variable drawn from a Gaussian distribution that follows from the presence of noise. In the following derivations, we will use operator E( ) for expectation, V( ) for variance, for sample variance and for sample covariance. The sample variance of the neural response r(t) in the presence of the noise is given by: and is distributed as a χ2 distribution [15]:

For a large number of degrees of freedom, i.e. M being large, the χ2 distribution can be approximated by a normal distribution using the following transformation [16]:

Using this transformation, Eq 5 becomes:

Experiments are usually designed in such a way that the studied neuronal responses are elicited by events, e.g. sensory stimuli or certain types of behaviour. However, the onset times of the neuronal responses are not known, since these are triggered internally by the brain. The time shift between an event t and the onset of the neuronal response t may not be constant and can be regarded as a stochastic process. In our model, the difference between these two time points is modelled by a Gaussian distribution: where μ and σ are the mean and the standard deviation of the distribution. If one has access to event times only, as it is the case in a real experiment where the time points t when the brain triggers a response are unknown, one can estimate the neuronal response by calculating the event-triggered sample mean : where is the average signal in the presence of the jitter but no noise (σ = 0). The sample variance of is given by:

The first term of Eq 10 is the sample variance in the absence of jitter and depends only on the noise η(t) in a way given by Eq 7. The second term is the variance contribution arising from the jitter in the absence of noise and depends only on the jitter. The third term is the sample covariance between the signal and the noise. We can rewrite Eq 10 as:

For normally distributed and small jitters t—t, we can use a Taylor series expansion to express :

For large M, we can use transformation (6) and obtain:

Since r(t + t—t) and the noise η are not correlated, the expectation of is zero while its variance depends on the shape of the neuronal response and, thus, cannot be precisely estimated.

Using Eqs 7 and 12, the expectation of for large M and small jitters can be expressed as:

Our intention was to design a quantitative measure of jitter that could be calculated without knowing the internal triggers of neural responses. Eq 14 provides a direct link between the jitter standard deviation σ and the sample variance of the event-triggered response which does not require internal neural response triggers. Specifically, Eq 14 shows how comparing values may be used to infer which one of any two σ values, e.g. σ’ and σ”, is larger than the other. If σ’ > σ”, the corresponding expectations of the , and , will follow this relationship:

To optimize the alignment, it is necessary to reduce the amount of jitter without the knowledge of the real neuronal response triggers t. A possible approach minimizes the variance of the stimulus triggered response (Eq 10) at one particular time point since a reduction of jitter may result in a decrease of and, hence, lead to a smaller variance (Eq 15). However, since the and terms depend on the noise in the signal and, thus, come from a stochastic process, their values might go up by chance and, therefore, mask the reduction of . Neighbouring time points of the neuronal response are correlated in time and so is the variance term arising from the jitter. On the other hand, the noise may be correlated on a smaller time scale. Therefore, the stochastic values of and may average out across time if the variance is averaged across a sufficiently large time window. A more reliable indicator of jitter reduction may, therefore, be the decrease of the time-averaged variance, TAV:

As discussed before, averaging over time can reduce the variance of the TAV, thereby increasing the reliability of the measure. On the other hand, integrating over periods of time where the neuronal response is small compared to the noise or completely absent may not increase the reliability of TAV but instead increase the variance of TAV. This trade-off means that the integration window across which the variance is averaged should neither be too long nor too short.

In this study, we used the difference of TAV, dTAV, as a measure of jitter reduction:

We used dTAV to optimize parameters of our realignment algorithm by assuming the following is true:

If Eq 18 is correct, we can choose those parameters of our realignment algorithm which lead to the strongest decrease of the time-averaged variability, i.e. we select the parameters for which dTAV is largest. While Eqs 15 and 16 indicate that Eq 19 is correct, analytical derivation of such relationship will depend on the neural responses and may not always hold. In the next section, we use numerical simulations to show that such relationship holds for two simulated neural responses that have a mono-phasic and a bi-phasic shape.

Numerical analysis of the reliability of dTAV as a measure of jitter reduction

To demonstrate the conditions for which the dTAV is a reliable measure of the reduction in jitter, we performed a set of simulations, each composed of 2000 repetitions of an experiment composed of 100 trials. The neuronal responses were simulated as mono-phasic and bi-phasic functions composed from Gaussian functions (Fig 2): where σ was taken to be 100ms. The shifts of the neural responses, t—t, were drawn from a Gaussian distribution with zero mean and standard deviation σ:

Fig 2

Simulated mono-phasic (left) and bi-phasic (right) neuronal responses.

Noise was modelled as white Gaussian noise with zero mean and standard deviation σ (Eq 1). Our simulations used discrete time with a time step of 1ms. The TAV was calculated for each combination of σ, σ and the integration time T; and for each simulation run k using the following equation:

Each simulation was performed by selecting a combination of σ, σ and T values. The used σ values ranged from 0ms to 120ms in steps of 1ms and the simulation for σ = 60ms was performed twice because a dataset with 60ms jitter was used as the starting point for the simulated realignment. T values ranged from 30ms to 990ms in steps of 30ms. The σ values were selected to model different signal to noise ratios (SNRs), defined as the ratio of the maximum absolute value of the neuronal response and the standard deviation of the noise σ. We used σ values that yielded SNR values of 0.03, 0.05, 0.08, 0.13, 0.20, 0.32, 0.50, 0.79, 1.26 and 2.00 for both mono and bi-phasic responses.

We used these simulations to emulate an experiment where the jitter standard deviation of the dataset was σ’ = 60ms before the realignment. This initial dataset was compared to datasets with jitter standard deviations σ” ranging from 0ms (no jitter) to 120ms (doubled jitter) which represented the dataset after the realignment. dTAV was then calculated for each combination of k, σ”, σ and T.

Ranges of dTAV values varied across different orders of magnitude for different σ and T values. We therefore normalized dTAV values by dividing them by the maximum of the absolute value of the 1st and 99th percentile. where is the X-th percentile operator acting over variables a and b; and max is the maximum value operator. The normalized dTAV (ndTAV) was binned in 50 equally wide bins spanning the space from -1 to 1. Binned values were used to calculate the probability of jitter reduction, , for different ndTAV values, while keeping σ and T constant. To show how the reliability of dTAV as a measure of jitter reduction depends on SNR and T, we calculated the contours in the space spanned by ndTAV and T for each value of SNR separately. We also calculated the joint probabilities for each combination of σ and T values, , in order to verify that the relationship in Eq 18 holds.

MaxCorr realignment algorithm

We tested whether dTAV-based optimization can improve realignment results obtained using previously published MaxCorr algorithm [9, 10]. The MaxCorr algorithm works by approximately maximizing crosscorrelations between each pair of trials in three steps. First, for N trials, N(N-1)/2 crosscorrelations CX for all possible trial pairs i and j and time lags (λ- λ) up to half of the trial length are calculated. Second, a parabolic function F(λ- λ) is fitted to the crosscorrelation between the i-th and the j-th trial around the time lag (λ- λ) for which the crosscorrelation is at maximum.

Parabolic functions F(λ- λ) are fitted using a neighbourhood (λ- λ) of several tens of milliseconds. In the last step, all parabolic functions are summed up to derive a new function F(λ2,…, λ), which is quadratic in all of its variables and, therefore, has a unique global maximum (λ2,…, λ) . The values of (λ2,…, λ) are calculated by solving the system of linear equations obtained by applying partial derivatives to the parabolic function F. (λ2,…, λ) are used to realign the trials relative to the first trial. MaxCorr algorithm was applied to the same signal used to evaluate the realignment of the dTAV algorithm and we evaluated jitter reduction metrics (σ’—σ”) / σ’ for both algorithms for comparison.

The implementation of the MaxCorr algorithm in the FIND Matlab toolbox [9] provides the user with the following options and parameters:

Maximum allowed correlation time lag (Δλ). Cross-correletaion is calculated only up to the provided value of the parameter thereby enforcing that all (λ- λ) values stay within the range of the parameter.

Processing of the cross-correlation coefficients. The user can choose whether or not to apply a natural logarithm to the cross-correlation coefficients before fitting the parabolic functions F(λ- λ).

Normalization of cross-correlation coefficients. The user can chose the type of cross-correlation coefficient normalization between: (i) no normalization, which leads to a linear drop of coefficients for larger time lags; (ii) unbiased normalization, which corrects for the linear drop of coeffiecints; and (iii) normalization to autocorrelation, which divides the coefficients with a squared root of a product of autocorrelations of both signals.

Iterative jitter reduction. The user can chose to repeat the whole MaxCorr realignment process multiple times (N), each time reducing Δλ by half.

Different settings of these parameters and options can improve the jitter reduction of a MaxCorr algorithm when applied to a specific dataset. However, a principal way of determining which parameter settings should be used for a given dataset has not yet been investigated. We tested whether dTAV can be used to select parameter values that improve the efficacy of the MaxCorr algorithm. I.e., we selected a set of parameters/options, applied the MaxCorr algorithm and calculated dTAV and jitter reduction (σ’—σ”) / σ’. After all permutations of the selected parameter values and option settings were used, we selected the (σ’—σ”) / σ’ for which the dTAV was the highest.

Simulated data

We used previously published MaxCorr realignment algorithm [9, 10] to investigate whether the dTAV–based optimization of parameters can improve its efficacy when performing realignment of simulated single trial neural responses.

We performed 16 types of simulated experiments that differed according to: (i) the number of trials in each experiment between 20, 50, 100 and 200 trials; (ii) the type of the neural response between mono-phasic (r) and bi-phasic response (r); and (iii) the distribution of the temporal jitters between a Gaussian distribution with zero mean and 0.1s standard deviation and an uniform distribution spanning -0.2s and 0.2s. For each type of simulated experiments, we performed 100 simulations. In each trial, the single channel neuronal response to an arbitrary stimulus was recorded at 1KHz.

Mono-phasic (r) and bi-phasic neural responses (r) were modelled to resemble reported neurophysiological responses as follows (Fig 2):

Neural responses were simulated in the following way:

Noise in the recordings was simulated as additive Gaussian noise with zero mean and different standard deviations σ: 31.62, 19.95, 12.59, 7.94, 5.01, 3.16, 2.00, 1.26, 0.79 and 0.50, which corresponded to SNRs of 0.03, 0.05, 0.08, 0.13, 0.20, 0.32, 0.50, 0.79, 1.26 and 2.00. Temporal distances between the stimulus times were drawn from a Gaussian distribution with a mean of 10s and a standard deviation of 10s. To keep neuronal responses from overlapping, temporal distances below 3s were redrawn. To avoid occasional very large jitters, all jitters with an absolute value above 0.3s were redrawn. The range of the uniform distribution was chosen to roughly contain 95% of the jitters of the Gaussian distribution.

We optimized parameters and options over the following range of parameter values / range of options: (i) Δλ = 50ms, 100ms, 200ms, 400ms and 800ms; (ii) processing of the cross-correlation coefficients: “lin” (no processing) and “log” (natural logarithm); (iii) normalization of cross-correlation coefficients: “none”, “unbiased” and “coeff” (normalization to autocorrelation); (iv) N = 1 and 3. dTAV was calculated using an integration time T of 350ms.

To correctly simulate the outcome of the experiment, we assumed that the person analysing the data would filter the data using a low-pass filter, given that low-frequencies dominate the simulated neuronal responses. We filtered the simulated recordings using 2nd order symmetric Savitzky-Golay filters [17, 18] with different time windows of 100ms, 250ms, 500ms or 1000ms. The time window length of the filter was therefore implemented as an additional parameter optimized in parallel with the MaxCorr algorithm parameters.

Results

dTAV as a measure of jitter reduction

Using the model of neuronal responses described in section 2.1, we calculated the dependence of dTAV on the reduction of the jitter for different SNR levels and integration times T by simulating 100-trial experiments 2000 times. Results are summarized in Fig 3 for the mono-phasic signal and in Fig 4 for the bi-phasic signal.

Fig 3

Reliability of dTAV as a measure of jitter reduction for mono-phasic neuronal responses.

Data obtained by simulating 2000 100-trial experiments. A: Expectation of dTAV as a function of the reduction of jitter standard deviation. Lines drawn only for values of SNR of 0.2 and higher. For lower SNR, 2000 repetitions were insufficient to provide a reliable estimate of the expected value of dTAV due to the high noise level. For the shown SNR range, the expected value of dTAV is independent of the SNR. B: The standard deviation (std) of dTAV as a function of the amount of jitter reduction for different SNRs. Standard deviations of dTAV for SNR of 0.2 and lower are above 10−6 and are, therefore, not shown. C: Probability of jitter reduction as a function of ndTAV for different SNRs. Panels A, B and C are shown for integration time T of 300ms. D: Values of jitter reduction and integration times for which the probability of correct dTAV prediction reaches 90%. For jitter reductions and integration times above the line, the probability for correct dTAV prediction, , is above 90%. For SNRs of 0.13 and lower, the probability of correct dTAV prediction never reached 90%.

Fig 4

Reliability of dTAV as a measure of jitter reduction for bi-phasic neuronal response.

See caption of Fig 3 for details.

For an integration time of T = 300ms, selected for presentation in Figs 3A–3C and 4A–4C, the expectation of the dTAV increased monotonously with the amount of jitter reduction and showed little or no dependence on SNR (Figs 3A and 4A). On the other hand, the standard deviation of dTAV (Figs 3B and 4B) depended strongly on the SNR and was comparable or larger than the expected value of the dTAV. For high SNRs, the expectation of dTAV surpassed the standard deviation of dTAV even for small reduction in jitter (see Figs 3A, 3B, 4A and 4B). In such cases, dTAV is a reliable measure of the amount of jitter reduction, even when the amount of jitter reduction is small.

To measure how well we can rely on the dTAV as a measure of jitter reduction, we calculated the probability of correctly predicting that the jitter was reduced based on the normalized dTAV values (Figs 3C and 4C). As the SNR was increased, the probability increased up to 1, even for the smallest jitter reductions. For low SNR values the probability never reached 1, even when the jitter was completely removed. To reach a substantial increase of the probability above 0.5, an SNR of about 0.20 or higher was needed.

To provide an insight into the dependence of the probability of correct dTAV prediction on the integration time, we calculated the values of integration times and reductions of jitter standard deviation for which the probability of correct dTAV prediction reached 90% (Figs 3D and 4D). For high SNRs the performance of the dTAV prediction was nearly independent on the integration time whereas for low SNRs the integration time had a stronger influence on the performance of the dTAV prediction. The performance of the dTAV prediction decreased faster for integration times below the optimal integration time, while it decreased more slowly for integration times bigger than the optimal integration time. Therefore, choosing a short integration time could be more disadvantageous than choosing a longer integration time.

So far we have computed the probability of jitter reduction given a reduction in dTAV. Next, we investigated whether the difference of dTAV is also predictor of the expected amount of jitter reduction, i.e. does a larger reduction in dTAV also indicate a larger amount of jitter reduction. If so, we can use dTAV to optimize parameters of our re-alignment algorithm by selecting parameter values that gave the highest dTAV values. To assess whether this is the case, we calculated the joint probability distribution of jitter reduction and dTAV for different SNR values for the mono-phasic neural response (Fig 5). For SNRs of 0.13 and lower, dTAV provided no or only very little information of the amount of jitter reduction but as the SNR increased (to values of about 0.2 and higher), the relation between dTAV and jitter reduction became less variable and dTAV became an increasingly good predictor of the amount of jitter reduction. For high SNRs (1.26 and higher), even small differences in dTAV indicated increased jitter reduction with high certainty. This suggests that dTAV is a good predictor of the amount of jitter reduction for sufficiently high SNRs (of about 0.2 and higher) and can be used to optimize the parameters of our realignment algorithm in such cases. The parameter selection based on dTAV will improve with increasing SNR.

Fig 5

Joint probability distribution of jitter reduction and ndTAV for different SNR values for the mono-phasic neuronal response.

An integration time window of T = 300ms was used. For low SNR, dTAV is uninformative as a measure of jitter reduction. As the SNR increases, dTAV becomes more informative of the jitter reduction, i.e the ability to differentiate different levels of jitter reduction based on dTAV improves substantially.

Re-alignment of simulated data using dTAV optimized MaxCorr algorithm

We used MaxCorr to realign mono-phasic or bi-phasic neuronal responses either for Gaussian or uniform distribution of jitters in 100 simulated experiments containing 20, 50, 100 or 200 trials for different levels of noise (Fig 6). MaxCorr parameters were determined using dTAV, which was calculated using an integration window starting at T = 0s and ending at T = 1s, both in respect to the stimulus times t. We intentionally did not want to use the results of the dTAV reliability analysis (Figs 3D and 4D) since a user of this method would not have the access to such results for his particular neural response. However, we did use the general finding of the dTAV reliability analysis that using an integration window wider than the response diminished the dTAV reliability less than using a window that is narrower than the response. In other words, we used a window that would certainly be wider than the optimal window, knowing that this does only weakly influence the reliability of dTAV as a measure of jitter reduction.

Fig 6

Jitter reduction obtained using dTAV-optimized MaxCorr.

Jitter reduction is shown for monophasic (left panels) and biphasic (right panels) simulated neural responses and for jitter distributed according to Gaussian (top panels) and uniform distributions (bottom panels). Different lines show jitter reduction for simulated experiments containing different numbers of trials. MaxCorr algorithm options and parameters were optimized over the following values: maximum allowed correlation time lag Δλ: 50ms, 100ms, 200ms, 400ms and 800ms; processing of the cross-correlation coefficients: “linear” and “logarithmic”; normalization of cross-correlation coefficients: “none”, “coefficient” and “unbiased”; number of consecutive iterations of the MaxCorr algorithm: 1 and 3; and the width of the filter window: 100ms, 250ms, 500ms and 1000ms. σ’ and σ”−jitter standard deviation before and after realignment, respectively. All results are averaged over 100 simulation repetitions; error bars depict the standard errors of the mean. For SNR of 1.3 and 2, standard errors are too small to be noticed on the plots.

For low SNRs (mono-phasic signal: SNR of 0.20 and lower; bi-phasic signal: SNR of 0.32 and lower) MaxCorr increased the amount of jitter, rather than decreasing it, except in 200-trial simulations with mono-phasic responses and SNR of 0.20, where the jitter was slightly reduced. For intermediate SNRs (mono-phasic signal: SNR of 0.32 and 0.50; bi-phasic signal: SNR of 0.50) MaxCorr reduced the jitter of responses, with its efficacy improving with the number of trials used for realignment. For high SNRs (SNR of 0.79 and higher), MaxCorr removed more than 83% of jitter from the recorded signal, regardless of the number of trials used, distribution of the jitter and the type of response used in the simulations.

To demonstrate that the dTAV optimization improves the efficacy of MaxCorr algorithm, we compared the jitter reduction obtained by dTAV-optimized MaxCorr with the maximum and the median of jitter reductions obtained using MaxCorr with all parameter sets used for dTAV-optimization, as specified in section 2.3. dTAV will contribute to the selection of parameter values (i) only if the jitter reduction is higher than the median jitter reduction, which is in some way analogous to random selection of parameter values, and (ii) only if the selected parameters can be used to obtain a large portion of the maximum jitter reduction value.

For all SNRs for which dTAV optimization of MaxCorr’s options and parameters reduced the jitter, jitter reduction was higher than the median jitter reduction (Fig 7). The proportion of the maximum jitter reduction recovered by the dTAV optimized parameters increased with the number of trials in the simulation and with the SNR. For SNR of 0.50 and higher, dTAV optimization of MaxCorr’s options and parameters recovered more than 85% of the maximum jitter reduction obtained by MaxCorr using the selected set of parameters, except for a bi-phasic response and 20-trial simulations, where the recovery was still above 50% (Gaussian distribution: 52%; uniform distribution: 58%). For SNR of 0.32 and a monophasic response, the recovery was 42% and 50% for 20-trial simulations for Gaussian and uniform distributions, respectively, and then increased with the increasing number of trials in the experiment up to 84% and 85% for 200-trial simulations for Gaussian and uniform distributions. For a biphasic response and SNR of 0.32 and for lower SNRs for both response types, the recovery is either low or there is even an increase of the jitter standard deviation after the realignment.

Fig 7

Maximum, median and dTAV optimized jitter reduction obtained by MaxCorr algorithm.

Main panels show maximum (black lines) and median (red lines) of jitter reduction obtained by using all permutations of options and parameter values used for dTAV optimization (see Fig 6 for the list of options and parameter values). Blue lines show reduction of jitter standard deviation obtained using dTAV optimized parameter values. Insets show the percentage of maximum jitter reduction recovered when using dTAV optimization (grey line with black dots); and the percentage of jitter reduction values, obtained when using all option selections and parameter values, that are lower than the jitter reduction obtained using dTAV optimization (purple line with black dots). Both are shown only for SNRs for which the jitter reduction obtained using dTAV optimization was positive. All results are averaged over 100 simulation repetitions; error bars depict the standard errors of the mean. For some SNRs standard errors may be too small to be noticed on the plots.

To further show that dTAV optimization selects efficient sets of parameter values, we counted the percentage of parameter value combinations for which the jitter reduction was lower than that obtained by using dTAV optimization, and called it dTAV percentile. For all SNRs for which dTAV optimization of MaxCorr’s options and parameters reduced the jitter, dTAV percentile was higher than 62% (Fig 7 insets). dTAV percentiles increase with increasing number of trials (mean dTAV percentile: 20 trials: 73%; 50 trials: 79%; 100 trials: 84%; 200 trials: 89%), confirming the property that more trials improve the robustness of the dTAV optimization. In addition, dTAV percentiles for SNR of 0.50 and higher were higher for biphasic neural response (minimum dTAV percentile: monophasic: 67%, biphasic: 77%), implying that dTAV optimization improves realignment to a larger extent for a more complex neural responses.

Discussion

We presented and evaluated a measure for the amount of temporal misalignment between the set of single trial neural responses, called difference of time-averaged variance (dTAV), that can be calculated without knowing the internal triggers of neural responses. We derived the relationship between dTAV and the precise amount of misalignment, as measured by the standard deviation of the temporal jitter between the trials. Using simulations, we showed that dTAV can be used as a measure of misalignment for a wide range of signal-to-noise ratios (SNRs). In addition, we showed that dTAV can be used to optimize the parameters of an established MaxCorr realignment algorithm designed to increase the temporal alignment of single-trial neural responses. MaxCorr realignment algorithm was applied to simulated sets of mono and bi-phasic neural responses. We showed that dTAV optimization of parameter values of the MaxCorr realignment algorithm leads to substantial efficacy improvements for noise levels up to two and three times larger than the neuronal response amplitudes for mono and biphasic neural response, respectively.

We used mono and bi-phasic potentials composed out of one and two Gaussian functions as examples of neuronal responses. In our simulations, the noise was modelled as Gaussian white noise independent of the simulated neural responses, and the amplitude and shape of the neural responses did not vary across trials. The success of the measure in selecting parameter values that effectively reduce the jitter for intermediate and low noise levels (SNR of 0.32 and higher) and for both mono and bi-phasic neuronal responses, suggests that the dTAV measure can be successfully applied in a large number of cases. At the same time, our results show differences between mono and bi-phasic response shapes, suggesting that the signal shape can affect the efficacy of the dTAV parameter optimization and, therefore, the performance of the realignment algorithms that use dTAV-based parameter optimization.

Our simulations assumed that the recorded neuronal signal is a continuously modulated signal, which is valid for example for local field potentials, electro-corticographic signals (ECoG), electro-/magnetoencephalographic recordings (EEG/MEG) as well as for optical calcium imaging, near-infrared spectroscopy (NIRS) and functional magnetic resonance imaging (fMRI) signals. Furthermore, our measure can also be applied to spike trains by estimating instantaneous neuronal firing from the spike times [19, 20] and using the instantaneous neuronal firing rates as the neuronal signal. In addition, the measure can be applied to continuously modulated signals which were derived from the aforementioned neuronal signals, such as time-resolved spectral amplitudes (e.g. extracted using short-time Fourier transform) or crosscorrelation measures. Misalignment of recordings that depend on other variables than time, such as space or frequency, can also be measured using dTAV by exchanging the time variable by the corresponding variable (e.g. spatial coordinate or frequency), thus enabling the effective use of parametric realignment algorithms.

In our study, dTAV was defined for one channel of neural signals only. Nonetheless, the measure can easily be extended to multiple channels by summing up dTAV of individual channels. MaxCorr algorithm has already been implemented for multiple channels [9, 10]. Therefore, using dTAV for MaxCorr parameter optimization in the case of multi-channel neural recordings can be easily achieved.

We show that the performance of MaxCorr algorithm can be improved by testing different parameter sets and calculating dTAV. This can be achieved by defining admissible values for each parameter (e.g. within an interval) and scan all possible combinations. If the number of admissible values for each parameter is very high, this process becomes computationally demanding and, therefore, may be time consuming. In our simulations we obtained good realignments within a reasonable timeframe (less than an hour), by optimizing the parameters across a limited number of values. Indeed, our results show (e.g. Figs 3 and 4) that the fine optimization of parameters may not be necessary and the required computational time for parameter optimization may therefore not be problematic.

In general, realignment algorithms can be used to improve the analysis of neurophysiological experiments by improving the estimation of neuronal responses. The obvious case is the estimation of the neuronal response by calculating trial averages. As shown in our example (Fig 1), even if neuronal responses are constant across trials and the noise is uncorrelated to the signal, averaging the jittered single-trial responses can lead to a distorted estimation of the response and the incorrect estimation of the noise. Removing the jitter can improve the estimation of the neuronal response and improves the estimation of the noise. A large number of neuroscience studies investigate neuronal responses related to sensory stimuli. When neuronal responses are well locked to the stimulus [21] realignment methods might be of limited use. On the other hand, neuronal responses may not be locked to the stimulus but, in addition to the stimulus, may also be affected by the internal neuronal state [22, 23] and, therefore be temporally jittered relative to the stimulus onset. Realignment algorithms can be used to find out whether the responses are locked to such internal events and compute the approximate timing of these events. Furthermore, neuronal responses related to behaviour may also be jittered. For example, neuronal responses related to movement planning [24, 25] may be jittered with respect to the times when the movements were initiated if the movements were either self-paced or triggered by another stimulus. Realignment algorithms could be used to align noisy individual trials in order to improve the accuracy of determining the underlying neuronal response (Fig 1). These more accurately determined neuronal responses would then also facilitate a comparison of responses between studies, allowing for identification of neuronal response parts shared between classes of responses. Computing the realignment times may give us also insights into the timing of internal events and into the variability of timings in internal cognitive processes.

Additionally, brain-machine interface systems that detect events based on neuronal recordings [26-33] may benefit from re-alignment algorithms. Such systems require a certain number of trials containing the neuronal responses to calibrate a template that is then used to detect the events from continuous neuronal recordings by matching the current neural activity to the template. If the jitter or the neuronal responses used to calibrate the template is reduced, the template will more closely resemble neural responses to events and, therefore, the detection may perform better.

In summary, we showed that the dTAV measure of misalignment can be used to improve the performance of realignment algorithms when applied on simulated neuronal responses for response waveforms commonly observed in neurophysiological recordings and noise levels higher than the neuronal response itself. Hence, the application of the dTAV measure and the realignment algorithms that use it to optimize their parameter values can improve analysis and interpretation of neuronal responses and improve the performance of asynchronous detection of events from neuronal recordings.