; ;

Single-trial Parameter Fits

This package contains MATLAB functions for two purposes:

  1. Fitting informative parameters to reduce time-series to two components
  2. Exporting the fit components along with fit information and relevant grouping metadata in order to model the fit parameters, highlighting group differences during recovery that occur on the scale of single-channel multi-unit activity.

Overview

This section indicates the fit parameters and generic database structure.

Motivation

The main fixed effect of interest is the interaction of Group, Electrode Location, and Post-Operative Day in predicting changes in the following response parameters:

Database Structure

These data are structured as summarized in the Data tab as described in Figures 1-3 here. Specifically, data in Figure 3 represent filtered views of rows from the main database table, T, wherein each row corresponds to the spike rate time-series estimate for a single recording channel, in alignment to a single behavioral event, from a single animal during one recording. In MATLAB, T contains the following variables (columns, listing most-relevant):

Parameter Fitting Procedure

Syntax

After recovering the master database, T, we recover the parameter fits using:

Where min_n_trials is an adjustable parameter that has been set to 5 for the exported table, meaning that there must be at least 5 (successful) trials from a given element of BlockID in order for that recording to be considered for further analysis. This is to avoid cases where we have a single trial for a given channel on some post-operative day, which could really bias the data if we have spurious fluctuation or a noisy channel for whatever reason on that particular instance. Ideally the restriction number could be higher, but given the nature of these data (the rehabilitation testing is restricted to a limited number of attempts to begin with), it's a compromise between causing spurious model associations and throwing away too much of the experiment.

Description

The bulk of the get_fitted_table function is in the parameterized built-in fmincon function provided by the Statistics and Machine Learning toolkit. The algorithm is as follows:

  1. Set a "hyperparameter" grid, sweeping a combination of relevant parameters.
  2. Evaluate at each point on the grid.
  3. Initialize the constrained optimizer with the "best" combination (in this case, determined by the minimum residual sum of squares).

The parameters are used by the gauspuls function (Signal Processing Toolbox) to fit a Gaussian-modulated cosine to the observed (smoothed) spike rate of each individual trial (data table row). Prior to the fmincon procedure, the individual trial spike rate is detrended using the following algorithm:

  1. Compute the trial median spike rate (b)

  2. Recover the least-squares optimal coefficient (a) of time (t) as a linear predictor of spike rate (Rate) as

  3. The de-trended rates are therefore the residual rate not explained by the model

Often, spike rate histograms in alignment to a motor event look like a series of peaks and valleys, indicating alternating periods where a spike is more likely (due to excitation related somehow to the behavior) or less likely (due to inhibition related to the behavior, or refractory period of having been recently likely to depolarize). This "shape" is captured well by gauspuls, which allows us to recover the following parameters:

tau

An offset parameter, which describes the time-shift of a zero-centered Gaussian-modulated cosine estimated at the relative times for which the spike rate was obtained.

omega

Frequency of the cosine component, which ranges on the order of a physiologically plausible behavior.

sigma

Fractional bandwidth of omega for the Gaussian pulse, which describes the bandwidth at half-maximum amplitude of the Gaussian "envelope."

Figures

Figure 1: Example of poor Gaussian pulse fit.

One of the most-common parameter-based exclusions (due to lack of model fit) was when the Envelope Bandwidth converged to a value where the product of sigma and omega exceeded 10-Hz, a value that is greater than the Nyquist Frequency given the bin width for generation of the initial spike rate histograms (60-ms). In this case, the fit (blue) captures a little bit of the modulation of the observed rate trajectory (black), but qualitatively doesn't capture the trajectory well, so such fits were excluded.

Figure 2: Distribution of fit errors.

After all exclusion filters are applied to the fit parameters table, the errors are an approximately overdispersed Poisson distribution with mean of 2.1 and variance of 6.8.

Figure 3: Example of a trajectory that is well-fit by the Gaussian pulse model.

In this example, the sum of square errors is less than 1. The reconstructed trajectory (blue) capturing a majority of the observed data (black).

Figure 4: Example of a trajectory that isn't perfectly fit by the Gaussian pulse model, but still acceptable.

The reconstructed trajectory (blue) still captures the qualitative components of the observed data (black) that we were hoping to model, even with an error term that is approximately 3 times the mean error value.

Figure 5: Even relatively large error terms still capture the desired data aspects.

Toward the tail of the error distribution of the retained observations, the model (blue) still captures the data well (black).

Figure 6: PCA summary of fit model parameters.

(Top): Cumulative percent of data explained by each principal component (individual stem lines). The first three components capture approximately 90 percent of the parameter space. (Bottom): Axes labels for this biplot correspond in color to the first three principal components from the top axes. 100 random observations are plotted as a scatter, showing that the majority of variance is explained by a component that is primarily weighted by the Peak Offset term tau while the third component may be most-heavily influenced by the frequency-related parameters (sigma and omega). The constant and coefficient for linear de-trending of individual trial rate series are included in the PCA, here labeled as Premotor Bias and Linear Time Bias, respectively.

Exported Table (statistics)

The table used to fit a Generalized Linear Mixed Effects model (Matlab function fitglme) is recovered from the table of estimated parameters using

Orthogonality of Estimated Parameters

This step initially performed a PCA to obtain components that qualitatively were orthogonal elements that when plotted using the built-in biplot function to view the PC coefficients turned out to correspond to parameters that related to the timing and frequency parameters independently. Instead of fitting the principal component scores, since the timing and frequency parameter estimates do not demonstrate much covariance I decided to fit the recovered parameters directly.

Exclusion Filters

The remove_excluded function also performs a series of filters that were included to reject poorly fit trajectories. Specifically, any time the product of sigma and omega is > 10 Hz (which is above the Nyquist frequency given the spike rate bin size), the data is excluded. Gr only contains rows in alignment to Grasp events, where the pellet was Successfully retrieved. Visual inspection of fitted data indicates that the gauspuls fits are robust across the majority of the sum-of-square error range; however, there was an outlier cluster that had very bad fits and so any rows for which the error exceeded 50 (a.u.) the trials were excluded. Finally, any time that the peak offset parameter tau exceeded 750-ms or was lower than -1000-ms relative to the alignment event, the trial was excluded, as values outside that range become exceedingly likely to be associated with some other behavior that is not under consideration in this experiment. After all exclusions, the number of observations was reduced from 211,424 rows to 20,220 rows. Categorical exclusions (total: 168,051) were performed first and breakdown (independently) as follows:

After categorical exclusions, removals from the remaining 43,373 observations due to parameter-fit rules (total: 23,152) breakdown (independently) as follows:

The remaining 20,221 rows were included in the following table (Gr from the code above) for statistical analysis contains metadata described in T as well as the following endpoint (response) variables:

In addition, the following covariates are included: