Temporal activity patterns of layer II and IV rat barrel cortex neurons in healthy and injured conditions

Abstract Neurons are known to encode information not just by how frequently they fire, but also at what times they fire. However, characterizations of temporal encoding in sensory cortices under conditions of health and injury are limited. Here we characterized and compared the stimulus‐evoked activity of 1210 online‐sorted units in layers II and IV of rat barrel cortex under healthy and diffuse traumatic brain injury (TBI) (caused by a weight‐drop model) conditions across three timepoints post‐injury: four days, two weeks, and eight weeks. Temporal activity patterns in the first 50 ms post‐stimulus recording showed four categories of responses: no response or 1, 2, or 3 temporally‐distinct response components, that is, periods of high unit activity separated by silence. The relative proportions of unit response categories were similar between layers II and IV in healthy conditions but not in early post‐TBI conditions. For units with multiple response components, inter‐component timings were reliable in healthy and late post‐TBI conditions but disrupted by injury. Response component times typically shifted earlier with increasing stimulus intensity and this was more pronounced in layer IV than layer II. Surprisingly, injury caused a reversal of this trend and in the late post‐TBI condition no stimulus intensity‐dependence differences were observed between layers II and IV. We speculate this indicates a potential compensatory mechanism in response to injury. These results demonstrate how temporal encoding features maladapt or functionally recover differently in sensory cortex after TBI. Such maladaptation or functional recovery is layer‐dependent, perhaps due to differences in thalamic input or local inhibitory neuronal makeup.

of depression (Kreutzer et al., 2001) and epilepsy (Liesemer et al., 2011;Yeh et al., 2013), likely caused by damage to and/or reorganization of inhibitory microcircuits (Cohen et al., 2007;Gaetz, 2004;Hunt et al., 2010;Povlishock & Katz, 2005;Reeves et al., 1997;Werner & Engelhard, 2007). The specific loss and damage of axons and cells after diffuse TBI will differ between individuals and progress at different rates. Subsequent "functional recovery" (i.e., gradual improvement in behavioral assessment outcomes) is typically slow and generally does not return to pre-injury functional levels, despite multifaceted and interacting regenerative processes being triggered post-injury Nudo, 2013). One hypothesis (Alwis et al., 2012Greer et al., 2013) is that maladaptations of the recovering cortex, like potential synaptic changes between particular cell types due to specific cell type(s) death and circuit remapping, cause disruption to microcircuit function in such ways that cause excitation-inhibition changes. This would partially explain damage and (functional) recovery differences observed between cortical layers and cell types (Cantu et al., 2015;Carron, Yan, et al., 2016). However, how such broad changes affect the encoding of sensory stimuli in sensory cortices during and after injury is not completely understood.
Multi-unit activity from extracellular physiological recordings from sensory cortex can be used to directly measure and compare neural activities in terms of mean amplitude over stimulus presentation, response onset delay, or the delay to maximum response (Allitt, Iva, et al., 2016;Maravall et al., 2007). Such methods have been used to study pooled neural population responses in TBI barrel cortex models (Allitt, Iva, et al., 2016;Allitt, Johnstone, et al., 2016;Alwis et al., 2012Alwis et al., , 2016Alwis & Rajan, 2014;Carron, Yan, et al., 2016;Johnstone et al., 2013Johnstone et al., , 2014Johnstone et al., , 2015Yan et al., 2013). Barrel cortex is an especially rich sensory area to study due to its ecological importance for rodents and the ability interpret and connect results to the significant number of past studies (Burns & Rajan, 2021). It is conveniently arranged into primary cortical columns which process the sensory information from individual whiskers from the rodent face, and the arrangement of these columns in barrel cortex maps to the arrangement of whiskers on the face. In this study, we re-analyze data from such studies on a unit-byunit basis, identifying and quantifying individual units' temporal activity patterns in response to one basic and one naturalistic stimulus in healthy and perturbed states.
For four reasons, we here focus only on layers II and IV:(1) layers II and IV have a very different proportions and types of interneurons (Markram et al., 2004), which are known to be affected differently by diffuse TBI ; (2) layer IV is the main input layer and layer II is one of the first and most sophisticated signal processing layers in barrel cortex (Burns & Rajan, 2021); (3) data from TBI-affected rats consistently show more superficial layers have medium-and long-term changes in their excitability compared to deeper layers (Allitt, Iva, et al., 2016;Johnstone et al., 2013); and (4) recordings of layers II and IV have been the most frequent in past studies of TBI barrel cortex. Future studies may extend our analyses to other layers. Although we do not perform waveform analysis or otherwise attempt to identify units as "fast spiking" or "regular spiking," in relation to reason (1) above, we expected units from layer II to show greater perturbations to their temporal responses properties than layer IV given interneurons' strong influence in various temporal coding schema in barrel cortex (Burns & Rajan, 2021;Markram et al., 2004).
Our analyses allow us to observe how temporal activity patterns transform between cortical layers IV and II and, when interpreted with knowledge of the neuroanatomy, indicate how sensory information is processed in these layers of barrel cortex in healthy and perturbed states. Such comparisons may help to establish the importance of specific interneuron subtypes, microcircuits, and the maladaptations relevant to TBI pathophysiology. We also present a novel graph theoretic method for studying temporal activity patterns of neurons.

| Electrophysiological database
We collated electrophysiological data from previous studies (Allitt, Iva, et al., 2016;Allitt, Johnstone, et al., 2016;Alwis et al., 2012Alwis et al., , 2016Carron, Yan, et al., 2016;Johnstone et al., 2013Johnstone et al., , 2014Johnstone et al., , 2015Yan et al., 2013). The database is composed of data from 36 male Sprague-Dawley rats, across three timepoints after traumatic brain injury (TBI) or sham surgery: four days (sham = 4, TBI = 6), two weeks (sham = 4, TBI = 9), and eight weeks (sham = 5, TBI = 8). TBI was caused by a modified version of the weight-drop impact acceleration method (Foda et al., 1994;Hellewell et al., 2010), which aims to make the primary injury a diffuse TBI. A hallmark of this type of injury is that it does not typically result in any sizeable lesion in cortex, making it difficult injury to simply diagnose for clinicians and grossly measure for researchers. An impact velocity of 6.15 m/s was used to model a severe injury . At different timepoints after the injury, animals underwent terminal surgery where they were anesthetized with 5% halothane and tracheomatized to maintain anaesthesia at 0.5%-3% and allow continuous ventilation. Body temperature was maintained at 37-38°C. Depth of anaesthesia was monitored using ECG/EMG recordings from forepaw musculature, hind paw pinch, and palpebral reflexes. A craniotomy 5 mm in diameter was performed over the right barrel cortex (2 mm caudal to bregma; 6 mm lateral to the midline) and a microelectrode (2-4 MΩ; FHC, Bowdoinham, ME) was used to record from a single cortical column with the same principal whisker (PW). This PW was identified by manually stimulating individual whiskers with a toothpick until consistent, strong drive was found for a single whisker. The PW was then threaded with a computer-controlled motor lever arm system to present precise stimuli protocols by deflecting the whiskers in the dorsal-ventral direction, perpendicular to the face. Output signals from the microelectrode were amplified and bandpass-filtered from 0.3 to 10 kHz. The microdrive controlling the depth of the inserted microelectrode was zeroed at the pia and advanced to different cortical depths from the pia. So as to limit the effects of the anaesthesia over time, some experimental recordings were from deep to superficial cortical depths whereas others were from superficial to deep (and in both cases moving in increments of 50 µm or 100 µm).
Some animals in the database underwent behavioural assessments for sensorimotor function and/or were studied using histological techniques post-mortem. We were not able to systematically analyze these behavioural and histological data in this larger, combined dataset due to differences or absences in the use of these techniques among the animals under study. This study therefore focussed on analysis of the electrophysiological data, which was consistent in method across all animals. All experiments were conducted in accordance with the Code of Practice for the Care and Use of Animals for Scientific Purposes (National Health and Medical Research Council, Australia) and were approved by the Monash University Standing Committee on Ethics in Animal Experimentation. Further details of the injury model and experimental techniques can be found in in the original studies (Allitt, Iva, et al., 2016;Allitt, Johnstone, et al., 2016;Alwis et al., 2012Alwis et al., , 2016Carron, Yan, et al., 2016;Johnstone et al., 2013Johnstone et al., , 2014Johnstone et al., , 2015Yan et al., 2013).
From these studies, we analyzed online-sorted units responding one of two stimuli (depicted in Figure 1a): a. "Basic": A trapezoid stimulus had an on-ramp of varying velocity (30, 60, 150, 250, or 400 mm/sec) to a distance of 3.6 mm from the rest whisker position, a hold period of 20 ms, and an off-ramp over 40 ms back to rest (Allitt, Iva, et al., 2016;Alwis et al., 2012); and b. "Contact": A naturalistic object contact stimulus which lasted a total of 100 ms and included an initial deflection followed by complex movement, reconstructed from high-speed video-tracking of awake behaving rats' whiskers contacting an object and then brushing past it (Hartmann et al., 2003).
All analyses from this database used online-sorted units, sorted by highly-experienced electrophysiologists using the Cambridge Electronic Design Spike2 software. The relative amplitudes of spikes were visualized and shape-based criteria (including rise time of the upstroke, the width of the action potential waveform, the size of the overshoot, etc.) were used to accept and reject units, ensuring units analyzed were highly reliable and of good quality. An example unit is shown in Figure 1b. Waveform analysis to functionally classify units (e.g., into categories of "fast spiking" and "regular spiking") was not performed, although this would be an interesting although this would be an interesting future study to consider in combination with other techniques (e.g., histological) to confirm cells' functional types and other anatomical information, and how this relates to their activity in TBI. The current database amounted to a total of 1210 units from layer II (150-300 µm from pia) or layer IV (750-1000 µm deep from pia) across both health and TBI conditions (summarized in Table 1). It should be noted that our classification of units into layer II and IV using depth estimation was not confirmed using histology, which represents a limitation of this data an analysis. Unit activity profiles were then quantified using established methods (Allitt, Iva, et al., 2016;Alwis et al., 2012).
Lastly, it is pertinent to note a few general limitations and confounding factors regarding this database. In general, diffuse TBI of the kind caused in these experiments causes a period of general period of hypoexcitation (especially in the supragranular layers) of barrel cortex. Many details are still under investigation, although we hope the current work helps to further address this question. Nevertheless, it is possible some of our results or interpretation thereof are confounded by unmeasured injury effects. Additional unmeasured effects may come from the anaesthesia. Furthermore, although rodent models are common animal model for biomedical and neuroscientific research, many of the structures in the rodent brain do not exist in humans, like barrel cortex, and the relationship of layers II and IV to primate granular layers is unclear. This means more translational work is necessary for the current work to have significant clinical impacts for human brain injury patients.

| Unit activity profiles
Unit activity profiles were produced by calculating the peristimulus time histograms (PSTHs) for each unit at each velocity (for Basic) or amplitude (for Contact) across all the repetitions of that stimulus at that layer for that recording. Following the convention of and to allow for easier comparison to previous studies, spikes were counted in 1 ms bins up until 50 ms post-stimulus onset. We then manually examined unit activity profiles in the form seen in Figure 1d-i. A diversity of response patterns existed. Most notably, some units appear to be non-responsive whereas other units contain one or more response components, that is, one or more times at which the unit firing rate is markedly higher than the apparent baseline. Furthermore, some of these components appeared to be relatively stationary in time whereas others are non-stationary and appear to occur earlier or later depending on the stimulus intensity-velocity (for Basic) or amplitude (for Contact).

| Independent neural response components analysis
To objectively identify and quantify response components, an independent neural response components analysis was developed by extending the traditional local maxima algorithm. This analysis was done on a unit-by-unit basis and involved first identifying the local maxima in each PSTH for a F I G U R E 1 (a) Stimulus waveforms for "Basic" and "Contact," the two tactile stimuli presented to the rodents which were studied in this project. (b) Representative waveform of an online-sorted unit. The dark line indicates the mean normalized amplitude over 5 ms for all spikes sorted into this unit, with the shaded region showing the range of all spikes. (c) Example unit's PSTH in response to Contact at a maximum amplitude of 2 mm and (d) the PSTH converted into mean firing rates (Hz), represented in a heatmap across all amplitudes for the same unit. This unit response consists of two major components, one occurring at ~10 ms post-stimulus onset and the second occurring at ~23 ms post-stimulus onset across most amplitudes. (e-i) Example units responding to the Contact stimulus identified as having 0 (e), 1 (f-g), or ≥2 components (h-i) by the independent neural response components analysis (extended local maxima algorithm). Units with 0 components are classed as non-responders and although they may show some tendencies to spike in general periods, for example, offset periods, these responses are not as reliable as those in units with ≥1 components. Units with ≥2 exhibit a wide diversity of component strengths and timings given unit and noting the times at which these local maxima occurred. Local maxima were discarded if the mean firing rate for that 1 ms bin was <25% of max firing rate and the neighboring 1 ms bins had no spikes (or 0% of max firing rate). Local maxima occurring within 5 ms of each other were combined and the median 1 ms time bin between these bins was taken as a single local maximum. These combinations were performed iteratively, starting with the two local maxima nearest each other. This left an approximation of the times where mean firing rates in the PSTH were highest and at least >25% of max firing rate. This analysis was conducted independently across all stimulus velocities (for Basic) or amplitudes (for Contact). We then used the same combination method to take the median times at which local maxima occurred and combine them across amplitudes to approximate, for the whole unit, when the major components of activity occurred. This resulted in units being identified as having either none or one or more neural response components. Kolmogorov-Smirnov (KS) tests were performed to determine the probability of component timing distributions being significantly different between the layers and conditions. For units with multiple components, linear regressions were conducted to determine if the times at which independent components occurred within the same unit were related. We then plotted circle charts showing the relative proportions, for each layer and animal condition (and again independently for both Basic and Contact), of units which were non-responders (having no response components) and for units which had one or more components. These relative proportions were then compared between layers and conditions using KS tests.

| "Peakiness" and significance of components analysis
A concern about the definition of components using the independent neural response components analysis described above was that some components, while being "local maxima," may or may not be statistically significant. Therefore, we used two methods to analyze how statistically valid the components we identified were: (i) a "peakiness" detection algorithm (Beveridge et al., 1989;Haralick & Shapiro, 1984;Jain et al., 1995); and (ii) a t test (which tests if the means of two populations are equal) comparing the components' firing rates with the surrounding, non-component firing rates. Both methods were used on a component-by-component basis for all units which had been identified as having at least one component. For these tests, we took the mean firing rates at the 5, 1 ms time bins occurring immediately before and after the identified component time. Together with the identified component time bin, this represented an 11 ms segment of mean firing rates for that unit, where the 6th ms bin represented the identified component time.
For the peakiness detection algorithm, we identified the local minima within the 11-element segment and called this g k . Then, peakiness = g i + 1 g k + 1 , where g i is the mean firing rate of the 6th ms bin. This provides us with a single number, "peakiness," which described every component. This number varied on a scale of 0 to 2, where 0 represented that the identified component time had infinitely less firing than the local minima within ±5 ms from the component time, where 2 represented it having infinitely more, and where 1 represented it being equal. In practice, this meant that a peakiness value of 1.25 indicates that the mean firing rate at the identified component time was >25% than the surrounding mean firing rate minima within ±5 ms from the component time; a value of 1.5 would indicate >50%, and so on.
For the t tests, we took the 5th to 7th ms bin mean firing rates and called this segment the 'component segment'. We then performed a one-tailed two-sample t test, comparing the mean firing rates of the component segment with the remaining 8, 1 ms bin mean firing rates of

| Graph-theoretic components analysis
The independent neural response components analysis is effective at identifying single, stationary components. However, some components noticeably shifted in time as a function of velocity (for Basic) or amplitude (for Contact). The previous analysis therefore was unable to capture information about how stationary or non-stationary these components were (e.g., Figure 2a). The non-stationarity of components could be an important feature to measure since, for units which had more than one component, some components appeared to interact, for example, while ascending in amplitude one component might increase its activity while the other simultaneously and proportionally decreases (e.g., left panel of Figure 2a). Such interactions, if sufficiently separated in time, could be evidence of feedforward inhibition or disinhibition, known to be present, for example, in layer IV of barrel cortex (Burns & Rajan, 2021). To identify these features, we developed a novel graph-theoretic analysis. Such analysis has the ability to more accurately and precisely capture timing information rather than aiming for a general, normalized analytic approaches such as latency, PSTH (peak-widths), response duration, and firing rates. PSTHs for each unit are treated as a matrix and transformed into weighted directed graphs. Vertices represent the neural activity at different velocities/amplitudes and times, and are connected to neighboring vertices via weighted, directed edges. Vertices receive weighted, directed edges from other vertices representing the neural activity at velocities/amplitudes of z ± 1 and at times of t ± 1 (including at combinations z ± 1, t ± 1), where z is F I G U R E 2 (a) Example unit activity profiles with red lines showing the component times approximated by the independent neural response components analysis. (b) Mean firing rates mapped in the amplitude-time domain for two healthy layer IV units responding to the Contact stimulus. Coloured stars are plotted over the amplitude-time bins along the three shortest, non-overlapping paths calculated via the graphic-theoretic method described in the text (black = primary shortest path, red = secondary shortest path, white = tertiary shortest path). (c) Visualisation of the graph-theoretic algorithm procedure to identify neural response components. Each black or orange vertex is equivalent to one ms bin at one amplitude in the PSTH matrix of a single unit. For example, for the Contact stimulus there are 50, 1 ms bins (only 10 are shown in this illustration) arranged in 10 rows of increasing amplitude. Each vertex is connected to its neighbours (in the time and intensity domains) via edges. Edges arriving at a vertex are weighted according to the firing rate at that vertex. A vertex with a high firing rate has a small edge weight and a vertex with a low firing rate will have a large edge weight. Using Dijkstra's shortest path algorithm, the shortest possible path (travelling along edges which accumulate the lowest weighted score, as determined by the sum of their weights) is found (orange path) along the path with high firing rates. The top green vertex and bottom red vertex are the start and goal vertices, a technical requirement of the Dijkstra's algorithm, and allow the path to start and finish at any vertex in the top and bottom rows, respectively the receiving vertex's velocity/amplitude and t is its time. These edges are made for all vertices that exist (Figure 2c). The edges are weighted according to the mean firing rate of the vertex they are directed towards so that edges leading to vertices with relatively high mean firing rates have lower weights than edges leading to vertices with relatively low mean firing rates, like so: unit max firing rate vertex mean firing rate . A start and goal vertex are connected to the top and bottom rows of the PSTH matrix graph, respectively, and used to calculate the shortest path between these vertexes through the PSTH matrix graph using Dijkstra's shortest path algorithm (Dijkstra, 1959). This path preferentially travels through the vertices representing the highest mean firing rates (since edges leading to them had the least weight) but will ensure such vertices are as contiguous as possible, owing to the graph's structure. Once a shortest path is calculated, we re-weight the edges leading to vertices in that shortest path with an arbitrarily large weight and re-run the shortest path algorithm to find the secondary shortest path, which does not re-visit any vertices from the primary shortest path. This process can be repeated as many times as desired and here we calculate only the first three paths. A copy of algorithm implemented in MATLAB is freely available at https://github.com/tfbur ns/graph -theor eticident ifica tion-of-PSTH-compo nents. Figure 2b shows an example of the three shortest, non-overlapping paths through two units' firing activity profiles, as mapped in the amplitude-time space for the Contact stimulus. Once identified, the following features of these paths were calculated: (i) the median time bin which the path vertices represent; (ii) the mean firing rates at each path step (and how these change as a function of mean firing rates at the equivalent steps of other components' within the same unit); and (iii) how stationary or non-stationary the path is in time.

| Component timings using the independent response components analysis
For both stimuli and layers there is a tendency for components to occur earlier, especially in layer IV, whereas layer II appears to have a more lognormal distribution of component times. The most common times for components to occur in layer II was typically a few ms after the most common times for components to occur in layer IV, however this difference was less pronounced for the Basic stimulus than for the Contact stimulus.
By overlaying the stimulus waveforms on to the histograms of component times (Figure 3a), it is possible to qualitatively identify discrete modes or groups of components occurring together which may be representative of particular features. For instance, the component times in layer IV for Contact appear to have three separate modes: 9-18 ms representing the initial whisker deflection; 22-31 ms representing a second 'bump' in the stimulus; and 39-44 ms representing where the stimulus dips to its lowest point since the original whisker deflection. Further, the vast majority of the components in the second mode and all of the components in the third mode belong to units with two or three components, meaning that many units which respond to these later features of the stimulus also respond to earlier features. In comparison to layer IV, the component timings for Contact in layer II do not appear as well separated into distinct modes. Component timings in response to the Basic stimulus are more difficult to interpret, perhaps due to amplitude differences between the on-ramp velocity conditions at particular time-points. However, in the case of layer IV, all components occurring after 26 ms come from units which have previously responded to earlier stimulus features. Such components could therefore be signaling velocity or amplitude accumulation over slower on-ramps, or the beginning or continuation of the holding period. Components which occur after 35 ms could also be associated with the beginning of the off-ramp for the 400 mm/s on-ramp velocity condition.
These results gave a helpful qualitative overview of the responses component timings relative to stimulus features. The most notable of these results being that singlecomponent responses typically occurred earlier in the stimulus windows, that is, as "stimulus on" or "stimulus start" responses. Figure 3b shows the proportions of units with different numbers of response components to the Basic and Contact stimuli across layers and conditions. Across conditions, KS tests showed significant differences between the proportions across conditions in layer II but none in layer IV. Across layers, the only significant difference between the proportions was in the 4 days TBI condition for the Contact stimulus. Post-TBI, layer 2 showed a dramatic reduction in the relative proportion of responsive units for the Contact stimulus. Responses to the Basic stimulus generally show a relative increase in the proportion of single-component units in 4 days and 2 weeks TBI conditions across layers II and IV before normalization to sham levels in the 8-12 weeks TBI condition, where the number of multi-component units become approximately as numerous as single-component units but where the majority

F I G U R E 3 (a) Histogram of times of identified components occurring in units from sham animals for the Basic (two left panels) and
Contact (right two panels) stimuli. Component times are shown from units with one component (grey), two components (green), and three components (orange). Plotted in blue is the relative amplitude over time for the Basic and Contact stimuli. These amplitudes have been aligned with the x-axis to match what was being presented to the animal at the times these components are responding. In the case of the Basic stimulus, since five different on-ramp velocities were tested the five respective amplitude trajectories have been plotted in dashed lines. The blue solid line at the top of the Basic stimulus depiction is the hold period, and the beginning of one off-ramp period (from the 400 mm/sec on-ramp velocity variation) is also depicted with a dashed line. Using these plots, it is possible to approximate what particular components might be representing in the stimuli. For instance, the component times in layer IV for Contact appear to have three separate modes: 9-18 ms representing the initial whisker deflection; 22-31 ms representing a second "bump" in the stimulus; and 39-44 ms representing where the stimulus dips to its lowest point since the original whisker deflection. (b) Proportions of units with different numbers of components (non-responders = 0 components) for the Basic stimulus (left) and Contact stimulus (right) in layers II and IV with different health statuses. KS tests were conducted to determine if any proportions were significantly from one another. *, p < 0.05; ~ and &, p < 0.05; # and **, p < 0.001 The same differences were more pronounced when looking at the t tests performed using the mean firing rates on the components and the surrounding ±5 ms. While a vast majority of components had significantly higher mean firing rates than the surrounding ±5 ms at α = 0.05 (91% of Basic and 87% of Contact components), far fewer reached significance at α = 0.01 (46% of Basic and 52% of Contact components). The most significant differences were found for units with a single component, whereas few units with two or three components reach similar levels of significance; at α = 0. These results demonstrated the systematic changes in the proportions of unit response categories. At a high level, these changes show quite dramatic changes in early post-TBI conditions followed by a gradual return to proportions which were functionally similar to healthy conditions in the late post-TBI conditions.

| Temporal regularity of components occurring within the same unit
For both Basic and Contact stimuli, a number of significant relationships were found between the timings of first and second components in multi-component units. These relationships were particularly significant in the sham condition (p < 0.001 for both stimuli in both layers II and IV). Relationships in the 8-12 weeks TBI condition were also highly significant, especially in layer IV (p < 0.0001 for Basic and p = 0.0014 for Contact). However, there were irregular and insignificant timing relationships for both stimuli at the 4 days and 2 weeks TBI conditions (Figure 4).
To see how the relationships between these timings might reflect changes in the intracortical information flow after TBI, we calculated the y-intercepts shown in Table 2. We notice that in both layers and for both stimuli there is a general increase in the standard errors (SEs) of the y-intercepts, reflective of the higher scatter seen in Figure 4 for these relationships. Interestingly, considering the y-intercepts from just the sham and 8-12 weeks TBI conditions shows that where the intercepts in one layer increase, intercepts in other layer decreases-and the relationships are opposite in the two stimuli. This may indicate temporal modulation in thalamic input to cortex via layer IV which is then being compensated for in layer II. The y-intercepts (for those relationships which were significant), also gradually lowered across the TBI conditions-with one exception: the Contact layer IV yintercept increased after TBI. This exception is unusual given the Contact layer II y-intercept decreased, where we might expect it to follow the layer IV trend.
Like in the proportions of unit response categories (i.e., at the population level), these results at the individual unit level showed quite dramatic changes in early post-TBI conditions followed by a gradual return to component timings which were functionally similar to healthy conditions in the late post-TBI conditions.

| Relationship between graphtheoretic and independent neural response components analyses
Comparisons between the graph-theoretic and independent neural response components analyses were only possible for units with one component. These relationships, where significant, are all positive. In the case of layer II, for the sham and 4 days TBI conditions the mapping is close to 1:1 between the times obtained for the same components from the two analyses. However, at the 2 weeks and 8-12 weeks TBI conditions this 1:1 mapping breaks down. In the case of layer IV, there is a weak positive correlation between the times obtained from the two methods, but this remains significant and relatively constant across animal conditions. Overall, the median path times tend to occur earlier than the local maxima component times calculated for the same unit components. In summary, this suggests the graph-theoretic analysis is more sensitive to accurate timing information within unit response profiles whereas the neural response components analysis was more sensitive to precise overall timing of the response components. Figure 5 shows the non-stationarity of components (how much components "shifted" in time across different stimuli variations) for the Contact stimulus, separated by layer and condition. We can see that very few components are stationary, with most components moving across one or more ms as the stimulus varies in amplitude, the most common degree being 2 and 3 ms. This indicates that the function of microcircuitry controlling these temporal features are amplitude-dependent.

| Non-stationarity of components
Across the layer II conditions, most components move 2-3 ms, with a few being stationary in the sham condition and a minority across all conditions moving >3 ms. However, this minority of less stationary components seems to proportionally increase in TBI conditions compared to sham, although this trend in insignificant. What is significant are the of the differences between layer II and IV for the sham and 2 weeks TBI conditions. In the case of sham, F I G U R E 4 Linear regressions between the timings of first and second components of multi-component units responding to the Basic (left) and Contact (right) stimuli, separated by layer and condition. Red lines show line of best fit where the relationship is significantly correlated. L2 and L4 = layers II and IV; X.1 and X.2 = first and second components. Linear regressions between the timings of first and second components of multi-component units responding to the Contact stimulus, separated by layer and condition Basic Contact layer IV components are much less stationary than layer II components, for example, only nine components had a movement of 3 ms in layer II, whereas in layer IV there were 43 components. Conversely, for the 2 weeks TBI condition, the significant difference between layers II and IV is due to a proportional increase in more stationary components-in layer II there were eight components with a movement of 2 ms, whereas in layer IV there were 28 components. This shift towards more stationary components also made the layer IV 2 weeks TBI condition distribution significantly different to the layer IV sham condition, which proportionally had more non-stationary components. These results suggest that the amplitude-dependent microcircuit properties controlling these temporal features may cause these changes in layers II and IV. Mechanistically, we speculate that we find components are more "mobile" in layer IV than layer II due to a layer IV having more direct thalamocortical inputs than layer II or layer II having more interneurons which serve to lessen these amplitude-dependent effects. Most likely, it is some combination of these mechanisms. Like results for the proportion of unit response types and intra-unit timing relationships of components, these results further point to a significant change in early post-TBI conditions before a return to a distribution which functionally similar to healthy conditions in the late post-TBI condition. The exception in this case, however, is that differences appear slightly later, at the 2-week timepoint, but not at the 4-day timepoint.

| Relationships between components' mean firing rates
Across all layers and conditions there were many significant, positive correlations between the mean firing rates along component paths. Many of these paths are likely tracing segments of the same components, however it is also possible that these paths trace separate components which are affected by stimulus variations in a similar manner, for example, Figure 1d-i, bottom-right panel. Unfortunately, we do not find any significant (at α = 0.05), negative correlations, as we might expect from a unit like that shown in Figure 2b, left panel. However, the majority of paths were not significantly correlated, positively or negatively, and this may be due to the low statistical power resulting in a higher likelihood of type II errors. Where we do see weak negative correlations, we notice that the sham condition generally contains more units with such relationships than the TBI conditions (especially for the 4 days TBI and 2 weeks TBI conditions). In the 8-12 weeks TBI condition, these weak negative correlations appear more frequently and comparable to the sham condition in this respect. This suggests that, if such results are indicative of cells with components which have oppositely-weighted amplitude-dependence (e.g., Figure 2b, right panel), such cells are less frequent in the early-to-mid functional recovery stages of TBI. This indicates such cells may be susceptible to maladaptation as a result of TBI. In summary, these results further added to evidence of functional changes in early post-TBI conditions which become less pronounced or otherwise statistically similar (functionally) to health conditions in late post-TBI conditions.

| DISCUSSION
We have shown in a large database of online-sorted units from layers II and IV of barrel cortex, diffuse TBI causes layerdependent functional maladaptations which affect the timing and number of response components. Multi-component responsive units in healthy and late post-TBI conditions have reliable inter-component timing relationships, but T A B L E 2 Y-intercepts (±SEs) for linear regressions between the timings of first and second components of multicomponent units responding to the Basic and Contact stimuli in layers II and IV, as shown in Figure 4 not in early post-TBI conditions. We introduced a graphtheoretic method to characterize the shifting towards earlier post-stimulus onset responses with increasing stimulus intensity, which was more pronounced in layer IV than layer II in the healthy condition. Immediately after TBI, this trend reversed and at 8-12 weeks post-TBI no stimulus intensitydependence differences were observed between layers II and IV. We speculate this potential compensatory mechanism in response to injury may be caused by differences in thalamic input or local inhibitory neuronal makeup.

| Temporal activity patterns in healthy units
In healthy conditions, layer II units largely reflect the same sorts of temporal activity patterns as are seen in layer IV and in almost the same proportions. This is slightly surprising given the layers serve markedly different functions (Burns & Rajan, 2021)-the main layer IV excitatory output is to layers II and III whereas the layer II excitatory output is laterally to other parts of layer II (typically over several barrel columns), to other cells in layers III and V, and to secondary somatosensory and motor cortices (Aronoff et al., 2010;Feldmeyer et al., 2006). However, although the temporal activity patterns may be similar between layers II and IV, such patterns are not necessarily (and in fact unlikely to be) processed similarly. For example, the diverse inhibitory makeup of layer II may make transformation of input timing easier than complete silencing or splitting of arriving layer IV afferents, which would also suggest sensory processing in layer II barrel cortex is significantly temporal. For both the Basic and Contact stimuli there were highly significant (p < 0.001) relationships between the times of components within individual layers. However, the y-intercepts of the slope of the regression line between component timings for each layer were dissimilar, and the direction of difference between the layers was opposite between the stimuli (Table 2). This suggests that the way(s) layer II transforms temporal activity patterns from layer IV is stimuli-dependent. In both cases, the layer II y-intercept moves in the direction closer to ~10 ms relative to the layer IV y-intercept. This might reflect an optimal inter-component timing for layer II for processing of higher features, whereas layer IV, not needing to perform such processing and primarily acting to amplify the thalamic signal and distribute it to other cortical layers (Cowan & Stricker, 2004;Feldmeyer, 2012;Staiger et al., 2004), is capable of producing a wider variety of inter-component timings (perhaps courtesy of its computationally simpler task). The initial inter-component timings are also likely different for the two stimuli due to inherent differences in the stimuli themselves, for example the Basic stimulus had a constant peak amplitude (but variable on-ramp time) whereas the Contact stimulus had varied peak amplitudes (but constant on-ramp time). Notably, however, a large majority of the components occurred during the initial stimulus onset period, which appears to be a common response characteristic in other sensory systems (Nencini & Ivanusic, 2017). For both stimuli, the first components generally occur slightly later in layer II than in layer IV, which matches with the classical view of sensory information arriving at layer IV and then being transmitted up to layers II and III. Mechanistically, we could speculate that differences in timings between the components within individual unit components may be due to differences in inhibitory makeup between the layers, as well as the difference in thalamic input.

| Temporal activity patterns in TBI units
Proportions of unit response profiles and their changes post-TBI were quite different between the two stimuli: for the Basic stimulus there were more non-responsive units and almost no changes post-TBI whereas for the Contact stimulus there were initially fewer non-responsive units and many changes post-TBI. Although, in layer II, there was a general increase in non-responsive units in early post-TBI for both stimuli and then a general decrease in late post-TBI, which may be reflective of the tendency of more superficial layers to experience greater functional and anatomical damage in TBI. The differences between the two stimuli may reflect the more complex, naturalistic nature of the Contact stimulus (Hartmann et al., 2003). Common between both stimuli, however, was a general reduction in the number of multi-component units in the post-TBI conditions, even after 8-12 weeks. If multi-component units are responsible for processing more complex sensory information than fewer-or singlecomponent units, it is possible this reflects a permanent reduction in sensory processing capacity and associated cognitive deficiencies caused by TBI (Cantu et al., 2015;Ding et al., 2011;Johnstone et al., 2014).
Although the proportions of unit response profiles may have changed due to TBI, component times returned at 8-12 weeks post-TBI to reliable inter-component timing relationships like those seen in the healthy condition (see y-intercepts in Table 2). Nevertheless, these results suggest stimuli-dependent changes occurring in the timing between components, which agrees with the previously discussed sham timing data. This suggests either some considerable changes in the timing of thalamic efferent activity arriving in layers II and IV or TBI-induced (possible compensatory) mechanisms. Given the primary input to layer II is via layer IV, it seems reasonable to expect if this was purely a cortical compensatory effect, the y-intercept change should be immediate for layer IV at the 4 days TBI condition for both stimuli and would remain constant (whatever the direction of change), whereas the layer II y-intercepts would slowly shift, indicating compensation in response to perturbed layer IV input. Instead, we find the y-intercept is not constant across the TBI conditions in layer IV and instead slowly changes (along with layer II), suggesting a possible mixture of changes in thalamic input and intra-cortical compensatory mechanisms (and it is possible this mixture is integrated via thalamocortical interactions in deeper layers). Possible intra-cortical compensatory mechanisms may include gross anatomical changes such as the severing of feedforward inhibitory or excitatory synapses in layer II or slight but critical changes in spike timing leading to failure of previously expected feedforward signals (or, oppositely, success of previously unexpected feedforward signals).
Units responding to the Contact stimulus in layer IV of the 2 weeks TBI condition distribution had significantly (p = 0.008) fewer non-stationary components compared to the sham condition. A similar, but less significant (p = 0.035) relationship was also found for the Basic stimulus. This further suggests a change in layer IV thalamic input. And while we do not find significant differences among the layer II non-stationarity distributions for either stimulus, we do find significant differences for the Contact stimulus between layers II and IV for both the sham and 2 weeks TBI conditions. In the case of sham, layer IV components, they were much less stationary than layer II components, for example, only nine components had a movement of 3 ms in layer II, whereas in layer IV there were 43 components. Conversely, for the 2 weeks TBI condition, the significant difference between layers II and IV is due to a proportional increase in more stationary components-in layer II there were eight components with a movement of 2 ms, whereas in layer IV there were 28 components. This mis-match between the direction of change for component timings agrees with our earlier results regarding multi-component unit timings and their regressions' y-intercepts. It also supports the speculative idea that thalamic input in layer II, cortical compensatory mechanisms, or some combination thereof, is attempting to work "against" changes in layer IV activity, as driven by thalamic input. However, it is also possible that what we have been referring to as "compensatory" mechanisms (whether in cortex, via thalamic input, or corticothalamic interactions) are actually the malfunctioning of inhibitory microcircuits in layer II and such malfunctioning is the driver of these mis-matches. If so, taken together with our results regarding the activity component distributions, it seems likely that inhibitory malfunction in TBI disrupts the timing of components more than the overall proportions and presence of such components. The fact that the non-stationarity was more pronounced in layer IV than layer II for the health conditions (but not the late post-TBI conditions) indicates, however, a permanent and significant deviation from healthy conditions. Mechanistically, we speculate this change is likely caused by damage or maladaptation in inhibitory microcircuits. We suggest future studies probe such circuits in the context of temporal response profiles of neurons in TBI conditions.

| CONCLUSION
Collectively, along with the various stimulus-dependent and layer-dependent effects mentioned, these analyses suggest TBI disrupts normal temporal coding mechanisms in individual units of barrel cortex. Along with the previously-established general reduction in inhibition after TBI , this emphasizes the importance of fine-tuned inhibition on normal temporal coding mechanisms applied to sensory signals in layers II and IV. Simultaneously, it demonstrates that some unit profiles can appear functionally unchanged despite global excitatory-inhibitory imbalance and microcircuit changes. Commutatively, however, these effects lead to an increase in non-responsive units, changes in distributions of temporal activity patterns, and misalignment of the different components in multi-component temporal response patterns. While some of these changes appear transient, others appear to persist long after TBI.
It is likely some of these changes are partially responsible for the short-and long-term behavioural and cognitive changes which occur after TBI. This warrants future work to more precisely identify the specific neuronal sub-types involved in leading to these temporal activity differences, and also whether and where they arise (at different layers, and/or as transformations from thalamic activity). Future work may also consider tracking the behavioural or awake electrophysiological correlates of these temporal activity differences.