Electrophysiology data analysis

Raw signals recorded with the 32-channel silicon probe were downsampled to 2 kHz for LFP analysis. Signals were imported into MATLAB software, version R2017a (MATLAB; RRID:SCR_001622) and processed using a combination of custom and modified routines from the freely available MATLAB packages FMAToolbox (FMAToolbox; RRID:SCR_015533) and Chronux (RRID:SCR_005547; Hazan et al., 2006; Mitra and Bokil, 2008).

For comparative analysis of neurophysiological signals between the injured and sham groups of animals, we first bandpass filtered (0.6–6 kHz) the LFP for each channel of the 32-channel HC probe, giving an estimate of the multiunit spiking activity. We then selected the channel with the maximum average power for each animal and designated that channel to be the pyramidal cell layer. HC probes from all animals were then aligned to this channel. This procedure resulted in 25 aligned channels of neurophysiological data common to all animals.

Current source density (CSD) of the LFP across all HC channels was calculated with the spline inverse CSD method of Pettersen et al. (2006), using the freely available MATLAB package CSDplotter version 0.1.1 (https://github.com/espenhgn/CSDplotter) visualized by convention as a heat map of current sinks (blue) and sources (orange; Mitzdorf, 1985).

To investigate any statistical differences in power between the two groups and among each of the 25 channels, we windowed each 15 min recording in 1 s epochs and calculated the power spectral density in the 1–300 Hz frequency band (Welch's method: frequency resolution, 0.5 Hz; overlap, 60%) for each epoch. As a result, we obtained a 25 (channels) × 600 (frequency bins) × 2250 (epochs) matrix of power values for each subject. We then calculated the channel–frequency map power difference distribution (PDD) relative to the sham group [(injured/sham − 1)%], averaging the aforementioned matrix along the epochs and groups. To evaluate where these differences were statistically significant, we calculated a significance mask through a repeated-measures ANOVA between the two groups using the epochs as independent variables, obtaining a channel–frequency map of p values: p values <0.05 were considered significant and converted to a value of 1, otherwise values were represented with a 0. Finally, we multiplied the significance mask by the PDD to identify where the two groups displayed statistically significant differences.

To investigate the interaction between oscillations in different frequency bands after injury, we used phase–amplitude coupling (PAC) as a form of cross-frequency coupling analysis in the HC. We investigated how the phase of lower-frequency oscillations drive the power of coupled higher-frequency oscillations by measuring the degree of synchronization of the amplitude envelope of faster rhythms with the phase of slower rhythms. To quantify the PAC between strata radiatum and pyramidale LFP, we used the Kullback–Leibler (KL)-based MI based on the joint entropy (Tort et al., 2010).

Briefly, MI was calculated as follows: components of the stratum pyramidale (SP) LFP for each frequency [1–300 Hz, denoted as Pyr(t,fpyr),fpyr=1,2,3...299,300 Hz] were extracted by wavelet transform (Morlet; 10 cycles). The phases of the bandpass filtered (1–20 Hz) stratum radiatum LFP (denoted as ΦRad(t,frad),frad=1,2,3...20 Hz) were obtained from the standard Hilbert transform. The Hilbert transform was also applied to extract time series of the amplitude envelopes (denoted as APyr(t,fpyr)) of Pyr(t,fpyr). The composite time–frequency series (ΦRad(t,frad), APyr(t,fpyr)) was then constructed, which gave the amplitude of the pyramidal layer LFP for each frequency bin oscillation (1–300 Hz) at each phase of the radiatum layer LFP for each frequency bin (1–20 Hz). Next, the phases ΦRad(t,frad) were binned and the mean value of APyr(t,fpyr) over each binned phase was calculated. We denoted [APyr(t,fpyr)]ΦRad(t,frad)(j) for the mean APyr(t,fpyr) value at the phase bin j. Afterward, the normalized mean amplitude of [APyr(t,fpyr)]ΦRad(t,frad)(j) was divided for each bin value by the sum over all bins. The phase–amplitude distribution for each combination (fpyr,frad) (P(j,fpyr,frad)) was obtained by the following equation, where N was the number of phase bins:

The joint entropy (H) of Rad(t,frad)−Pyr(t,fpyr) coupling was calculated as follows:

When the joint distribution of ΦRad(t,frad) − APyr(t,fpyr) was uniform, the maximum value of joint entropy H0 was obtained, as follows:

We then calculated the KL distance by subtracting the value H ₀ from H, and measured how much the ΦRad(t,frad) – APyr(t,fpyr) joint distribution deviated from the uniform joint distribution. Finally, MI was obtained by dividing the KL distance by the constant factor H ₀, which made the measure assume values between 0 and 1, as follows:

Thus, a larger MI value denoted a stronger coupling between the stratum radiatum and stratum pyramidale layers of the HC.

To investigate any statistical differences in coupling between the two groups, we calculated the t test of MI(fpyr,frad) for each combination of fpyr and frad.

Automated offline spike detection and sorting were performed on the wideband signals using the Klusta package (KlustaKwik; RRID:SCR_014480), developed for higher-density electrodes, and manually refined using the KlustaViewa (https://github.com/klusta-team/klustaviewa), or phy (https://github.com/kwikteam/phy) software package. The Klusta routines are designed to take advantage of the spatial arrangement of and spike timing from all probe channels simultaneously (32 channels for HC probe, 4 channels for MS tetrode) in constructing putative clusters (Rossant et al., 2016). Putative single units were required to have at least 100 spikes and a refractory period violation rate of <2% (interspike interval of <3 ms). For putative units isolated on the 32-channel HC probe, we further required a quality index of >0.95, a measure of cluster uniqueness relative to all other potential clusters in the dataset computed by the KlustaViewa software (Rossant et al., 2016). For putative units isolated on the MS tetrode, we additionally required the isolation distance (Harris et al., 2001) of each putative single unit to be ≥20 (Harris et al., 2003). The resulting curated putative single-unit clusters were then imported into MATLAB software, version R2017a, for visualization and further analysis using custom and built-in routines (MATLAB; RRID:SCR_001622). Wideband signals were bandpass filtered (0.6–6 kHz) for average waveform visualization only.

We created a three-step process for classifying isolated HC putative single units on the 32-channel probe as either pyramidal cells or interneurons. In the first step, putative units were assigned to an anatomic layer (stratum oriens, stratum pyramidale, or stratum radiatum). We defined a five-channel span (100 μm) centered on the channel of maximum multiunit activity for each rat and defined this as stratum pyramidale. For each putative unit, if the channel of maximum amplitude fell within this band, it was assigned to stratum pyramidale. If the maximum amplitude was above this band, it was assigned to stratum oriens, and if below, it was assigned to stratum radiatum. If the maximum amplitude of a putative unit fell on a channel adjacent to the stratum pyramidale band, it was classified temporarily as a border unit. All nonborder putative units outside of stratum pyramidale were then classified, by definition, as interneurons, regardless of firing rate. In the second step, we modified the method of Csicsvari et al. (1999) of putative unit classification into pyramidal cells, or interneurons and applied it to all units assigned to stratum pyramidale or classified as border units. This approach computes the following three features of each putative unit: overall firing rate, a measure of averaged spike waveform width, and a measure of “burstiness” in the spike train of the unit. Spike width was calculated by extracting the waveform for each spike from the unfiltered wideband signal, which was then averaged and upsampled by a factor of 5. The peak amplitude of the averaged waveform was then detected, and its height was calculated by subtracting a baseline amplitude 1 ms before the peak. The spike width was then calculated as the time interval between the points along the averaged waveform that fell at 25% of the spike height, centered on the peak. Burstiness was defined as the first moment of the autocorrelogram, with smaller values indicating higher burstiness. The distributions of these three features across all putative stratum pyramidale units were plotted and inspected for natural clusters. Each feature revealed two clusters, and we thus made a simple division into two clusters for each feature using threshold values (see Fig. 4B). Putative pyramidal cells were defined as those units with a spike width above the determined threshold, a first moment of the autocorrelogram below the determined threshold (equating to high burstiness), and a firing rate below the determined threshold (division was at 3.5 Hz). Similarly, putative interneurons were defined as those units with a spike width below threshold, a first moment above threshold (equating to low burstiness), and a firing rate >3.5 Hz. In a final third step, those putative units classified as both interneurons and border cells were reassigned to stratum pyramidale if there were also isolated putative pyramidal cells with maximum amplitudes on the same channel, otherwise they were assigned either to stratum oriens or stratum radiatum, depending on which border they were located.

Neuronal firing rates and putative neuron number are preserved in CA1 following TBI. A, B, Division of isolated CA1 single units within SP into putative pyramidal cells and putative interneurons. A, Dark gray bars indicate histograms of the first moment of the autocorrelogram (x-axis), with lower numbers indicating higher burstiness, and spike width (at 25% maximum height) of the unfiltered waveforms (y-axis). B, Firing rate histogram of all isolated single units. Each feature reveals two distributions of single units, defining the following two clusters: magenta circles in the right bottom quadrant (narrow, tonic, fast spiking cells—putative interneurons); and green circles in the left top quadrant (wide, bursting, slower-spiking cells—putative pyramidal cells). Gray circles represent units that could not be subtyped. C, Example of single isolated putative pyramidal cell (left) and putative interneuron (right). Top, Mean (±SD) waveforms. Bottom, Autocorrelograms expressed as firing probability as a function of time bin (bin size, 2 ms). Note in the putative pyramidal cell the wider waveform as well as the multiple short-latency peaks in the autocorrelogram that correspond to burst activity compared with the putative interneuron. D, Overall mean (±SEM) firing rates of putative pyramidal cells (left) and interneurons (right) demonstrating no change in baseline activity between sham and injured animals. E, No change in the median (with 95% CI) number of isolated putative pyramidal cells (top) or interneurons (bottom) per rat after injury. F, Pie charts depict the relative distribution of putative neuronal subtypes across CA1 layers. G, Laminar distribution of all isolated putative pyramidal cells (circles) and interneurons (stars) as a function of recording channel normalized to SP. Dotted lines indicate SO–SP and SP–SR borders. H, Plots compare median (with 95% CI) number of isolated putative interneurons per rat between groups in each layer. There is no significant dropout of putative interneurons within any layer. Blue and orange represent sham and injured rats, respectively. Circles and stars represent putative pyramidal cells and interneurons, respectively. inj, Injured; pyr, pyramidals; int, interneurons.

To classify isolated single units in the MS as either putative GABAergic/glutamatergic (GABA/Glu) or putative cholinergic, we modified a method first described using intracellular and juxtacellular recordings in slice preparations (Matthews and Lee, 1991) and subsequently extended to in vivo recordings (King et al., 1998). This method relies solely on waveform shape by capturing an observed “inflection” of the repolarization phase of the action potential found to be characteristic of cholinergic neurons, but not GABAergic/glutamatergic neurons. Spectral power analysis of the waveforms from 1 to 16 kHz reveals increased power in higher-frequency bins in the “inflected” waveform group. We found that the distribution of the ratio of the fourth to second frequency bin revealed two natural clusters in both the sham and injured groups. The lower ratio, or “noninflected,” waveform cluster was then classified as putative GABAergic/glutamatergic cells. The higher ratio, or inflected, waveform cluster was then further refined by firing rate, as we expect cholinergic units to have firing rates generally <4 Hz. The firing rate distribution of the inflected subgroup revealed a low-firing cluster that was then classified as the putative cholinergic subgroup.

For entrainment of putative HC or MS single units to the low-frequency oscillations, we bandpass filtered [8–20 or 2–5 Hz (“anesthesia theta”)] the LFP of the channel 425 μm below the center of stratum pyramidale, which corresponded to the channel of maximal power difference between groups within the 8–20 Hz range, and extracted the phase angles for each spike within each putative single unit using the Hilbert transform. Directional statistics were performed using MATLAB Circular Statistics Toolbox (CircStat; RRID:SCR_016651). For each putative unit, the mean resultant vector was calculated, giving the mean phase angle of entrainment for that unit, as well as the resultant vector length, a measure of spread of the individual phase angles about the mean. Putative units were then considered significantly entrained if the resultant vector length was greater than the threshold of nonuniform distribution calculated by the Rayleigh test (Berens, 2009). The mean angles of entrainment between significantly entrained sham and injured putative units as well as their distributions were compared using the two-sample Watson–Williams test and the two-sample Kuiper test, respectively.

For entrainment of putative HC single units to fast oscillations in stratum pyramidale we first extracted components of the stratum pyramidale LFP by wavelet transform (Morlet, 10 cycles) for each frequency (1–400 Hz, denoted as Pyr(t,f),f=1,2,3,...399,400 Hz). The phases for each frequency bin (denoted as ΦPyr(t,f),f=1,2,3,...399,400) were obtained from the standard Hilbert transform. Next, the phases ΦPyr(t,f) were binned (18 bins, 20°), and the number of units per bin [#units(j)] was calculated. We denoted the normalized unit distribution over the phase bin j as follows:

The joint entropy H was calculated as follows:

We then calculated the MI for putative single unit–stratum pyramidale LFP entrainment, as per methods detailed above (see Phase–amplitude coupling section). A larger MI value denoted a stronger entrainment. As above, we used the Student’s t test to determine frequencies of significant differences in MI.

To investigate the influence of HC oscillatory activity on MS neuron firing patterns, we computed MS spike-triggered averages of the HC LFP recorded simultaneously with spiking activity in the MS. For each putative MS neuron, we averaged 1 s windows of the unfiltered HC LFP centered on each putative MS cell spike. Spike-triggered averages of putative units meeting criteria for “CA1 modulated” (see below) were then averaged within the sham and injured groups. To quantify differences in the averaged spike-triggered averages between the sham and injured groups, we calculated the absolute difference between the maximum positive and maximum negative voltage deflections in the averaged tracings for each unit, making between-group comparisons using Welch’s t test.

For each putative CA1-modulated MS unit the concurrently recorded HC LFP within stratum radiatum was bandpass filtered from 8 to 20 Hz or from 2 to 5 Hz (anesthesia theta). The troughs of the resulting oscillation were detected, and the times of these events were used as a trigger for the creation of the peristimulus time histogram (PSTH; 100 ms window; bin size, 2 ms) for each putative CA1-modulated MS unit. Because individual MS units had different overall firing rates for each bin, the number of spikes was divided by the total number of spikes for that putative unit to derive a firing probability that could be compared across units.

Custom MATLAB scripts written for electrophysiology data analysis are described above. Custom scripts and data presented in this article are available on request.