Conditioning by Subthreshold Synaptic Input Changes the Characteristic Firing Pattern of CA3 Hippocampal Neurons

The action potential discharge dynamics of neurons has conventionally been measured using somatic injection of step currents. Different morphological types of neurons are considered to have characteristic firing responses and therefore, these responses have been used for phenotypic classification. However, neuronal excitability is subject to homeostasis at both network and single neuron level, suggesting that discharge patterns should also be subject to change. In this study, we show that the firing patterns of neurons in the CA3 field of the rat hippocampus in vitro change significantly after only a few minutes of low frequency subthreshold stimulation of the neuron’s afferents. This effect, which was long-term, could be reproduced by subthreshold somatic depolarizing pulses and was blocked by kinase inhibitors, indicating that discharge dynamics are modulated locally. Cluster analysis of the firing patterns before and after conditioning revealed systematic transitions towards adapting and intrinsic burst behaviours, irrespective of the patterns initially exhibited by the cells. Using a conductance-based model and subsequent pharmacological blocking, we demonstrate that the experimental transitions can be mediated by a recruitment of calcium and M-type potassium conductances. We conclude that CA3 neurons adapt their conductance profile to the statistics of subthreshold activity in their embedding circuits, making their intrinsic firing pattern not a constant signature, but rather the reflection of their history of activity. Significance Statement Different anatomical types of neuron express a characteristic action potential discharge pattern in response to intra-somatic injections of step currents. Together with the cell’s morphology and molecular markers, these patterns have been used to classify neuronal phenotypes. However, we show that in the case of hippocampal CA3 neurons, this discharge pattern is not as characteristic as generally assumed. Instead, they modify their fundamental processing according to the subthreshold signals that they receive from their embedding circuit. This result implies that CA3 neurons collectively adapt their network processing, and also that their discharge response patterns cannot be used for phenotypic classification.

Three examples of neurons in the CA3 area presenting different morphologies and different firing patterns in control conditions. The discharge patterns were measured by injection of step currents of increasing amplitude. Control measurements (gray traces, left) were followed by stimulation of the mossy fibers. The upper trace shows all voltage traces elicited upon different levels of current injection on that cell. Two sample traces of this set are shown below. EPSPs (middle panel) were evoked in response to a stimulation with double current pulses, separated by 20 ms and repeated 500 times at 1 Hz. The series of repeated pulses are shown superimposed. A sample trace is highlighted in red. The inset shows the configuration of recording and stimulating electrodes (on the CA3 region of the hippocampus and on the dentate gyrus, respectively). Below, the morphology obtained by labeling the cells with biocytin is shown. After the conditioning, patterns were measured again (blue traces, right). A) Pyramidal cell switches from non-adapting burst to intrinsic burst firing. B) Pyramidal cell switches from delay accelerating to intrinsic burst continuous pattern. C) Bipolar cell switches from non-adapting continuous to adapting continuous firing (scale bars = 50µm). D) Mean Fraction of Spikes for the population in the first and second half of the voltage trace (green and yellow rectangle below the trace, respectively) for both control and conditioned cases. A significant redistribution on the fraction of spikes is observed after the conditioning, where the fraction of spikes on the first half is increased while it decreases in the second half (n=50, p=1.92e-6, two-sided Wilcoxon signed rank test). E) Empirical Cumulative Distribution Function for the data shown in D. Every individual cell, for both control and conditioned cases, is represented as the number of spikes for the first half of the trace minus the spikes for the second half (n=50) Whole-cell patch clamp recordings of CA3 neurons were performed in rat hippocampal 48 organotypic cultures. The intrinsic firing patterns of the neurons were recorded before 49 and after conditioning by extracellular stimulation of the mossy fibers originating in the 50 dentate gyrus. The conditioning stimuli consisted of paired pulses (0.1 ms duration 51 pulses, interval 10 -20 ms) applied at 1 Hz, and repeated 500 times for a total period of 52 approximately 8 minutes. The amplitude of the pulses was adjusted for each recorded 53 cell to elicit only subthreshold excitatory post-synaptic potentials (EPSPs). This mossy 54 fiber stimulation protocol is a modification of that described by Brandalise and 55 Gerber [27,28], which has been previously shown to elicit heterosynaptic subthreshold 56 plasticity in CA3 pyramidal-pyramidal synapses. The firing patterns of neurons were 57 assessed with a sequence of constant current injections. For convenience, we used the 58 terminology of the Petilla classification [15] to label these patterns. Interestingly, we 59 observed that after the conditioning protocol, the Petilla discharge label had to be 60 adapted for most of the cells, independently of their initial firing type. For example, the 61 pyramidal cell shown in Fig 1A had a non-adapting burst pattern before stimulation 62 (gray traces). After conditioning (blue traces), this response changed to intrinsic burst. 63 The same transition was observed for the pyramidal cell on panel 1B, whose initial 64 pattern was delayed accelerating. The bipolar cell on panel 1C switched from 65 non-adapting continuous to adapting continuous firing. We observed that the most 66 common transition performed by the cells was towards adapting and intrinsic burst 67 patterns. Indeed, the quantification of the mean fraction of spikes in the first half versus 68 the second half of the voltage for the population of recorded cells showed a distribution 69 of the spikes in favor of the first half (Figs 1D, 1E) (n=50). This result supports our 70 observations that the main pattern transitions are towards adapting and intrinsic burst 71 behaviors after the conditioning. These changes in firing pattern were present in most 72 cells immediately after the stimulation protocol, and were stable at least 15 minutes 73 after the stimulation. The mossy fiber conditioning was followed by a significant 36 MΩ 74 (25%) decrease in input resistance (Rin), (from 144.8 ± 73.0MΩ to 108.4 ± 65.3MΩ, 75 two-sided Wilcoxon signed rank test, p=1.1e-5). There was also a significant 5 mV (7%) 76 depolarization of the resting membrane potential (Vm) (-65.3 ± 5.0mV) with respect to 77 resting level (-70.4 ± 5.7mV, two-sided Wilcoxon signed rank test, p=2.3e-5, n = 50). 78 However, the firing pattern changes could not be induced neither by simply holding the 79 resting membrane potential at different values (see Fig S1,n = 10), nor by the 80 step-currents used to measure the discharge patterns (see Fig S1,  We attempted to resolve if synaptic input was necessary to elicit the changes, or 91 whether they could be induced directly at the soma. To this end, we used intra-somatic 92 injection of paired step current pulses whose parameters were chosen to elicit a similar 93 Fig 2. CA3 firing pattern transitions occur upon somatic conditioning and are blocked by kinase inhibitors. A) Example of an intrasomatic conditioned cell that switch from delay accelerating (gray traces) to intrinsic burst firing (blue traces). The conditioning protocol is shown in the middle column. EPSPs were evoked by injection of paired current steps, of 50 ms in duration and separated by 20 ms. The double steps were repeated 500 times at 1 Hz. The series of repeated pulses are shown superimposed. A sample trace is shown in red. B) Mean Fraction of Spikes for the population in the first and second half of the voltage trace for both control and conditioned cases. A significant redistribution on the fraction of spikes occurs after the conditioning. The fraction of spikes on the first half is increased while it decreases in the second half (n=12, p=0.0024, two-sided Wilcoxon signed rank test). C) Empirical Cumulative Distribution Function for the data shown in B. Every individual cell is represented as the number of spikes for the first half of the trace minus the spikes for the second half (n=12). D) Example of a mossy fiber conditioned cell (as described in Fig 1) under the presence of H-89 and Go 6983 (PKA and PKC inhibitors) on the intracellular pipette. The cell presents a delay accelerating pattern in control conditions and remains under such pattern after the conditioning protocol is applied. E) Mean Fraction of Spikes for the population in the first and second half of the voltage trace for both control and conditioned cases. The redistribution of the fraction of spikes was not significant after the conditioning (n=13, p=0.266, two-sided Wilcoxon signed rank test). F) Empirical Cumulative Distribution Function for the data shown in D. Every individual cell is represented as the number of spikes for the first half of the trace minus the spikes for the second half (n=13). somatic voltage response compared to that generated by the mossy fiber stimulation 94 (Fig 2). This direct subthreshold somatic stimulus evoked changes in discharge pattern 95 that were similar to those elicited by the indirect mossy stimulation. The cell in Fig 2A 96 displayed a delay accelerating firing pattern in control conditions and underwent a 97 transition towards intrinsic burst pattern after somatic conditioning. The population 98 data showed a significant redistribution in the fraction of spikes in favor of the first half 99 of the trace versus the second half after the conditioning (Figs 1B and C) (n=12). In 100 this result we observed the same tendency of neurons to become adapting and intrinsic 101 burst after conditioning. Furthermore, due to the nature of the conditioning at the 102 soma, this result also suggests that the mechanism inducing the firing pattern change is 103 not localized to synapses, but rather acts at a more central, probably somatic or 104 proximal dendritic level. We next sought to identify what internal mechanism could be 105 responsible for the firing pattern transitions. The firing pattern of the cell depends on 106 the distribution of membrane ion channels that the cell presents at its membrane [18]. 107 A possible mechanism would act upon this distribution. Due to the time scale of the 108 response (on the order of minutes) we ruled out protein synthesis of new channels on 109 the membrane. An alternative would be channel phosphorylation, a mechanism known 110 to affect the conductance on a relatively short timescale [29]. We reproduced the  We observed that the conditioning induced firing pattern changes from more regular 123 patterns towards early bursting and adapting patterns. We sought to quantify these 124 changes using hierarchical clustering methods [16,30,31] to establish more objectively 125 which discharge type to associate to every response, and to quantify the frequencies of 126 transitions between them. Previous studies have used clustering methods to quantify the 127 similarity between vectors of features extracted from the voltage traces, such as action 128 potential (AP) amplitude, firing rate, or accommodation index [16,30,31]. However, 129 those metrics are not suitable for of our dataset, because several features commonly 130 used in those methods are unaffected by the conditioning. For example, AP amplitude, 131 width and afterhyperpolarization (AHP) showed no difference before and after the 132 stimulation (AP amplitude: 78.63 ± 14.95mV, 75.60 ± 9.77mV, paired t-test, p=0.11, 133 AP half width: 1.11 ± 0.26ms, 1.10 ± 0.24ms, paired t-test, p=0.74, AHP: 13.62 ± 134 3.76mV, 12.66 ± 4.15mV, paired t-test, p=0.12, n = 50). Consequently, we chose to use 135 Dynamic Time Warping (DTW) as a comparison measure, because it operates directly 136 on the action potential sequence rather than relying on a pre-defined set of features (see 137 Methods for a detailed explanation). Feature vectors of the instantaneous firing rate of 138 the voltage traces were compared pairwise using the DTW algorithm. As an internal  Two main families can be identified: one containing adapting and bursting traces, together with delayed spiking patterns (left branch); and another branch containing regular and accelerating traces (right branch) (n=50). B) Representative traces from each cluster. Below, average instantaneous firing rate over all traces belonging to the same cluster. Middle lines indicate the mean; light outer lines indicate standard deviations. The instantaneous firing rate (in Hz) is normalized to 1. C) Transitions observed between firing patterns before and after conditioning. Each cell is assigned to a single cluster (represented as a box) for both the control and conditioned cases. Arrows indicate transitions between types whenever a cell changed cluster. Self-loops indicate that the firing pattern was retained after conditioning. Numbers indicate percentages of observed transitions, and the number of cells in each category under control conditions is displayed next to each pattern type. Cells tend to transition towards adapting and bursting patterns following conditioning (n = 43). Seven cells were assigned as unclassified.
interpret the hierarchy in terms of recognized response types [15]. Representative traces 144 of each family are shown in Fig 3B. The average of the firing rate vectors of every 145 cluster is depicted beneath each representative trace. The clustering algorithm captures 146 well the typical time courses of the firing patterns. The right branch of the cluster tree 147 contains accelerating and non-adapting firing patterns, while the other contains 148 adapting and intrinsic bursting patterns together with a smaller group of traces that 149 have delayed spiking profiles (Fig 3A). The consistency of the algorithm was confirmed 150 by its successful clustering of independent feature vectors derived from the same set of 151 current injections (same cell under the same conditions) into a single cluster. Indeed, in 152 86% of cases (43 of the 50 cells) the algorithm successfully allocated the majority of 153 vectors from the same set of current injections into single clusters. Vectors from the 7 154 remaining cells were not consistently classified. For 50% of the cells all of their voltage 155 traces fell into the same cluster, and for 90% of the cells at least 50% did (see Fig S3). 156 The allocation of some responses from the same cell into more than a single cluster does 157 however follow a biological logic. For example, for cells classified as accelerating, some 158 of their voltage traces could reasonably fall into the non-adapting cluster because 159 acceleration may vanish at high current injections. A similar reasonable misclassification 160 is possible for adapting traces. In this case low current injections may be classified as 161 non-adapting because the currents are not high enough to elicit adaptation (see Fig S4). 162 In particular, many of the traces belonging to the delayed spiking cluster come from 163 cells whose traces at low current injections were assigned to the accelerating cluster, or 164 belonged to non-adapting cells with spiking delay. The transitions between cluster types 165 induced by the stimulation protocol are shown in Fig   The consistent transition towards adapting and intrinsic bursting behaviors suggests a 183 common underlying mechanism for most cell types. Our results showing that 184 phosphorylation inhibition blocks firing pattern change after conditioning (Fig 2) 185 support the hypothesis that the prime candidate for this mechanism is a change in the 186 profile of active conductances contributing to action potential discharge dynamics. We 187 explored this possibility using simulations of action potential discharge in a 188 conductance-based single compartment neuron model containing 9 voltage and calcium 189 gated ion channels (see Methods). The densities and kinetics of these channels were 190 derived from experimental measurements of CA3 pyramidal neurons [12]. We tuned  reported in the literature [12]. In order to explain the experimental transitions, we 194 compared the performance of the clustering procedure on the model and the 195 experimental data. In a first step, the maximal conductance densities of the model were 196 tuned to match the various experimentally observed firing patterns. This tuning was 197 performed manually, and the match to the traces was qualitative. The absolute values 198 for the conductances required to match the main experimental categories (Fig 3: 199 adapting, intrinsic burst, delay spiking, accelerating and non-adapting) are reported in 200  Table S1. We also were able to reproduce the experimental traces in the 201 morphologically realistic model described by Hemond et al. [12] (see Fig S5). Although 202 the maximal conductance values had to be adjusted to satisfy the different impedance 203 of the more detailed morphology, the same key channels are responsible for each were swept through ranges that would likely encompass the experimentally observed 209 patterns (see Table S2 for the exact ranges). In this way a total of 861 conductance 210 profiles were generated. We obtained the discharge response to different levels of current 211 injection for each conductance profile, giving a total of 5166 voltage traces with their 212 associated conductance profiles. Every single experimental trace (coming from both, 213 control and conditioned cases) was matched against the collection of traces in the model 214 database using the DTW algorithm. The best fit was then selected, allowing us to 215 obtain an estimate of the conductance profile likely to be present in the experimental 216 neuron. These estimates also define the subset of model traces that best represent their 217 experimental counterparts. This subset was then fed to the same hierarchical clustering 218 procedure that was previously performed for the experimental data (Fig 3).  reported extensively, although the attention has been restricted primarily to the 261 modulation of firing rates for homeostatic plasticity [20,[32][33][34]. Regarding the dynamics 262 of the discharge, plasticity has been reported in lobster, with activity isolation being a 263 crucial component in shaping the patterns [35]. Modulation of the delay spiking pattern 264 in the hippocampus [36,37] and in the cortex [26] have been shown to be induced by This study was performed on organotypic cultures, derived from brain slices of 274 newborn rats that are incubated for three weeks using the roller-tube technique [40]. 275 Organotypic cultures have been used extensively to characterize electrophysiological 276 properties of hippocampal neurons and it is know that the tissue preserves the 277 anatomical organization of the adult hippocampus, as well as its connectivity and 278 characteristic spontaneous activity [41,42]. Most of the studies cited in this chapter 279 were done in cultures or juvenile acute brain slices, indicating that the plasticity of the 280 patterns is not unique to the organotypic preparation. It would be interesting to know 281 however whether this type of plasticity is also prominent in the adult brain and if it also 282 happens, at the same time scale, in other brain areas such as the cerebral cortex. Activity dependent changes of conductance have been extensively studied, and shown to 285 be triggered even by learning paradigms [39,43,44]. The work of Turrigiano et al. [35] 286 suggested that a calcium dependent mechanism could modulate the neural conductances 287 in STG lobster neurons, and that this would translate into changes in the cells' firing 288 patterns. Later work showed that depolarizing pulses at 1Hz could alter the density of 289 the calcium-dependent outward current ICaK and the transient outward current IA in 290 the STG [45]. These studies led to theories of homeostatic plasticity [20,33], which 291 propose that cells maintain both the turnover of ion channels, and a stable level of 292 activity, to compensate for changes in synaptic strength. However the time scale of such 293 mechanisms typically extends over hours, and presumably involves processes of gene 294 expression [46], whereas in our experiments the changes were observed immediately 295 after conditioning. Aizenman and Linden [47] observed rapid changes of excitability of 296 cerebellar cells after synaptic stimulation, and proposed a calcium-dependent 297 modification though phosphorylation of gCaT and gCaK to account for the observed 298 changes. Interestingly, these are the same candidate channels that we have identified as 299 underlying the discharge pattern changes in this study. Supporting these lines, rapid up-300 or down-regulation of ion channel conductance via phosphorylation or vesicle 301 modulation due to calcium signaling has been extendedly demonstrated [22,29,43] and 302 it has been shown that ion channels possess a complex of scaffold proteins containing 303 certain protein kinases that could selectively regulate channel conductance through 304 phosphorylation [29]. This mechanism could provide a link between the activity of the 305 network and the specific conductance recruitment. An alternative explanation to the 306 conductance recruitment is that continuous stimulation of the neuron may alter the ion 307 concentrations in the cellular environment; for example, by altering intracellular 308 potassium and calcium concentrations [48,49]. However, our simulations show that the 309 decay time constant of the intracellular calcium is too short to allow significant 310 accumulation over the period of conditioning (see Fig S6A-C). Even if the time constant 311 were greatly increased, the accumulation of calcium during conditioning would be 312 insufficient to elicit a significant change in firing pattern (see Fig S6D). Regarding 313 potassium, our extracellular concentration was less than that required [48] for the  gCaK. We are aware that alternative channels could elicit a similar dynamical response. 321 The effect on the spike delay mediated by a slow inactivating hyperpolarizing current, 322 such as gKd can also be elicited by a slow non-inactivating depolarizing current such as 323 gN ap. Thus, it is possible that different cells recruit different set of conductances 324 depending on their initial conductance profile. However, the candidates we propose have 325 been previously reported to shape the spiking response of the cell via activity dependent 326 mechanisms. For example, it is well established in the epilepsy literature that gCaT is 327 strongly associated with the switch to bursting mode in hippocampal cells [38,50] while 328 gKd in the hippocampus and similar potassium conductances in the cortex have been 329 shown to be up-or down-regulated according to network activity and modulate the 330 delay firing response of the cell [26,36,37]. Modulation of the M-type current upon 331 activity has also been shown in the hippocampal region CA1 [51] and in CA3 [52], with 332 the latter group reporting that transient subthreshold depolarizing pulses are more 333 effective in the modulation of the current.

334
The conditioning protocol elicited stereotypic transitions of pattern towards 335 adapting or intrinsic burst patterns. However, it was not equally likely for all cell types 336 to perform such transitions. For example, accelerating cells moved towards regular 337 patterns with higher probability than the rest of patterns (Fig 3C). We speculate that 338 either the initial density of channels favors the different likelihood of transitions, or that 339 a cell on such initial state must necessarily become regular during the transition to any 340 PLOS 12/31 . CC-BY-NC-ND 4.0 International license peer-reviewed) is the author/funder. It is made available under a The copyright holder for this preprint (which was not . http://dx.doi.org/10.1101/084152 doi: bioRxiv preprint first posted online Mar. 26, 2017; other pattern. An alternative is that there may be some cell types that obey distinct 341 rules. For example, we noticed that 4 cells from the non-adapting cluster had high firing 342 rates under control conditions (see Fig S7). Two of these had smooth cell morphologies. 343 The other two cells correspond to very densely spiny cells, with stellate morphologies. 344 Interestingly, although transitions towards bursting or classic adapting behaviors were 345 not observed on these cells, there was a modulation on the delay of the first spike in 346 both cell types, suggesting that the stimulation protocol had a differential effect on this 347 particular neural population.

348
One of the typical transitions that we observe in our dataset is the switch of cells 349 towards bursting behaviors. We emphasize that this is not the only transition that is 350 induced, but special attention should be given to the burst mechanism. It is known 351 from the literature that different types of cells can present this dual behavior. For 352 example, relay cells on the thalamus become bursty upon hyperpolarization because of 353 T-type conductance inactivation [53]. In our case, after the induction protocol, the cells 354 depolarized 5 mV in average, so we rule out this hyperpolarization mechanism. The 355 main form of discharge of CA3 cells have been known to be either regular or 356 bursting [12]. Although the firing pattern transitions were abolished in the presence of 357 PKA and PKC inhibitors, 2 cells out of 13 showed still transitions to intrinsic burst.

358
This could be likely due to failure of diffusion of inhibitors from the electrode, but we 359 cannot exclude a different mechanism for this type of transition (for example, through 360 different kinase pathways). 361

Functional implications of firing pattern modulation 362
The fact that neurons possess the internal machinery to mediate the observed transitions 363 raises questions about the computational consequences of such behavior. As proposed 364 by Shin et al. [54], a neuron that can dynamically adapt its output firing in response to 365 its input statistics would have important advantages. If such neuron could adjust its 366 threshold and dynamic range upon activity, it could respond to stimuli over a broad 367 range of amplitudes and frequencies without compromising the sensitivity and dynamic 368 range of the cell. Spike frequency accommodation has the characteristics of a high-pass 369 filter [55]. Since our conditioning stimuli occurred at constant frequencies, the cells may 370 have recruited a specific set of conductances that shift their integration properties so as 371 to gain sensitivity in the new spectrum range. Differences in filtering properties of brain 372 stem neurons have also been shown to facilitate the extraction of spatial information 373 from natural sounds [56] and most of the conductances that we identify in this study 374 have been shown to be frequency resonance candidates [57][58][59]. These resonance 375 properties of cells may have important functional implications for neural activity and 376 brain rhythms [60,61]. In addition, modeling studies have shown that a neuron able to 377 adapt to its own input statistics is able to maximize the mutual information between its 378 input and output firing rates [62]. This type of effect can emerge following firing rate 379 homeostasis rules and promote metaplasticity [63]; on the other hand it can be their 380 cause [64]. Finally, this fast adaptability of the firings may also be important for 381 specific memory acquisition on the hippocampus [39,65]. Further studies will be needed 382 in order to unravel the role that such firing pattern transitions may have for 383 computations in neural circuits. A first step towards this goal must be to explore more 384 generally how the form and frequency spectrum of somatic input signals on the long 385 time scale affect the distinct firing patterns that neurons exhibit on the short scale. We have shown that hippocampal neurons in rat organotypic cultures can rapidly adapt 388 their supratheshold action potential discharge patterns in response to subthreshold 389 paired pulse conditioning stimuli delivered to their somata either by activation of their 390 synapses, or directly by intrasomatic current injection. We propose that these changes 391 are mediated via phosphorylation by recruitment of calcium and M-type potassium 392 conductances, conditional on the statistics of their somatic input currents. Such a 393 mechanism would allow the neuron to adapt its output behavior to the requirements of 394 the network in which it is embedded. Our results also imply that the discharge   Current-voltage relationships were determined by step command potentials and had 417 duration of 1 s to ensure steady-state responses. Data were recorded using an Axopatch 418 200B amplifier (Molecular Devices). Series resistance was monitored regularly, and was 419 typically between 5 and 15 M Ω. Cells were excluded from further analysis if this value 420 changed by more than 20% during the recording. Junction potential and bridge was not 421 corrected.

422
Mossy fibers were stimulated with a bipolar tungsten electrode. The intensity of the 423 stimulus was adjusted to evoke subthreshold post-synaptic potential responses of 15 mV 424 on average in the recorded neuron (minimal stimulation + 20% stimulation intensity). 425 Action potential discharges were evoked by injected current steps (-0.08 up to 1.   Similarity distances between pairs of traces were calculated using the Dynamic Time 483 Warping (DTW) measure [66]. DTW takes into account that two similar signals can be 484 out of phase temporarily, and aligns them in a non-linear manner through dynamic 485 programming [67]. The algorithm takes two time series Q = q 1 , q 2 , . . . , q n and 486 C = c 1 , c 2 , . . . , c m and computes the best match between the sequences by finding the 487 path of indices that minimizes the total cumulative distance where w k is the cost of alignment associated with the k th element of a warping path 489 W . A warping path starts at q 1 and c 1 respectively, and finds a monotonically 490 increasing sequence of indices i k and j k , such that all elements q i in Q and c j in C are 491 visited at least once, and for the final step of the path i end = n and j end = m holds.

492
The optimal DTW distance is the cumulative distances y(i, j), corresponding to the 493 costs of the optimal warping path q 1 , . . . , q i and c 1 , . . . , c j . This distance can be 494 computed iteratively by dynamic programming: where d(q i , c j ) is the absolute difference between the elements of the sequence. The 496 optimal warping path is obtained by backtracking from the final element y(n, m), and 497 finding which of the three options (increasing i only, increasing j only, or increasing i 498 and j simultaneously) led to the optimal warping distance, until i = 1, j = 1 is reached. 499 A warping window constraint of 10% of the vector size was chosen [67].

500
The pairwise DTW distances were used to perform hierarchical clustering by Ward's 501 algorithm [68]. The number of classes increases with the level of the hierarchy. We 502 choose to cut the tree at a level that provided sufficient structure to interpret the 503 hierarchy in terms of recognized response types (for example, Ascoli et al. [15]).

504
Every recording for a given cell was treated as an independent observation, and 505 could in principle be assigned to any cluster. If the electrophysiological state of the cell 506 is expressed in all of its responses, then we expect that all the independent observations 507 derived from that cell should be assigned to the same cluster. However, traces derived 508 from current injections to the same cell in different conditions (pre-or post-stimulation) 509 are expected to be assigned to different clusters if there is significant change in the 510 underlying electrophysiological state.

511
In fact the independent traces did not cluster perfectly. Instead, the majority of 512 independent observations derived from a given state clustered together and there were a 513 few that fell into other clusters. Therefore, we chose to label the electrical type of each 514 cell according to the cluster that contained the mode of the traces for one set of current 515 injections. Cells for which no clear dominant cluster could be identified, e.g. because 516 half of the traces fell into one cluster, and half of them into another, were labeled as 517 unclassified. A cluster transition was recognized whenever the cell was assigned to 518 different clusters before and after the stimulation protocol.

519
The analysis was performed using custom-written software in MatlabR2011b. The  avoid spatial discretization problems in a single compartment [69,70]. The passive 527 properties associated with the model were obtained from Hemond et al. [12]. We set the 528 length and diameter of our compartment to 50 µm. The active properties were modeled 529 by including appropriate voltage and calcium gated ion channels whose density and 530 kinetics were obtained from experimental recordings performed in CA3 neurons [12].

531
The simulations were performed using NEURON [71]. We choose an integration step of 532 25 µs, which was approximately 1% of the shortest time constant in the model. The Each current I x is described by the equation whereḡ is the maximal conductance, m and h are activation and inactivation terms, 537 V is the membrane potential, and E the reversal potential of the channel. The reversal 538 potentials for N a+ and K+ were E N a = 50 mV and E K = -85 mV, respectively. The 539 equations describing the different channel kinetics (m, h) for every current were 540 obtained from Hemond et al. [12]. Following this reference, the three calcium 541 conductances (T, N and L) were incorporated into a single parameter g Ca .

542
The set of maximal conductance values that are consistent with all our 543 experimentally observed firing patterns are shown in the Fig S1. The intracellular 544 calcium dynamics were modeled [12], as follows: The first term of the above equation describes the change caused by Ca 2+ influx 546 into a compartment with volume v. F is the Faraday constant, I Ca is the calcium 547 current and τ Ca is the time constant of Ca 2+ diffusion.

548
The occasional decrease in spike amplitude seen in some of the experimental traces is 549 probably due to sodium inactivation. We choose not to include this feature in the 550 model, because it does not affect the overall dynamics of the spike discharge itself.     Protocol applied to the model cell: 1 Hz current stimulation by double current pulses that elicited a depolarization of 10 mV, repeated 500 times. B) Comparison of model pulses with those elicited in the soma of experimental cells. C) Due to kinetics of calcium decay, the ion does not accumulate over period of stimulation (black trace). Decay must be much longer for calcium to accumulate significantly (green trace). D) Hypothetical increase in intracellular increase has little effect on pattern of discharge, even when increased 1000 fold (from left to right).