Experimentally-constrained biophysical models of tonic and burst firing modes in thalamocortical neurons

Somatosensory thalamocortical (TC) neurons from the ventrobasal (VB) thalamus are central components in the flow of sensory information between the periphery and the cerebral cortex, and participate in the dynamic regulation of thalamocortical states including wakefulness and sleep. This property is reflected at the cellular level by the ability to generate action potentials in two distinct firing modes, called tonic firing and low-threshold bursting. Although the general properties of TC neurons are known, we still lack a detailed characterization of their morphological and electrical properties in the VB thalamus. The aim of this study was to build biophysically-detailed models of VB TC neurons explicitly constrained with experimental data from rats. We recorded the electrical activity of VB neurons (N = 49) and reconstructed morphologies in 3D (N = 50) by applying standardized protocols. After identifying distinct electrical types, we used a multi-objective optimization to fit single neuron electrical models (e-models), which yielded multiple solutions consistent with the experimental data. The models were tested for generalization using electrical stimuli and neuron morphologies not used during fitting. A local sensitivity analysis revealed that the e-models are robust to small parameter changes and that all the parameters were constrained by one or more features. The e-models, when tested in combination with different morphologies, showed that the electrical behavior is substantially preserved when changing dendritic structure and that the e-models were not overfit to a specific morphology. The models and their analysis show that automatic parameter search can be applied to capture complex firing behavior, such as co-existence of tonic firing and low-threshold bursting over a wide range of parameter sets and in combination with different neuron morphologies.

Introduction Thalamocortical (TC) neurons are one of the main components of the thalamus and have been extensively studied in vitro and in computo, especially in first order thalamic nuclei in different species [1]. One of these nuclei, namely the ventral posterolateral nucleus (VPL), relays somatosensory, proprioceptive, and nociceptive information from the whole body to the somatosensory (non-barrel) cortex [2]. The VPL is located close to ventral posteromedial nucleus (VPM), which transmits information from the face to the barrel cortex. The VPL and VPM nuclei constitute the ventrobasal (VB) complex of the thalamus [3].
Despite its key role in sensory functions, a systematic characterization of the cellular properties of the VB complex is still missing. The morphologies of VPL neurons in adult rats were described in early anatomical studies but were limited to two-dimensional drawings of Golgiimpregnated cells [4]. The general electrical properties of TC neurons maintained in vitro are known and similar in different thalamic nuclei and species with respect to the generation of two distinct firing modes, called tonic firing and low-threshold bursting [5][6][7][8]. However, a systematic description on the electrical types in the VB thalamus in the rodents is still missing.
Collecting morphological and electrophysiological data, by following standardized experimental procedures, is essential for the definition of cells types and it is the first step to constrain computational models of single neurons [9,10]. Although models of TC neurons have been published previously, they were typically aimed at studying specific firing properties and their parameters were hand tuned to achieve the desired result [11][12][13][14][15].
The purpose of our study is to systematically define the morphological and electrical types by collecting in vitro experimental data and to constrain biophysically detailed models of VB TC neurons of the juvenile rat. To the best of our knowledge, automatic parameter search has not been applied, thus far, to capture complex firing behavior in thalamic neurons, in particular low-threshold bursting and tonic firing. We defined the electrical and morphological types of TC neurons through in vitro patch-clamp recordings and 3D morphological reconstructions. We then extended an existing method [16] to account for their distinctive firing properties. These electrical models (e-models) were constrained by the electrical features extracted from experimental data [9,10,17,18]. Other experimental data were used to assess the generalization of the models to different stimuli and morphologies. We further performed a sensitivity analysis by varying each parameter at a time by a small amount and recording the resulting electrical features. This analysis provided an assessment of the robustness of the models and a verification that the selected features provided sufficient constraints for the parameters.

Physiological and morphological characterization
We characterized TC neurons in slices of the rat VB thalamus, by combining whole-cell patchclamp recordings, biocytin filling and 3D Neurolucida (MicroBrightField) reconstruction, along with anatomical localization in a reference atlas [19] (Fig 1).
Visual inspection of 50 reconstructed morphologies (24 from the VPL, 26 from the VPM nuclei) revealed variability in the number of principal dendritic trunks and their orientation, in agreement with previous anatomical studies [4]. The maximum radial extent of the dendrites ranged between 120 and 200 μm and they started to branch between 20 and 50 μm from the soma (S1 Fig). We then analyzed the morphologies with two methods in order to quantitively classify different morphological types. We used algebraic topology to extract the persistent homology of each morphology and to visualize the persistence barcode [20] (Fig 2A, see Methods). Each horizontal bar in the persistence barcode represents the start and end point of each dendritic component in terms of its radial distance from the soma. The barcodes of all Third row is a low-threshold burst response from a hyperpolarized holding potential, V hold = −84 mV (burst mode), the other responses are elicited from a depolarized holding potential, V hold = −64 mV (tonic mode). Two different holding currents (I hold -tonic, I hold -burst) are injected to obtain the desired V hold . The vertical scale bar applies to all the traces, the first horizontal scale bar from the top refers to the first two rows, the second applies to the last four rows. (C) Analysis of adaptation index (AI) from recordings in tonic mode. Solid line is a non-parametric estimation of the distribution, dashed lines are two the morphologies followed a semi-continuous distribution of decreasing length. To quantify the differences between the barcodes, we computed the pairwise distances of the persistence images (see Methods and S1 Fig). We found that they were in general small (<0.4, values expected to vary between 0 and 1). These findings indicate that the morphologies cannot be grouped in different classes based on the topology of their dendrites. Furthermore, we performed Sholl Analysis [21] to compare the complexity of the dendritic trees ( Fig 2B). We observed that all the morphologies had dense dendritic branches, with a maximum number of 50-100 intersections between 50-80 μm from the soma. When comparing the Sholl profiles for each pair of neurons we could not find any statistically significant difference (S1C Fig). Considering the results of topological and Sholl analyses, we grouped all the morphologies in one morphological type (m-type) called thalamocortical (TC) m-type.  We used an adaptive stimulation protocol, called e-code, consisting of a battery of current stimuli (see Methods for details), where the stimulation amplitude was adapted to the excitability of different neurons. This standardized protocol has previously been used to build biophysically-accurate models of cortical electrical types (e-types) [16]. However, TC neurons from different thalamic nuclei and species fire action potentials in two distinct firing modes, namely tonic firing, when stimulated from a relatively depolarized membrane potential or low-threshold bursting, from a hyperpolarized membrane potential [5][6][7][8]. We thus extended the e-code to include two different holding currents. All the neurons recorded in this study displayed tonic and burst firing, when stimulated with the appropriate holding current (Fig 1B). Moreover, we were able to classify different e-types by considering the voltage traces recorded in tonic mode in response to step current injections (Fig 1B). The majority of the cells (59.3%) showed a non-adapting tonic discharge (continuous non-adapting low-threshold bursting, cNAD_ltb e-type) while others (40.7%) had higher adaptation rates (continuous adapting lowthreshold bursting, cAD_ltb e-type), as reflected by the adaptation index ( Fig 1C). We followed the Petilla convention [22] for naming the tonic firing discharge (cNAD or cAD), extending it to include "_ltb" for the low-threshold bursting property. In some rare examples, we noticed acceleration in the firing rate with decreasing inter-spike intervals (ISIs) towards the end of the stimulus. Similar adapting and accelerating responses have already been described in the VB thalamus of the cat [7]. We also observed stereotypical burst firing responses within the same cell, with variation of the number of spikes per burst in different cells, but the burst firing responses alone were insufficient to classify distinct e-types.

Constraining the models with experimental data
Multi-compartmental models comes with the need of tuning a large number of parameters [23], therefore we constrained the models as much as possible with experimental data. We first combined the morphology and the ionic currents models in the different morphological compartments (soma, dendrites and axon). Given that the reconstruction of the axon was limited, we replaced it with a stub representing the initial segment [16]. We used previously published ionic current models and selected those that best matched properties measured in rat TC neurons (see Methods). The kinetics parameters were not part of the free parameters of the models.
The distribution of the different ionic currents and their conductances in the dendrites of TC neurons is largely unknown. The current amplitudes of the fast sodium, persistent and transient (A-type) potassium currents were measured, but only up to 40-50 μm from the soma [24]. Indirect measures of burst properties [15] or Ca 2+ imaging studies [25] suggest that the low-threshold calcium (T-type) channels are uniformly distributed in the somatodendritic compartments. We thus assumed different peak conductance in the soma, dendrites and axon for all the ionic currents, except for I CaT , which had the same conductance value in the soma and dendrites. We then extracted the mean and standard deviation (STD) of different electrical features in order to capture the variability of firing responses from different cells of the same etype [9,10,17] (Fig 3). We observed that some features extracted from tonic firing responses had distinct distributions between the cAD_ltb and cNAD_ltb e-types ( Fig 3A).
For optimizing the models' parameters, we chose features that quantified passive (input resistance, resting membrane potential), burst and tonic firing properties (number of spikes, inverse of inter-spike intervals, latency to first spike, adaptation index), action potentials shape (amplitude, half-width, depth of the fast after-hyperpolarization). We aimed at finding the minimal set of features that captured the most important properties in the two firing modes. This set was a trade-off between comprehensively describing the experimental data (i.e. extracting all possible features), which can lead to over-fitting and loss of generalizability, and a too small set that would miss some important characteristics. For the tonic firing responses, we used three stimulation amplitudes (150%, 200%, 250% of firing threshold) which have been shown to reproduce the complete input-output function of the neurons [16,17]. Responses to two hyperpolarizing steps of different amplitudes (−40% and −140% threshold) constrained the input resistance (conductance of the leak current) and the conductance of currents activated in hyperpolarization, for example the h-current, I h (sag_amplitude feature). We included baseline voltage values to the optimization objectives to ensure that the models were in the right firing regime and spike count to penalize models that were firing in response to the holding currents. To constrain the low-threshold burst we used features (such as number of spikes) which are influenced by specific ionic currents, for example the low-threshold calcium current, I CaT .
The average value and STD of each feature were used to calculate the feature errors ( Fig  4C). Each error measured how much the features of the models deviated from the experimental mean, in units of the experimental STD. We used a multiobjective optimization approach (MOO), where each error was considered in parallel. To rank the resulting models after optimization, we considered model A better than model B if the maximum error of all the features of A was smaller than the maximum error of all the features of B.
By applying this MOO procedure, we generated multiple models with distinct parameter combinations for each of the twenty-two free parameters ( Fig 4B). The free parameters were allowed to vary between the upper and lower bounds shown in Fig 5B. The models reproduced well the key firing dynamics observed in the experimental recordings. They showed a lowthreshold burst when stimulated from a hyperpolarized membrane potential, crowned by a variable number of sodium spikes (Fig 4B). In the tonic firing regime, they reproduced adapting and non-adapting firing discharges as observed in the two e-types. These results indicate that the ion channels included in the models were sufficient to reproduce the experimental Experimentally-constrained models of thalamocortical neurons firing properties and that different e-types in TC neurons could be generated by different ion channel densities.

Model and experimental diversity
We found that different sets of parameter values reproduced the target firing behavior ( Fig  5B). We further analyzed models that had all the feature errors below 3 STD. Models' voltage responses reflected the characteristic firing properties of TC neurons (S3 Fig), indicating that the selected set of features and ion channels were sufficient to capture the two firing modes, in both the adapting and non-adapting e-types. The voltage traces from different models showed  small differences in spike amplitude, firing frequency, and depth of the after-hyperpolarization, as reflected by the variability of features values (Fig 5C), arising from differences in ion channel densities between models.
Spike-shape related features (e.g. AP. amplitude) in the different models covered the space of the experimental variability, while for some features (e.g. input resistance, R input ), all models tended to cluster on one of the tails of the experimental distribution. R input relates to the neuron passive properties and depends both on the number of channels open at rest (inverse of the leak conductance in the model) and the size of the cell. Given that all the models for a given e-type were constrained on a single morphology, this result is not surprising. Other features, such as sag amplitude were less variable in the models compared to experiments. We hypothesized that this depended on the variable stimulation amplitudes applied to different experimental cells, while all the models were stimulated with the same current amplitudes.
Some other features were systematically above or below the experimental values in both etypes. We suggest that this depend on the exact dynamics of some specific ion channels. For example, the amplitudes of the first and second spikes in the burst tended to be similar or above and below the experimental values, respectively. This can depend on the specific activation/inactivation properties of some ionic currents, for example the transient sodium current (I NaT ) and delayed potassium current (I Kd ). During the rising phase of the low-threshold spike, I NaT in the model is readily activated and generated a first spike with higher amplitude, but is not repolarized enough by the activation of I Kd . At higher potentials, reached towards the peak of the low-threshold spike, the availability of I NaT and other depolarizing currents seemed reduced and generated a spike with smaller amplitude. Sensitivity analysis confirmed that I NaT and I Kd had an influence on the amplitude of the first and second spike in the burst. Furthermore, these two currents operate together with currents that generate the burst, such as the low-threshold calcium current (I CaT ) and the I h in shaping the amplitude of the second spike in the burst. Interestingly, the models also tended to have lower instantaneous frequency of the first two spikes in the burst (Inv. 1 st ISI) and this feature had similar sensitivity (but of opposite signs) to the amplitude of the second spike in the burst.
Another possible explanation is the lack of some ionic currents in the model, for example some specific subtype of potassium channels that promote higher firing rates (Kv3.1 and Kv3.3). While neurons of the thalamic reticular nucleus are known to express this channel subunit [26], the expression in TC neurons has not been confirmed yet. The dynamics of I Kd could also explain why the after-hyperpolarization (AHP depth) tended to be smaller in the models compared to the experimental values. AHP depth is also influenced by other ionic currents, such as high-threshold calcium current (I CaL ), calcium-activated potassium current (I SK ) and the intracellular calcium dynamics. The number of action potentials (Num. of APs) in different conditions (No stim, I hold ) ensured that the models did not spike in the absence of a stimulus or in response to the holding current. We examined the diversity of the parameter values with respect to the initial parameter range (Fig 5B). Most of the optimized parameter values spanned intervals larger than one order of magnitude. On the other hand, some parameter values were restricted to one order of magnitude, for example the permeability of the low-threshold calcium current P CaT . This result is in agreement with experiments showing a minimum value of I CaT is critical to generate burst activity and this critical value is reached only at a certain postnatal age [27]. The value of P CaT was constrained by features measuring burst activity (such as number of spikes, frequency, etc.).

Assessment of model generalization
We used different stimuli for model fitting (current steps) and for generalization assessment (current ramps and noise). We simulated the experimental ramp currents, by stimulating the models with the appropriate holding currents for the two firing modes and a linearly increasing current. We first compared visually the model responses with the experimental recordings ( Fig 6A). In burst mode, the models reproduced the different behaviors observed experimentally: absence of a burst, small low-threshold spike, burst (S4A Fig). Moreover, the latency of burst generation substantially overlapped with the experimental one. However, a small fraction of models (1.2%) generate repetitive burst that we have never observed in the experimental recordings (S4B Fig). These models were quantitatively rejected by considering the number of spikes and the inter-spike intervals. In tonic mode, the latency to first spike, the voltage threshold, the shape of the subsequent action potentials and the increase in firing frequency were comparable with the experimental recordings ( Fig 6A). In addition, we quantified the Experimentally-constrained models of thalamocortical neurons generalization error to ramp stimuli (Fig 6C), by considering the latency to first spike, firing frequency increase over time (tonic mode) or number of spikes (burst mode).
Although conductance-based models can be fit by using step and ramp currents, these stimuli are different from synaptic inputs, which can be simulated by injecting noisy currents. To test the response to such network-like input, we used a noisy current varying accordingly to an Ornstein-Uhlenbeck (OU) process [28] to compare models' responses with the experimental data. Each experimentally recorded cell was stimulated with the same OU input, scaled by a factor w. Experimentally, w was calculated by evaluating the responses to previous stimuli. We developed a similar approach to generate the noise stimuli in silico (see Methods). The noise current was injected on top of the holding currents used during the optimization. We found that the models reproduced well the subthreshold potential, spike times and the distribution of single spikes and bursts (Fig 6B). Moreover, we quantitatively evaluated the generalization to the noise stimulus by extracting features (e.g. number of spikes) and comparing them with the experimental mean.
We computed generalization errors for each model, which were calculated similarly to the optimization errors (Fig 6C). We considered a model acceptable after generalization if it had all generalization errors <3 STD and we found that the majority of the models (>90%) passed the generalization test (Fig 6D).

Sensitivity of electrical features to small parameter perturbations
We assessed the robustness of the models to small changes in their parameter values. To that end, we varied each parameter at a time by a small amount (± 2.5% of the optimized value) and computed the values of the features. A sensitivity value of 2 between parameter p and feature y means that a 3% change in p caused a 6% change in f. We ranked the parameters from the most to the least influential and the features from the most sensitive to the least sensitive.
Some features resulted to be more sensitive to parameter changes, both in term of magnitude of the sensitivity and number of parameters (e.g. adaptation index, inverse of inter-spike intervals, ISIs, AHP depth). Most of these features describe the model firing pattern, which depend more on the interplay between the different ionic currents than on the specific activation/inactivation dynamics. Conversely, spike shape-related features were less sensitive to parameter changes (e.g. AP half-width, AP amp.) because they depend more on specific ionic current dynamics (e.g. I Kd , I L , I NaT, ). Some features were very weakly influenced by small parameter changes, e.g. baseline voltage, which depend more on the holding current amplitude, than on the model parameters ( Fig 7A).
The conductance of the leak current g leak emerged as the most influential parameter ( Fig  7A). An increase in g leak caused a decrease in firing frequency (inverse of ISIs) in both the tonic and burst firing modes. These results are easy to interpret when considering Ohm's law: increasing g leak means decreasing the input resistance of the model, so that for the same input current the voltage response becomes smaller. The second most influential parameter was the conductance of the persistent sodium current g NaP in the dendrites, which increased the tonic firing rate as expected from a depolarizing current. Interestingly, g NaP had an effect on the late phase of the low-threshold burst (inverse last ISI-burst), suggesting that the low-threshold burst is initiated by the activation of I T and modulated by I NaP . An increase in the permeability of the low-threshold calcium current P CaT , known to be one of the main currents underlying low threshold bursting, enhanced burst firing responses (it increased the inverse of ISIs, i.e. the firing frequency) and had effects on some of the tonic features. Increasing the somatic permeability of the high threshold calcium current P CaL decreased the tonic firing rate, despite being a depolarizing current. Increasing P CaL means higher Ca 2+ influx and higher activation of the Ca 2+activated potassium current (I SK ). The parameter g SK had indeed a similar effect on the features and thus clustered together with parameters regulating the intracellular calcium dynamics γ Ca and τ Ca (Fig 7B). Sag amplitude, that is known to depend on the activity of I h , was mainly influenced by change in g leak , P CaT and g h . In summary, each parameter influenced at least one feature. These results indicate that the model ability to generate tonic and burst firing is robust to small changes in parameter values and that all the parameters were constrained during the optimization by one or more features.
We then analyzed which features depended similarly on parameter changes, as they may add superfluous degrees of freedom during parameters search. Fig 7B shows the same sensitivities as in Fig 7A, clustered by their similarities (see Methods). Features clustered together if they were sensitive to similar parameter combinations and parameters clustered based on their similar influence on the features. Not surprisingly, the same tonic features measured at different level of current stimulation clustered together (e.g. AP amplitude and half-width, AHP depth, latency of the first ISI) and tonic firing features belonged to a cluster that was different from burst features. Some features measured in tonic mode (such as AP half-width and AP amp.) clustered together because they depended mainly on the dynamics of I NaT and I Kd : increasing the conductance of I NaT increased the amplitude of the APs and decreased its duration. This was also true for the amplitude of the 1 st AP in the burst. Features measured in burst Experimentally-constrained models of thalamocortical neurons mode had similar sensitivities because they depend on currents that are active at relatively hyperpolarized potential (such I H and I CaT ).

Preservation of model firing properties with different morphologies
We optimized the parameters for the adapting and non-adapting e-models in combination with two different experimental morphologies selected at random and then tested them with the other 48 morphologies. Considering that morphologies could not be classified in different m-types based on topological analysis of their dendrites and that TC neurons have been shown to be electrically compact [15], we expected the electrical behavior to be conserved when changing morphology. Nonetheless, different neurons vary in their input resistance R input and rheobase current I thr due to variation in the surface area. Variation in R input and I thr made the current amplitude applied during the optimization inadequate to generate the appropriate voltage trajectories. We thus devised an algorithm to search for the holding current to obtain the target holding voltage (for example −64 mV or −84 mV for tonic and burst firing, respectively) and I thr from the desired holding voltage. The different e-model/morphology combinations (me-combinations) were evaluated by computing the same feature errors calculated during optimization (Fig 8A). For each morphology, we selected the e-model that generated the smallest maximum error. We chose the value of 3 STD as the threshold to define which me-combinations were acceptable [29], yielding 48 acceptable me-combinations out of the 48 tested ( Fig 8A). All me-combinations reproduced burst and tonic firing (Fig 8B).
Given that the generalization of the electrical models to the other 48 morphologies worked well, we can conclude that the morphological properties of the modeled neurons are very similar, at least for properties that have an impact on the electrical models (e.g. surface area, diameters of the compartments).

Discussion
Our objective was to apply and extend an existing data-driven pipeline to identify the cell types and build models of VB thalamocortical neurons that reproduced the multiple firing modes that have been experimentally observed. We successfully modelled these novel firing types, by including additional stimulation protocols and features to constrain the low-threshold burst.
Our morphological and electrical data were used to define the properties of VB TC neurons in the rat. We found two electrical types (e-types) of TC neurons, but no objectively different morphological types (m-types) were revealed either using Sholl analysis [21] or topological analysis of dendritic branching [20]. We cannot exclude that refinements to these methods will reveal different m-types similar to the ones described in the visual thalamus of the mouse [30]. We also showed that automatic parameter search can be applied to build biophysically and morphologically detailed models. This method was already applied to model canonical firing behavior in cortical [9,10,16,17], hippocampal [31], cerebellar granule neurons [32] and corticospinal neurons [33]. To the best of our knowledge, such an automatic parameter search has not been used previously to capture different firing modes and complex firing behavior such as low-threshold bursting in thalamic neurons. Standardized electrophysiological protocols allowed us to identify for the first time in juvenile rat adapting and non-adapting e-types of TC VB neurons that were previously observed in other species [7]. This finding suggests that the intrinsic properties of TC neurons contribute to adaptation, a key phenomenon for filtering out irrelevant stimuli, before sensory information reaches the neocortex. Further experiments are needed to elucidate the relative contribution of intrinsic mechanisms and network properties to adaptation in somatosensory systems. We named the two main e-types continuous non-adapting low-threshold bursting (cNAD_ltb) and continuous adapting lowthreshold bursting (cAD_ltb) by following and extending existing conventions [16,22,31].
In this study, we improved upon previous morphologically and biophysically detailed models of tonic and burst firing in TC neurons [12,13,15] by explicitly constraining the parameters with experimental data, without hand-tuning of their values. Unlike previous models, we chose a multi-objective optimization for a methodological and a scientific reason: it is more time-efficient, reproducible, and it approximates the variability in ionic channel expression of biological neurons [31,[34][35][36], as shown by the family of acceptable solutions we found. However, experiments aimed at quantifying ion channel conductances are essential to assess if these solutions fall between biological ranges. Furthermore, we tested the generalization capability of the models and found that more than 90% of the models were comparable with the experimental data.
Nonetheless, we noticed some inaccuracies when comparing the voltage traces with the experimental data when assessing the generalization of the models. For instance, some models tended to generate small transient oscillations in response to ramp stimuli in burst mode. This result is not surprising, considering that the exact kinetics for all the ionic currents are not available and that there are known limitations in models of ionic channels derived from the literature or from other models [37,38]. In particular, modifications to the kinetics of the lowthreshold calcium current was shown to explain the propensity to generate oscillatory bursts in TC neurons of other nuclei and species [39]. More generally, we included ion channels that were used in previous models and that were validated with experimental data whenever possible. We undertook an extensive literature review to use channel kinetics derived from recordings in rat TC neurons from the ventrobasal (VB) thalamus or other first-order thalamic nuclei, whenever the data was available (see Methods). Moreover, we cannot exclude that some ionic currents were missing from our models and that they could have improved their fitness.
TC neurons have been shown to be electrically compact [15] and could, in principle, be modeled as a single compartment. However, active mechanisms need to be located in the dendrites in order to ensure synaptic integration and amplification [40]. Information regarding specific conductances or firing properties in the dendrites of TC neurons is limited. For this reason, dendritic parameters in our models may be underconstrained. However, the sensitivity analysis revealed that dendritic parameters did not appear to be the least constrained because they influenced different tonic and burst-related features.
We included in the model fitting and validation pipeline a sensitivity analysis, which is often neglected in computational neuroscience [41]. Although we cannot use our simple univariate approach to explore multidimensional parameter correlations and principles of co-regulation of ion channels expression, it is useful to find better constraints for parameters optimization. The selection of the features is indeed a step that still requires care and experience by modelers. Furthermore, this type of sensitivity analysis allows to identify parameters that can be traded-off during the optimization and that can be removed in order to reduce the dimensionality of the problem. In our study, parameters related to the calcium dynamics were shown to influence the features in a very similar fashion. This type of analysis is of particular importance in future work aimed at using the full diversity of ion channels that can be inferred from gene expression data. Gene expression data could also provide additional constraints on the choice of ion channels and indicate the ones that are missing in our models. More in detail, we propose that sensitivity analysis should be a fundamental tool in selecting which conductances are successfully optimized by the available experimental constraints. The example we showed is a local approach, applied to a specific solution to the optimization problem, which showed that our models are robust to small parameter changes. This analysis can be extended to study how the sensitivities vary in the neighborhood of different solutions.
In conclusion, we systematically studied the morphological and electrical properties of VB TC neurons and used these experimental data to constrain single neuron models, test their generalization capability and assess their robustness. Further work will validate these models in response to synaptic activity, in order to include them in a large-scale model of thalamocortical microcircuitry [16].

Experimental procedures
Experimental data were collected in conformity with the Swiss Welfare Act and the Swiss National Institutional Guidelines on Animal Experimentation for the ethical use of animals. The Swiss Cantonal Veterinary Office approved the project following an ethical review by the State Committee for Animal Experimentation.
Membrane potentials were sampled at 10 kHz using an ITC-18 digitizing board (Instru-TECH, USA) controlled by custom-written software operating within IGOR Pro (Wavemetrics, USA). Voltage signals were low-pass filtered (Bessel, 10 kHz) and corrected after acquisition for the liquid junction potential (LJP) of −14 mV. Only cells with a series resistance <25 MO were used.
After reaching the whole-cell configuration, a battery of current stimuli was injected into the cells and repeated 2-4 times (e-code). During the entire protocol, we defined offset currents in order to keep the cell at −50 mV (tonic firing) or −70 mV (burst firing) before LJP correction and applied them during the entire protocol. The step and ramp currents were injected with a delay of 250 ms in the experiment. In the models, the stimuli were injected with a delay of 800 ms, to allow for the decay of transients due to initialization. Each stimulus was normalized to the rheobase current of each cell, calculated on-line as the current that elicited one spike (stimulus TestAmp, duration 1350 ms). The stimuli used in the experiments, for fitting and testing the models were: • IDRest: current step of 1350 ms, injected at different amplitude levels in 25% increments (range 50-300% threshold). IDRest was renamed to Step in the model.
• SponNoHold: the first 10 seconds of this stimulus was used to calculate the resting membrane potential. No holding or stimulation currents were applied.
• SponHold: the first 10 seconds of this stimulus was used to calculate the holding current applied to keep the cells at the target potential.
• PosCheops: ramps of current from 0 to 300% and from 300 to 0% having progressively shorter durations (4000 ms, 2000 ms, 1250 ms). To test the models in tonic mode we used the first increasing ramp in the stimulus, while we used the last one in the bursting firing mode. We chose the last one because the biological cells were more likely to generate a burst, while they were silent during the first two ramps.
• NOISEOU3: the original wave was scaled and offset for each cell based on the spike frequency in response to the IDRest protocol. The scaling factor w was extracted from the frequencycurrent curve and corresponded to the current value that made the cell fire at 7.5 Hz.
Neurons that were completely stained and those with high contrast were reconstructed in 3D and corrected for shrinkage as previously described [16]. Reconstruction used the Neurolucida system (MicroBrightField). The location of the stained cells was defined by overlaying the stained slice and applying manually an affine transformation to the Paxinos and Watson's rat atlas [19].

Electrical features extraction
Electrical features were extracted using the Electrophys Feature Extraction Library (eFEL) [42]. We calculated the adaptation index (AI) from recordings in tonic mode (Step 200% threshold) and classified TC VB neurons into adapting (AI> = 0.029) and non-adapting (AI<0.029) electrical types. AI was calculated using the eFEL feature adaptation_index2 and corresponded to the average of the difference between two consecutive inter-spike intervals (ISI) normalized by their sum. The cut-off value was calculated after fitting a Gaussian mixture model to the bimodal data, using available routines for R [43,44]. In order to group data from different cells and generate population features, we normalized all the stimuli by the rheobase current I thr of each cell. To calculate I thr , we used IDRest and IDThresh and selected the minimal amplitude that evoked a single spike. Along with the voltage features, we extracted mean holding and threshold current values for all the experimental stimuli. Description of the features and the details on their calculation are available on-line [42]. Current stimuli applied during the optimization and generalization were directly obtained from the experimental values or automatically calculated by following the experimental procedures (e.g. noise stimulus).

Morphology analysis
Reconstructed morphologies were analyzed to objectively identify different morphological types. The Sholl profiles of each pair of cells was statistically tested by using k-samples Anderson-Darling statistics. This test was preferred to the most common Kolmogorov-Smirnov test, because it does not assume that the samples are drawn from a continuous distribution. The different Sholl profiles are indeed an analysis of the intersections with discrete spheres.
To compare the topological description of each morphology we transformed the persistence barcodes into persistence images and calculated their distances as in [20]. Briefly, we converted the persistence barcode, which encodes the start and end radial distances of a branch in the neuronal tree, into a persistence diagram. In the persistence diagram, each bar of the barcode is converted into a point in a 2D space, where the X and Y coordinates are the start and end radial distances of each bar. The persistence diagram was then converted into a persistence image by applying a Gaussian kernel. We used the library NeuroM [45] to perform Sholl and morphometrics analyses.

Ionic currents models
We used Hodgkin-Huxley types of ionic current models, starting from kinetics equations already available in the neuroscientific literature. Along with kinetics of the ionic currents, we stored information on the experimental conditions, such as temperature and LJP, by using the software NeuroCurator [46]. Whenever the data was available, we compared simulated voltage-clamp experiments to experimental data from juvenile rats. Ionic currents I i were defined as functions of the membrane potential v, its maximal conductance density g i and the constant value of the reversal potential E i : i h y i ðv À E i Þ m ion and h ion represent activation and inactivation probability (varying between 0 and 1), with integer exponents x and y. Each probability varied according to: n 0 ðvÞ ¼ ðn 1 ðvÞ À nÞ=t n ðvÞ where n 1 (v) is a function of voltage that represents the steady-state activation/inactivation function (normally fitted with a Boltzmann curve) and τ n (v) is a voltage-dependent time constant. Exceptions to this formalism are ionic currents that do not inactivate (y = 0) and ionic currents with (in)activation processes mediated by two or more time constants. Calcium currents (I CaT and I CaL ) were modeled according to the Goldman-Hodgkin-Katz constant field equation and had permeability values instead of conductance [47]. Fast transient sodium current I NaT and delayed potassium current I Kd . I NaT and I Kd were taken from a previous models of rat TC neurons from the VB nucleus [12], available on SenseLab ModelDB (accession no. 279). I NaT was compared with recordings of transient sodium currents in P7-11 rat neurons from the dorsolateral geniculate (dLGN) nucleus [48].
Low-threshold activated (T-type) calcium current I CaT . I CaT model was taken from [12] and available on-line (ModelDB, accession no. 279). This model was based on data recorded from VB neurons of Sprague-Dawley rats (P7-12) at room temperature and corrected for −9 mV LJP [11].
Hyperpolarization-activated cationic current I h . The steady-state activation for I h was derived from VB thalamic neurons in P10-20 Long-Evans rats and was already corrected for −10 mV LJP in the original publication [49], The equation used was: The time constant of activation was modeled as in [50], which derived a mathematical description of I h based on data from the dLGN in adult guinea pigs [51]. The equation describing the time dependence of activation was: t m ¼ 1=½expðÀ 14:59 À 0:086vÞ þ expðÀ 1:87 þ 0:0701vÞ� The equilibrium potential of the channel E H was −43 mV. Persistent sodium current I NaP . We modeled I NaP as in [17] which based their model on recordings from entorhinal neurons of Long-Evans rats (P25-P35) [52]. The model is available in ModelDB, accession no. 139653. The steady-state activation was modified according to [48] and the steady-state inactivation according to [14]. The original steady-state activation data were recorded at room temperature (22-24˚C) and corrected for −6/−7 mV LJP.
Fast transient (A-type) potassium current I KA . The mathematical formulation of I KA was based on data recorded from VB neurons in Sprague-Dawley rats (P7-15), recorded at room temperature (22-24˚C) [11]. A Q 10 = 2.8 was experimentally determined and used for simulations at different temperatures. In the original experiments a small LJP (<−4 mV) was measured and not corrected. The current had a rapid and a slow component, represented by two activation and two inactivation variables. The model of this current was provided by the authors of [14].
High-threshold (L-type) calcium current I CaL . I CaL model is the same as TC neurons model previously published [14]. The model was based on data from isolated guinea-pig hippocampal neurons, recorded at room temperature (20-22˚C) with modifications to the Boltzmann curve parameters of activation contained in the correction to the original models. A small LJP (<3 mV) was not corrected [50]. A Q 10 of 3 was used for simulations at different temperatures.
Calcium-activated potassium currents. TC neuron express genes for BK-type [53] and SK-type calcium-activated potassium channels [54]. Models of BK-type currents, similar to the I C current, have already been used to model TC neurons [14,50,53]. However, data characterizing this current in mammalian neurons are not available. We thus included only a model of I SK (available on ModelDB, accession no. 139653) based on rat mRNA expression data in Xenopus oocytes [55]. Intracellular calcium dynamics. A simple exponential decay mechanism was used to model the intracellular calcium dynamics (ModelDB, accession no. 139653). Both I CaT and I CaL contributed to the intracellular calcium concentration.
In addition, we included a voltage-insensitive membrane current I leak . The equilibrium potential was −79 mV and corresponded to the average resting potential from our experimental recordings. NEURON 7.5 software was used for simulation [56]. We used NEURON variable time step method for all simulations. For the sake of spatial discretization, each section was divided into segments of 40 μm length. The following global parameters were set during optimization and generalisation: initial simulation voltage (−79 mV), simulation temperature (34˚C), specific membrane capacitance (1 μF/cm 2 ), specific intracellular resistivity 100 Ocm for all the sections, equilibrium potentials for sodium and potassium were 50 mV and −90 mV, respectively.

Simulation and parameters optimization
BluePyOpt [18] with Indicator Based Evolutionary Algorithm (IBEA) were used to fit the models to the experimental data. Each optimization run was repeated with three different random seeds and evaluated 100 individuals for 100 generations. The evaluation of these 300 individuals for 100 generations was parallelized using the iPython ipyparallel package and took between 21 and 52 h on 48 CPU cores (Intel Xeon 2.60 GHz) on a computing cluster. Each optimization run typically resulted in tens or hundreds of unique acceptable solutions, defined as models having all feature errors below 3 STD from the experimental mean.

Sensitivity analysis
We performed a sensitivity analysis of an optimization solution by varying one parameter value (p m ) at a time and calculating the electrical features from the voltage traces (y + and y -). We defined the sensitivity as the ratio between the normalized feature change and the parameter change, which for smooth functions approximates a partial derivative [57,58]. The features changes were normalized by the optimized feature value. For small changes of parameter values, we assumed that the features depend linearly on its parameters. We could thus linearize the relationship between the features and the parameters around an optimized parameter set and calculate the derivatives. The derivatives were calculated with a central difference scheme [57]. @y n @p m � y þ n À y À n 2Dp m We collected the derivatives (sensitivities) in the N X M Jacobian matrix, with N representing the number of features and M the number of parameters.
To rank parameters and features we computed their relative importance by calculating their norms (the square root of the summed squared values) from the Jacobian columns and rows, respectively. To cluster parameters based on similar influences on the features and to cluster features that were similarly dependent on the parameters, we used angles between columns (or rows) to compute distances D between parameters (or features): Features where thus considered similar if they depended in a similar manner on the parameters, independent of sign or magnitude.
Supporting information S1 Fig. Morphological analysis. (A) Morphometrics of the thalamocortical (TC) morphological type. Each histogram shows basic morphometrics at the level of the neuron (first row) or at the level of dendritic trees (second row). (B) Distance matrix between persistence barcodes for all TC morphologies. Related to Fig 2A). (C) P-values for the k-samples Anderson-Darling statistics. It tests the null hypothesis that the Sholl profiles of each pair of morphologies are drawn from the same population. The p-values are not corrected for multiple comparisons and show that we cannot reject the null hypothesis for most of the morphology pairs (at 0.05 significance level). Related to Fig 2B. (TIF)