Analysis of Spike Counts as Proxy for Excitability

A common argument regarding neurophysiological activity changes during reorganization following ischemia is that deficits observed during post-stroke motor paresis relate to unbalanced "aggregate excitation" as evidenced by the total counts of unit spiking in layers with presumptive descending projections to the spinal cord. One hypothesis is that this imbalance is caused by the loss of inhibition contralateral to the injury that was previously supplied by the callosal fibers which degrade when cell bodies in the ischemic tissue are lost. In turn, this may cause hyper-excitability of the intact hemisphere, which strengthens inhibition of the injured hemisphere via its own callosal projections.
We did not test this hypothesis explicitly. Rather, we simply wanted to test the data that we collected against existing hypotheses in the field. One can imagine that many confounding covariates would be present in a study designed to test unit dynamics during motor behavior, with respect to a simple comparison of aggregate multi-unit rates between the two hemispheres. Therefore, this script is provided wherein the first part allows the user to manipulate the defined epoch window as well as the alignment behavior, to generate the statistical results as needed and actively compare the resulting output.
% "Run All" button (not shown in html export)

Load Data

clc;
clearvars -except R G
if exist('R','var')==0
fprintf(1,'Loading raw rates table...');
R = getfield(load(defaults.files('raw_rates_table_file'),'R'),'R');
fprintf(1,'complete\n');
end

Ad hoc model generator

This section allows the user of the live-script to generate a statistical model by inserting fixed timing epoch and alignment parameters as desired. If you are reading the html-exported version, or just want to skip to the results part, click here.
In the batch script analyze_raw_rates.m, we used some fixed periods around behavioral alignments, based on the general observation that task-related modulation of unit activity fluctuates on a timescale of a few hundred milliseconds near motor events of interest. In this case, we selected TASK_EPOC = [-450,150] for the task epoch in each of the three alignments, and TASK_EPOC = [-1350,-750] in alignment to Reach only as a "non-behavior" epoch. When modifying the values of TASK_EPOC below, the first value must be lower than the second value. Click Run Section when parameters are configured as desired.
mdl = "Spikes~Group*Area+(1|ChannelID:Day)+(1+Duration|AnimalID)+(1|SupportLimbMovement)+(1|Trial_ID)+(1|ICMS)"; % Wilkinson notation model
For statistics, we will fit a generalized linear mixed-effects model for a Poisson-distributed variable (spike counts in some epoch of interest) using a log link function. We will maximize the model log-likelihood using the Laplace fit method. In addition, any exclusion information is included with the glme variable in the output table G at the end of this script.
The model is defined using Wilkinson Notation, which includes main effects and interaction effects when terms are specified using the * operator. The only fixed terms in the model are:
The rest of the model includes two groupings of Random effects, which are actually pretty crucial in fitting the model well (the main thing being that each channel should have a random intercept term that has a nested grouping by day). There is also a grouped random effect for animal, which includes Duration of the total reach, the Extension duration, and Retraction duration, since each of these things influence how "close" temporally the behavioral events of interest are relative to one another through the execution of the behavior.
ALIGN = 'Grasp'; % Parameter-1: Event
TASK_EPOC = [-990 ,450]; % Parameter-2: Relative epoch (ms)
% "Run Section" button (not shown in html export)
adhoc_glme = analyze.stat.fit_spike_count_glme(R,ALIGN,'Successful',TASK_EPOC,'Model',mdl);
============================================== GLME: Grasp trials (Successful) -> For spikes occurring -990.0-ms to 450.0-ms relative to Grasp...complete Generalized linear mixed-effects model fit by ML Model information: Number of observations 29565 Fixed effects coefficients 4 Random effects coefficients 4525 Covariance parameters 7 Distribution Poisson Link Log FitMethod Laplace Formula: Spikes ~ 1 + Group*Area + (1 | ChannelID:Day) + (1 + Duration | AnimalID) + (1 | SupportLimbMovement) + (1 | Trial_ID) + (1 | ICMS) Model fit statistics: AIC BIC LogLikelihood Deviance 756253.679688356 756344.917499946 -378115.839844178 756231.679688356 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF pValue {'(Intercept)' } 4.34527019438396 0.0675546180242805 64.3223264591946 29561 0 {'Group_Ischemia' } -0.00410171631846071 0.0652532883321384 -0.0628583849687855 29561 0.949879698248617 {'Area_RFA' } -0.0748940444617274 0.0113569059017951 -6.59458175574823 29561 4.33640886524715e-11 {'Group_Ischemia:Area_RFA'} -0.000157556181335499 0.0107870827787706 -0.0146060046600899 29561 0.988346607306514 Lower Upper 4.21286015458847 4.47768023417945 -0.132001048111583 0.123797615474662 -0.0971540824338147 -0.05263400648964 -0.0213007156245716 0.0209856032619006 Random effects covariance parameters: Group: ChannelID:Day (3227 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.483299011830188 Group: AnimalID (10 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std' } 0.185482124302949 {'Duration' } {'(Intercept)'} {'corr'} -0.195847049361466 {'Duration' } {'Duration' } {'std' } 0.134132881791443 Group: SupportLimbMovement (2 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.0142011887470843 Group: Trial_ID (1271 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.149194569835617 Group: ICMS (5 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.0288144183760502 Group: Error Name Estimate {'sqrt(Dispersion)'} 1 R-squared: Ordinary: 0.463301344504417 Adjusted: 0.463246877606596 ==============================================
analyze.stat.panelized_residuals(adhoc_glme);
R-squared Ordinary: 0.463301344504417 Adjusted: 0.463246877606596
analyze.stat.surf_predicted(adhoc_glme,'Area','Group','Tag',sprintf('%s: [%d-ms to %d-ms]',ALIGN,TASK_EPOC(1),TASK_EPOC(2)));
Ad hoc model. This will change depending on what parameters are entered by the user of the live-script. Note that we are using counts of spikes instead of spikes normalized by duration, so depending on the length of time that is specified in TASK_EPOC, there will be a bias toward more or less spikes.

GLME model export for Test Statistic Endpoints

Values generated here were used in the corresponding Results discussion section of the associated manuscript.
% Use fixed pre-defined values as shown here
Align = {'Reach';'Grasp';'Grasp';'Support';'Reach'}; % Pellet-retrieval behavioral event alignment
Outcome = {'Successful';'Successful';'Unsuccessful';'Successful';'All'}; % Pellet-retrieval Outcome filter
tStart = [-450; -450; -450; -450; -1350]; % Relative epoch start times (ms)
tStop = [ 150; 150; 150; 150; -750]; % Relative epoch end times (ms)
glme = cell(size(Align));
% "Run Section" button (not shown in html export)

1. Export model for Successful Reach trial spike counts during behavior

glme{1} = analyze.stat.fit_spike_count_glme(R,Align{1},Outcome{1},[tStart(1) tStop(1)],'Model',mdl);
============================================== GLME: Reach trials (Successful) -> For spikes occurring -450.0-ms to 150.0-ms relative to Reach...complete Generalized linear mixed-effects model fit by ML Model information: Number of observations 35064 Fixed effects coefficients 4 Random effects coefficients 4674 Covariance parameters 7 Distribution Poisson Link Log FitMethod Laplace Formula: Spikes ~ 1 + Group*Area + (1 | ChannelID:Day) + (1 + Duration | AnimalID) + (1 | SupportLimbMovement) + (1 | Trial_ID) + (1 | ICMS) Model fit statistics: AIC BIC LogLikelihood Deviance 560045.658457738 560138.772690401 -280011.829228869 560023.658457738 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF pValue {'(Intercept)' } 3.63702951184947 0.0733101178955625 49.6115627181337 35060 0 {'Group_Ischemia' } 0.023317061349299 0.0717131560947712 0.325143427218358 35060 0.745074462935193 {'Area_RFA' } -0.0266545769548263 0.0116183382669815 -2.29418152082702 35060 0.0217859574588479 {'Group_Ischemia:Area_RFA'} -0.0139655315569984 0.0109752521340026 -1.27245655830827 35060 0.203219388336961 Lower Upper 3.49333936050745 3.78071966319149 -0.117242994319298 0.163877117017896 -0.049426887678752 -0.00388226623090074 -0.0354773731065963 0.00754630999259948 Random effects covariance parameters: Group: ChannelID:Day (3375 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.505443694835514 Group: AnimalID (10 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std' } 0.287131713974417 {'Duration' } {'(Intercept)'} {'corr'} -0.999999999744363 {'Duration' } {'Duration' } {'std' } 0.0857966184735757 Group: SupportLimbMovement (2 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 1.7456327177593e-05 Group: Trial_ID (1272 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.13783561172523 Group: ICMS (5 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.0311411796472692 Group: Error Name Estimate {'sqrt(Dispersion)'} 1 R-squared: Ordinary: 0.411485098553609 Adjusted: 0.411434740746868 ==============================================
analyze.stat.surf_predicted(glme{1},'Area','Group','Tag','Reach Onset','ZLim',[29.5,46]); % Keep constant z-limits
Model 1: multi-unit spike counts during Reach onset. There is a statistically significant effect of Area (p = 0.022; DF = 35060; tStat = -2.29), but neither fixed-effect involving Group has a significant impact in the model. In terms of the predicted model output for differences in per-group spiking, we can see that the model-predicted mean difference between the groups is relatively small (on the order of one spike).

2. Export model for Successful Grasp trial spike counts during behavior

% "Run Section" button (not shown in html export)
glme{2} = analyze.stat.fit_spike_count_glme(R,Align{2},Outcome{2},[tStart(2) tStop(2)],'Model',mdl);
============================================== GLME: Grasp trials (Successful) -> For spikes occurring -450.0-ms to 150.0-ms relative to Grasp...complete Generalized linear mixed-effects model fit by ML Model information: Number of observations 29514 Fixed effects coefficients 4 Random effects coefficients 4523 Covariance parameters 7 Distribution Poisson Link Log FitMethod Laplace Formula: Spikes ~ 1 + Group*Area + (1 | ChannelID:Day) + (1 + Duration | AnimalID) + (1 | SupportLimbMovement) + (1 | Trial_ID) + (1 | ICMS) Model fit statistics: AIC BIC LogLikelihood Deviance 471610.735223518 471701.954043584 -235794.367611759 471588.735223518 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF pValue {'(Intercept)' } 3.47635389806183 0.0791185806269178 43.9385270882767 29510 0 {'Group_Ischemia' } 0.0168880486850163 0.0767594953883083 0.220012502682353 29510 0.82586293939458 {'Area_RFA' } -0.0760182994916975 0.0120638401739863 -6.30133509689711 29510 2.9925599137014e-10 {'Group_Ischemia:Area_RFA'} 0.0302887344475609 0.0115717636376634 2.61746915992807 29510 0.00886296850884949 Lower Upper 3.321277969027 3.63142982709666 -0.133563968594181 0.167340065964214 -0.0996639615837718 -0.0523726373996233 0.00760756420331645 0.0529699046918053 Random effects covariance parameters: Group: ChannelID:Day (3225 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.508703194283732 Group: AnimalID (10 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std' } 0.219138456805056 {'Duration' } {'(Intercept)'} {'corr'} -0.251550168795761 {'Duration' } {'Duration' } {'std' } 0.221320625802525 Group: SupportLimbMovement (2 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.0200667381712226 Group: Trial_ID (1271 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.161945228968784 Group: ICMS (5 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.0253520152368285 Group: Error Name Estimate {'sqrt(Dispersion)'} 1 R-squared: Ordinary: 0.435660016146892 Adjusted: 0.435602645087876 ==============================================
analyze.stat.surf_predicted(glme{2},'Area','Group','Tag','Grasp Onset','ZLim',[29.5,46]);
Model 2: multi-unit spike counts during Grasp onset. The shape of the expected output manifold for the generalized linear mixed-effects model with the same structure as Model 1 changes somewhat as the temporal region-of-interest shifts from Reach to Grasp. In addition to a significant effect of Area (p << 0.001; DF = 29510; tStat = -6.30), there is now a significant interaction effect for Group*Area (p = 0.009; DF = 29510; tStat = 2.62). The model actually predicts more spikes in Ischemia*RFA than Ischemia*CFA; by comparison, the intact ranking is the opposite. This is interesting, because it goes against the Interhemispheric Inhibition hypothesis. An alternative is that there is still interhemispheric inhibition, but the reorganization of activity related to the "redistribution of function" necessary for the remaining cells that survive in the hemisphere with the lesion leads them to become more active for some reason around the Grasp. Another possibility is that there is some other feature of behavior present during this time in Intact rats, which causes a higher likelihood of CFA spiking in the hemisphere ipsilateral to the grasping forelimb, for example, movement of the supporting forelimb.

3. Export model for Unsuccessful Grasp trial spike counts during behavior

% "Run Section" button (not shown in html export)
glme{3} = analyze.stat.fit_spike_count_glme(R,Align{3},Outcome{3},[tStart(3) tStop(3)],'Model',mdl);
============================================== GLME: Grasp trials (Unsuccessful) -> For spikes occurring -450.0-ms to 150.0-ms relative to Grasp...complete Generalized linear mixed-effects model fit by ML Model information: Number of observations 38537 Fixed effects coefficients 4 Random effects coefficients 5825 Covariance parameters 7 Distribution Poisson Link Log FitMethod Laplace Formula: Spikes ~ 1 + Group*Area + (1 | ChannelID:Day) + (1 + Duration | AnimalID) + (1 | SupportLimbMovement) + (1 | Trial_ID) + (1 | ICMS) Model fit statistics: AIC BIC LogLikelihood Deviance 654803.864015983 654898.017131058 -327390.932007992 654781.864015983 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF pValue {'(Intercept)' } 3.55534086955833 0.0716644428786991 49.6109468900245 38533 0 {'Group_Ischemia' } -0.0416413012798485 0.0662979543477154 -0.628093305284365 38533 0.529946540844898 {'Area_RFA' } -0.0904766198034605 0.00778681102295259 -11.6192135056019 38533 3.71014599573949e-31 {'Group_Ischemia:Area_RFA'} 0.0167804902142262 0.00778679830331555 2.15499227803052 38533 0.031168613816828 Lower Upper 3.41487673041031 3.69580500870635 -0.171586985787422 0.0883043832277247 -0.10573896837008 -0.0752142712368405 0.00151816657841979 0.0320428138500325 Random effects covariance parameters: Group: ChannelID:Day (4008 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.46973335341334 Group: AnimalID (10 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std' } 0.291152705436787 {'Duration' } {'(Intercept)'} {'corr'} -0.756431378601529 {'Duration' } {'Duration' } {'std' } 0.196623109985754 Group: SupportLimbMovement (2 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.0366121590854065 Group: Trial_ID (1790 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.170278117018333 Group: ICMS (5 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 7.50878955813196e-05 Group: Error Name Estimate {'sqrt(Dispersion)'} 1 R-squared: Ordinary: 0.440582123413158 Adjusted: 0.440538569741506 ==============================================
analyze.stat.surf_predicted(glme{3},'Area','Group','Tag','Unsuccessful Grasps','ZLim',[29.5,46]);
Model 3: multi-unit spike counts during Unsuccessful grasp attempt onset. A similar pattern to Model 2 for the relative significance of activation by Area and Group is seen when only Unsuccessful Grasp trials are used for alignment. Notably, in Ischemia*RFA, there is now a reduced likelihood of spiking during the Grasp (specific to unsuccessful outcomes).

4. Export model for Successful Support trial spike counts during behavior

% "Run Section" button (not shown in html export)
glme{4} = analyze.stat.fit_spike_count_glme(R,Align{4},Outcome{4},[tStart(4) tStop(4)],'Model',mdl);
============================================== GLME: Support trials (Successful) -> For spikes occurring -450.0-ms to 150.0-ms relative to Support...complete Generalized linear mixed-effects model fit by ML Model information: Number of observations 19627 Fixed effects coefficients 4 Random effects coefficients 3406 Covariance parameters 7 Distribution Poisson Link Log FitMethod Laplace Formula: Spikes ~ 1 + Group*Area + (1 | ChannelID:Day) + (1 + Duration | AnimalID) + (1 | SupportLimbMovement) + (1 | Trial_ID) + (1 | ICMS) Model fit statistics: AIC BIC LogLikelihood Deviance 335843.160504743 335929.891780675 -167910.580252372 335821.160504743 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF pValue {'(Intercept)' } 3.47340581315441 0.0800237696288746 43.4046762513562 19623 0 {'Group_Ischemia' } 0.0706326076750676 0.0784074413472154 0.900840614888602 19623 0.367684125056784 {'Area_RFA' } -0.0884776280995819 0.0126929032681669 -6.97063754684704 19623 3.25524846159652e-12 {'Group_Ischemia:Area_RFA'} -0.0105135620806888 0.0120475339610416 -0.872673371553605 19623 0.382851813796652 Lower Upper 3.31655243192389 3.63025919438494 -0.0830526329225872 0.224317848272722 -0.113356795932805 -0.0635984602663585 -0.0341277512902548 0.0131006271288773 Random effects covariance parameters: Group: ChannelID:Day (2560 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.527521690285567 Group: AnimalID (10 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std' } 0.276090526591347 {'Duration' } {'(Intercept)'} {'corr'} -0.999999999984288 {'Duration' } {'Duration' } {'std' } 0.0656303169025966 Group: SupportLimbMovement (1 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 5.36182720889897e-05 Group: Trial_ID (820 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.145541629613485 Group: ICMS (5 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.0317934629566562 Group: Error Name Estimate {'sqrt(Dispersion)'} 1 R-squared: Ordinary: 0.370722085728685 Adjusted: 0.370625880574386 ==============================================
analyze.stat.surf_predicted(glme{4},'Area','Group','Tag','Support Onset','ZLim',[29.5,46]);
Model 4: multi-unit spike counts during Support onset. The model for Support onset is similar to the model for Reach (Model 1), and makes sense; we would expect that during movement of the supporting forelimb, M1 in the intact hemisphere contralateral to the supporting forelimb would be more active than premotor cortex in the hemisphere ipsilateral to it. When we align specifically to a supporting forelimb movement, coherence to that behavioral event is stronger than to the other elements of behavior that are presumably co-occurring. Therefore, we expect this activity to "outweight" any coherent averaging of activity related to the Reach or Grasp that might bias an increase in RFA expected counts during this epoch.

5. Export model for All Reach trial spike counts removed from behavior

% "Run Section" button (not shown in html export)
glme{5} = analyze.stat.fit_spike_count_glme(R,Align{5},Outcome{5},[tStart(5) tStop(5)],'Model',mdl);
============================================== GLME: Reach trials (All) -> For spikes occurring -1350.0-ms to -750.0-ms relative to Reach...complete Generalized linear mixed-effects model fit by ML Model information: Number of observations 80183 Fixed effects coefficients 4 Random effects coefficients 7743 Covariance parameters 7 Distribution Poisson Link Log FitMethod Laplace Formula: Spikes ~ 1 + Group*Area + (1 | ChannelID:Day) + (1 + Duration | AnimalID) + (1 | SupportLimbMovement) + (1 | Trial_ID) + (1 | ICMS) Model fit statistics: AIC BIC LogLikelihood Deviance 1516295.56039534 1516397.77313015 -758136.780197668 1516273.56039534 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF pValue {'(Intercept)' } 3.6010286817735 0.0573601583332636 62.7792667665151 80179 0 {'Group_Ischemia' } -0.0553734180058561 0.0542197955847619 -1.0212767755513 80179 0.307126407844185 {'Area_RFA' } -0.0424288769110821 0.00804742755639877 -5.27235276288328 80179 1.35033775014239e-07 {'Group_Ischemia:Area_RFA'} -0.041366955262837 0.00723686182720332 -5.71614551314755 80179 1.09352678561916e-08 Lower Upper 3.48860314014176 3.71345422340524 -0.16164386883626 0.0508970328245482 -0.0582017831941125 -0.0266559706280516 -0.0555511579268246 -0.0271827525988495 Random effects covariance parameters: Group: ChannelID:Day (4651 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.462566193956428 Group: AnimalID (10 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std' } 0.221151649340413 {'Duration' } {'(Intercept)'} {'corr'} -0.707441463807171 {'Duration' } {'Duration' } {'std' } 0.134754114384249 Group: SupportLimbMovement (2 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.00774331629956484 Group: Trial_ID (3065 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.177481032737534 Group: ICMS (5 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.0374179701331303 Group: Error Name Estimate {'sqrt(Dispersion)'} 1 R-squared: Ordinary: 0.342094377384128 Adjusted: 0.342069761002434 ==============================================
analyze.stat.surf_predicted(glme{5},'Area','Group','Tag','Pre-Reach | All-Outcomes','ZLim',[29.5,46]);
Model 5: multi-unit spike counts prior to Reach onset (all outcomes). This "control" model is potentially also very interesting. There are statistically significant fixed effects for both Area (p << 0.001; DF = 80179; tStat = -5.27) and the interaction Area*Group (P << 0.001; DF = 80179; tStat = -5.72), which indicate that on average, RFA has a reduction in the expected number of spikes during this epoch that is "removed" from the behavior, but with the partial-dependence plot indicating that this effect may be primarily driven by the interaction term, which sees a substantially lower number of expected spikes in RFA of Ischemia rats (35.06 expected spikes) by comparison to any of the other combinations (39.68 expected spikes in Ischemia*CFA; 40.26 expected spikes in Intact*RFA; and 41.66 expected spikes in Intact*CFA). The contrast between the relatively similar expected spike counts for both areas in the Intact group, compared to the large reduction in expected spikes in the Ischemia group during this epoch suggests that, at least during pre-movement epochs, the Interhemispheric Inhibition hypothesis remains viable. However, in general, during the execution and "read-out" of the motor behavior, the comparative increase of spiking in RFA of Ischemia rats suggests that Inhibition may play less of a roll during the actual behavioral execution itself.

General Comment

Notice that the shape of the output manifold for the recovered generalized linear mixed-effects models changes substantially depending upon what times we select relative to the behavioral alignment event. One reason that it is difficult to study this system purely in terms of unit excitation is because defining critical temporal epochs that are used for such analyses has a tremendous impact on the statistical model end-point. A dynamical systems perspective, in which the cortical areas under control are considered as part of a control system, allows a broader perspective in order to discern qualitative trends.

Fit full model

We can fit a model that incorporates each of the Successful (during-behavior) alignments, but due to the problems stated above, it is unlikely that the fit for this model (without including even more terms) would be very good.
full_glme = analyze.stat.fit_full_spike_count_glme(R);
============================================== GLME: All (successful) motor events -> For spikes occurring -450.0-ms to 150.0-ms relative to given event...complete Generalized linear mixed-effects model fit by ML Model information: Number of observations 115426 Fixed effects coefficients 16 Random effects coefficients 4908 Covariance parameters 7 Distribution Poisson Link Log FitMethod Laplace Formula: Linear Mixed Formula with 10 predictors. Model fit statistics: AIC BIC LogLikelihood Deviance 2070694.64594067 2070916.74279363 -1035324.32297034 2070648.64594067 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF {'(Intercept)' } 3.57572346737713 0.0733927302748854 48.7204039689571 115410 {'Group_Ischemia' } -0.0195157027956764 0.0715499145860975 -0.272756479285419 115410 {'Alignment_Reach' } 0.00528995465089975 0.00145654620937481 3.63184814656196 115410 {'Alignment_Grasp' } -0.0182196556354626 0.00157593730685186 -11.5611551019493 115410 {'Alignment_Support' } 0.00840461636147413 0.00152407741103524 5.51456002209577 115410 {'Area_RFA' } -0.0681624120627805 0.00976507158863496 -6.98022655994771 115410 {'Group_Ischemia:Alignment_Reach' } -0.000841397803468946 0.00145654612248933 -0.577666433268137 115410 {'Group_Ischemia:Alignment_Grasp' } 0.0178731018064334 0.00157593719920992 11.3412525672939 115410 {'Group_Ischemia:Alignment_Support' } 0.00414584843906602 0.00152407725430637 2.72023509789395 115410 {'Group_Ischemia:Area_RFA' } -0.00117509136133283 0.00909834158655024 -0.129154456353883 115410 {'Alignment_Reach:Area_RFA' } 0.0584730775415553 0.00145573192062348 40.1674763829536 115410 {'Alignment_Grasp:Area_RFA' } 0.00397686347368532 0.00157473972210449 2.52541002037507 115410 {'Alignment_Support:Area_RFA' } -0.0135693896674162 0.00151737632511617 -8.9426659970969 115410 {'Group_Ischemia:Alignment_Reach:Area_RFA' } -0.0377072034252674 0.00145573131529096 -25.9025845148704 115410 {'Group_Ischemia:Alignment_Grasp:Area_RFA' } 0.0255273235491251 0.0015747400454985 16.210500026399 115410 {'Group_Ischemia:Alignment_Support:Area_RFA'} 0.0113326316689496 0.00151737618770691 7.46857091917045 115410 pValue Lower Upper 0 3.43187485069629 3.71957208405797 0.785040903890067 -0.159752429216582 0.120721023625229 0.000281520480224908 0.00243514659885504 0.00814476270294446 6.74257212311319e-31 -0.0213084683927656 -0.0151308428781595 3.50400151747125e-08 0.00541744819821318 0.0113917845247351 2.96285387871772e-12 -0.0873018014076568 -0.0490230227179043 0.56349047250788 -0.00369620568521946 0.00201341007828157 8.49983580927189e-30 0.0147842892601069 0.02096191435276 0.0065245283052763 0.00115868058299123 0.0071330162951408 0.897235657909965 -0.0190077002098263 0.0166575174871607 0 0.0556198654828742 0.0613262896002364 0.0115576644812606 0.000890397963972331 0.0070633289833983 3.85322313086563e-19 -0.0165434238058752 -0.0105953555289573 1.64991718000032e-147 -0.0405604142975061 -0.0348539925530287 4.94133662755614e-59 0.0224408574055648 0.0286137896926853 8.16380531920693e-14 0.00835859779981064 0.0143066655380885 Random effects covariance parameters: Group: ChannelID:Day (3609 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.449333360975598 Group: AnimalID (10 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std' } 0.209582475597 {'Duration' } {'(Intercept)'} {'corr'} 0.0752307877613288 {'Duration' } {'Duration' } {'std' } 0.0427497614205649 Group: SupportLimbMovement (2 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 2.12216162032972e-06 Group: Trial_ID (1272 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.108919270346587 Group: ICMS (5 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.0342556487479956 Group: Error Name Estimate {'sqrt(Dispersion)'} 1 R-squared: Ordinary: 0.298127582313874 Adjusted: 0.298036358968711 ==============================================
analyze.stat.surf_predicted(full_glme,'Area','Group','Tag','Full Model','ZLim',[29.5,46]);

Organize output into table

G = table(Outcome,Align,tStart,tStop,glme);
clear glme;
G.Outcome = string(G.Outcome);
G.Align = string(G.Align);
G.Properties.VariableUnits = {'','','ms','ms',''};
G.Properties.VariableDescriptions = {'Trial result',...
'Behavioral alignment event','Relative start time',...
'Relative end time','Generalized Linear Mixed-Effects Model'};
G.Properties.UserData = struct('type',"SpikeCountsGLME",'modelspec',mdl);
G.Properties.Description = ...
'Raw spike counts Generalized Linear Mixed Effects models';
disp(G);
Outcome Align tStart tStop glme ______________ _________ ______ _____ _________________________________ "Successful" "Reach" -450 150 {1×1 GeneralizedLinearMixedModel} "Successful" "Grasp" -450 150 {1×1 GeneralizedLinearMixedModel} "Unsuccessful" "Grasp" -450 150 {1×1 GeneralizedLinearMixedModel} "Successful" "Support" -450 150 {1×1 GeneralizedLinearMixedModel} "All" "Reach" -1350 -750 {1×1 GeneralizedLinearMixedModel}

Check model residuals

In general, the models are not completely accurate and are likely missing some covariates, but they are still pretty good.
One noted point is that we are analyzing multi-unit spiking (not sorted single-unit activity); therefore our assumption of a Poisson process does not entirely hold since it's more likely that counts on each channel are the result of a hyper- or hypo-exponential process, depending on the architecture of nearby detected units. We chose to use multi-unit activity because the density of the electrode arrays used was not high enough for us to make good use of spatial features in combination with temporal features of the waveforms, as might be done with tetrodes or newer high-density arrays such as Neuropixels.
Another likely source of error that our model does not capture is the relative temporal proximity of the behavioral events of interest; the Reach and Grasp are typically only separated by a few-hundred milliseconds, and this temporal "smearing" of activity may also account for the skewness of the residual distributions.
% "Run Section" button (not shown in html export)
analyze.stat.panelized_residuals(G.glme{1});
R-squared Ordinary: 0.411485098553609 Adjusted: 0.411434740746868
analyze.stat.panelized_residuals(G.glme{2});
R-squared Ordinary: 0.435660016146892 Adjusted: 0.435602645087876
analyze.stat.panelized_residuals(G.glme{3});
R-squared Ordinary: 0.440582123413158 Adjusted: 0.440538569741506
analyze.stat.panelized_residuals(G.glme{4}); % Support-aligned successes
R-squared Ordinary: 0.370722085728685 Adjusted: 0.370625880574386
analyze.stat.panelized_residuals(G.glme{5}); % Pre-reach epoch: All Outcomes
R-squared Ordinary: 0.342094377384128 Adjusted: 0.342069761002434
In general, we notice that the residuals are right-skewed (top row panels) and slightly unbalanced when plotting the Fitted Value vs. Residuals (although my intuition is that the apparently linear lower-bound on the residuals vs. fitted values for lower fitted values probably relates to the model dealing with counts that can be "close" to zero, without any formal explanation for why I think that). We do not see an indication of trends being driven by any particular observation with exceedingly high leverage, which is reassuring; the histogram and normal probability plots are not perfect, but I think given the information we are working with, this is a pretty reasonable model fit.

Fit model with fixed effect of PostOpDay

Rather than using Day as a nested grouping variable by ChannelID, we can use PostOpDay as part of a fixed-effect interaction in order to estimate the significance of these trends over the time-course observed in the experiments.
mdl_time = "Spikes~Group*Area*PostOpDay+(1|ChannelID)+(1+Duration|AnimalID)+(1|SupportLimbMovement)+(1|Trial_ID)+(1|ICMS)";
time_glme = analyze.stat.fit_spike_count_glme(R,'Grasp','Successful',[-450 150],'Model',mdl_time);
============================================== GLME: Grasp trials (Successful) -> For spikes occurring -450.0-ms to 150.0-ms relative to Grasp...complete Generalized linear mixed-effects model fit by ML Model information: Number of observations 29514 Fixed effects coefficients 8 Random effects coefficients 1570 Covariance parameters 7 Distribution Poisson Link Log FitMethod Laplace Formula: Linear Mixed Formula with 9 predictors. Model fit statistics: AIC BIC LogLikelihood Deviance 561149.979437611 561274.368737702 -280559.989718806 561119.979437611 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF pValue {'(Intercept)' } 3.69911026342569 0.0577814215262802 64.0190248996081 29506 0 {'Group_Ischemia' } 0.143189605450723 0.0574256503145869 2.49347816988241 29506 0.0126552699532911 {'PostOpDay' } -0.00939026385266123 0.00194646913883828 -4.82425519382529 29506 1.41228249190328e-06 {'Area_RFA' } 0.000875524195701153 0.0144347834774415 0.0606537809915482 29506 0.951635352083969 {'Group_Ischemia:PostOpDay' } -0.00776263853583118 0.00194647043118978 -3.98805880194454 29506 6.67766500050809e-05 {'Group_Ischemia:Area_RFA' } 0.175225994567408 0.014302548648283 12.2513825246414 29506 2.00268660749331e-34 {'PostOpDay:Area_RFA' } -0.00265123745214568 0.00025642409082784 -10.3392682161275 29506 5.16497609446972e-25 {'Group_Ischemia:PostOpDay:Area_RFA'} -0.00891735582962774 0.000256424168160704 -34.7758009457172 29506 9.99296598417254e-260 Lower Upper 3.5858561124669 3.81236441438447 0.0306327818586155 0.255746429042831 -0.0132054297634955 -0.00557509794182694 -0.0274172921426441 0.0291683405340465 -0.0115778069797318 -0.00394747009193059 0.147192364363774 0.203259624771042 -0.00315384005216785 -0.0021486348521235 -0.00941995858122576 -0.00841475307802972 Random effects covariance parameters: Group: ChannelID (272 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.20694540072534 Group: AnimalID (10 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std' } 0.142525974513132 {'Duration' } {'(Intercept)'} {'corr'} -0.363684925544796 {'Duration' } {'Duration' } {'std' } 0.150184815824692 Group: SupportLimbMovement (2 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 9.82667238649771e-06 Group: Trial_ID (1271 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.251184558612127 Group: ICMS (5 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.0112153511306786 Group: Error Name Estimate {'sqrt(Dispersion)'} 1 R-squared: Ordinary: 0.246824655009931 Adjusted: 0.246645971778895 ==============================================
Model with fixed-effect of PostOpDay. Each interaction term (and the main effect) containing PostOpDay tests as statistically significant. The three-way interaction Group*PostOpDay*Area indicates that as the time from injury increases, the likelihood of observing spikes in RFA of Ischemia rats decreases during alignment to the Grasp, relative to RFA activity in Intact rats or in CFA of Ischemia rats (p < 0.001; tStat = -34.8; DF = 29506). The Fixed terms from the previous model that used Day as a nested Grouping variable (the Random-Effects model) now have a different significance: there is a statistically significant effect of the main fixed-effect of Group (p = 0.013; tStat = 2.49; DF = 29506); Area, which was a statistically significant effect in the previous model, is no longer statistically significant.

What's going on?

We expect the R-squared ordinary and adjusted values to be lower by comparison to the Random-Effect model, since we are no longer accounting for some of the spread in variance as a grouped random variable. Indeed, R-squared is considerably lower (by approximately 0.2). The significant difference in model R-squared indicates that there is considerable variance in the "patterns" of fit by ChannelID, when it comes to shift in Channel-specific trends of activity by day. We can use this fixed-effect model to study the effect of post-operative day on trends by Group and Area, but when it comes to differentiating between group means that are specific to Group*Area combinations, we are better off using the previous Random-Effect model. We could also test for the effect of PostOpDay using a random effects design matrix and the previous Random-Effect model, but in Matlab using the Fixed-Effect model will be more straightforward.
analyze.stat.panelized_residuals(time_glme);
R-squared Ordinary: 0.246824655009931 Adjusted: 0.246645971778895
The distribution of residuals is right-skewed; we could use a square-root transformation to help account for variability in spike counts that may arise due to the use of multi-unit activity, but this ultimately shouldn't make much of a difference in the end result, especially since we are using a log link already in the glme.
analyze.stat.surf_predicted(time_glme,'Area','Group','Tag','Fixed-Effect of PostOpDay (Area*Group)','ZLim',[29.5,46]);
We expect the optimization procedure to work similarly for the Fixed-Effect model by comparison to the Random-Effect model; indeed, the predicted model outputs are identical.

Partial-dependence plots for time-course of effects

analyze.stat.qPDP(time_glme,{'PostOpDay','Group'});
A partial-dependence plot for PostOpDay*Group shows us that the increase in mean spike counts around the Grasp for Ischemia rats is highest in the earliest days after the injury.
analyze.stat.qPDP(time_glme,{'PostOpDay','Area'});
A partial-dependence plot for PostOpDay*Area shows us that the increase in mean spike counts around the Grasp for Ischemia rats at early time-points is* related to an increase in the expected number of spikes in RFA at early points. So the overall model prediction, which indicates that the ipsilesional cortex premotor area and contralesional cortex primary motor area have a roughly equal number of expected spikes around the Grasp (and a much-different ratio by comparison to the Intact model) is due to an initial increase in RFA spiking around the Grasp (during Successful pellet retrievals only).