Degeneracy in the robust expression of spectral selectivity, subthreshold oscillations and intrinsic excitability of entorhinal stellate cells

Biological heterogeneities are ubiquitous and play critical roles in the emergence of physiology at multiple scales. Although neurons in layer II (LII) of the medial entorhinal cortex (MEC) express heterogeneities in their channel properties, the impact of such heterogeneities on the robustness of cellular-scale physiology has not been assessed. Here, we performed a 55-parameter stochastic search spanning 9 voltage- or calcium-activated channels to assess the impact of channel heterogeneities on the concomitant emergence of 10 electrophysiological characteristics of LII stellate cells (SCs). We generated 50,000 models and found a heterogeneous subpopulation of 155 valid models to robustly match all electrophysiological signatures. We employed this heterogeneous population to demonstrate the emergence of cellular-scale degeneracy in LII SCs, whereby disparate parametric combinations expressing weak pairwise correlations resulted in similar models. We then assessed the impact of virtually knocking out each channel from all valid models and demonstrate that the mapping between channels and measurements was many-to-many, a critical requirement for the expression of degeneracy. Finally, we quantitatively predict that the spike-triggered average of LII SCs should be endowed with theta-frequency spectral selectivity and coincidence detection capabilities in the fast gamma-band. We postulate this fast gamma-band coincidence detection as an instance of cellular-scale efficient coding, whereby SC response characteristics match the dominant oscillatory signals in LII MEC. The heterogeneous population of valid SC models built here unveils the robust emergence of cellular-scale physiology despite significant channel heterogeneities, and forms an efficacious substrate for evaluating the impact of biological heterogeneities on entorhinal network function. KEY POINTS Stellate cells (SC) in layer II (LII) of the medial entorhinal cortex express cellular-scale degeneracy in the concomitant manifestation of several of their unique physiological signatures. Several disparate parametric combinations expressing weak pairwise correlations resulted in models with very similar physiological characteristics, including robust theta-frequency membrane potential oscillations spanning several levels of subthreshold depolarization. Electrophysiological measurements of LII SCs exhibited differential and variable dependencies on underlying channels, and the mapping between channels and measurements was many-to-many. Quantitative predictions point to theta-frequency spectral selectivity and fast gamma-range coincidence detection capabilities in class II/III spike-triggered average of LII SCs, with the postulate for this to be an instance of cellular-scale efficient coding. A heterogeneous cell population that accounts for both channel and intrinsic heterogeneities in LII SCs, which could be employed by network models of entorhinal function to probe the impact of several biological heterogeneities on spatial navigation circuits.

Neurons in layer II (LII) of the rodent medial entorhinal cortex (MEC) have been implicated in spatial navigation, especially since the cells in LII MEC are known to act as grid cells that elicit action potentials in a grid-like pattern as the animal traverses an arena (Hafting et al., 2005;Moser et al., 2008;Buzsaki & Moser, 2013;Moser et al., 2014;Ray et al., 2014;Tang et al., 2014;Moser et al., 2015). Ever since the discovery of grid cells, several theoretical and computational models have been proposed for their emergence, and have been tested from several different perspectives with varying degrees of success (O'Keefe & Burgess, 2005;Burak & Fiete, 2006;Burgess et al., 2007;Jeewajee et al., 2008;Welinder et al., 2008;Burak & Fiete, 2009;Giocomo et al., 2011b;Sreenivasan & Fiete, 2011;Navratilova et al., 2012;Couey et al., 2013;Domnisoru et al., 2013;Schmidt-Hieber & Hausser, 2013;Yoon et al., 2013;Bush & Burgess, 2014;Rowland et al., 2016;Schmidt-Hieber et al., 2017). Although these models and associated experiments have provided significant insights into entorhinal function, a lacuna common to these models relates to the systematic assessment of the impact of the different biological heterogeneities in the medial entorhinal cortex. Specifically, a systematic evaluation of the role of different forms of network heterogeneities, including those in channels, structural, intrinsic and synaptic properties and in afferent connectivity, with reference to entorhinal physiology has been lacking.
A first and essential step in addressing these and other related questions on the impact of biological heterogeneities on entorhinal network function is to assess the robustness of cellular physiology in the presence of well-established heterogeneities in channel properties. Specifically, measurements of channel properties, including kinetics, voltage-dependent gating and conductance values, from entorhinal neurons are known to exhibit significant variability across neurons (Bruehl & Wadman, 1999;Magistretti & Alonso, 1999;Dickson et al., 2000;Fransen et al., 2004;Castelli & Magistretti, 2006;Dudman & Nolan, 2009;Pastoll et al., 2012;Schmidt-Hieber & Hausser, 2013). Despite this, entorhinal neurons exhibit signature electrophysiological characteristics that robustly fall into specific ranges for each physiologically relevant measurement. How do these neurons achieve such cellular-scale robustness in their physiology despite widely variable conductances and channel properties? Is there a requirement for individual channels or pairs of channels to be maintained at specific levels with specific properties for signature electrophysiological properties to emerge? Is there a one-to-one mapping between individual channels and the physiological properties that they regulate?
In this study, we build a conductance-based intrinsically heterogeneous population of LII stellate cell (SC) models of the medial entorhinal cortex that satisfied several of their unique electrophysiological signatures. We employed this heterogeneous population of LII SC models to demonstrate the expression of cellular-scale degeneracy (Edelman & Gally, 2001) in the concomitant emergence of these measurements.
Specifically, we showed that LII SC with very similar electrophysiological characteristics emerged from disparate channel and parametric combinations expressing weak pairwise correlations. We employed these models to demonstrate the differential and variable dependencies of measurements on underlying channels, and showed that the mapping between channels and measurements was many-to-many. Finally, we employed this electrophysiologically validated model population to make quantitative testable predictions that the spike triggered average (STA) of LII SCs should be endowed with theta-frequency spectral selectivity and coincidence detection capabilities in the fast gamma-band. We postulate this fast gamma-band coincidence detection to be an instance of cellular-scale efficient coding (Narayanan & Johnston, 2012) that matches neuronal response properties to the dominant oscillatory band in the superficial layers of MEC (Colgin et al., 2009;Colgin & Moser, 2010;Colgin, 2016;Trimper et al., 2017). The heterogeneous population of models built here also forms an efficacious substrate for constructing network models of the entorhinal cortex, towards assessing the impact of cellular and channel properties and associated heterogeneities on emergent behavior such as grid cell formation.

METHODS
We employed a single compartmental cylinder model of 70-µm diameter and 75-µm length (Fig. 1A). The choice of a single compartmental model was largely driven by the absence of direct and detailed electrophysiological characterization of dendritic intrinsic properties or of ion channels that express in LII SCs. As a consequence, morphologically precise models with specific channel expression profiles and matched intrinsic properties were infeasible. On the other hand, as the somatic channel properties and intrinsic physiological measurements of LII SCs are well characterized, we employed a singlecompartmental model that wouldn't have to make explicit or implicit assumptions about dendritic physiology. Additionally, as a goal of this study was to develop an intrinsic heterogeneous model population of LII SCs towards their incorporation into network models, it was essential to ensure that the computational complexity of single neurons was minimal. A single compartmental conductance-based model population that was endowed with the different ion channels and systematically reflects intrinsic heterogeneities in LII SCs served as an efficacious means to achieve this goal as well.

Passive and active neuronal properties
Passive properties were incorporated into the model as RC circuit that was defined through a specific membrane resistance, R m and a specific membrane capacitance, C m .
All channel models were based on Hodgkin-Huxley formulation (Hodgkin & Huxley, 1952) except for the SK channel, which was modeled using a six-state Markovian kinetics model. Sodium, potassium and HCN channel currents were modeled using the Ohmic formulation and calcium channels followed the Goldman-Hodgkin-Katz (GHK) formulation (Goldman, 1943;Hodgkin & Katz, 1949) for current computation.
The reversal potentials for Na + , K + and HCN channel were 50, -90 and -20 mV, respectively. Calcium current through voltage gated calcium channels contributed to cytosolic calcium concentration ([Ca] c ), and its decay was defined through simple first order kinetics (Destexhe et al., 1993;Poirazi et al., 2003;Carnevale & Hines, 2006;Narayanan & Johnston, 2010;Honnuraiah & Narayanan, 2013): Channel models were directly adopted in cases where they were explicitly based on direct measurements from LII SC channels. In cases where such explicit models were not available, model formulations were taken from other cell types and were explicitly adapted to match direct electrophysiological measurements. As channel models were either adopted from different studies or were adapted to match experimental observation, in what follows, we provide details of the models that we employed for gating and kinetics of each channel. The parameters that define these channel models are described in Table 1, along with the base values of the 55 passive and active parameters that govern these models. In channels that employed the Hodgkin-Huxley formulation, the model evolved by modifications to one or two gating particles, with each gating particle following first-order kinetics as follows: where m ∞ and τ m respectively defined the steady-state value and the time constant of the state variable that governed the gating particle. Channel gating and kinetics were appropriately adjusted for temperature dependence from corresponding experimental measurements.

The fast sodium channel
The NaF model was adopted from (Dudman & Nolan, 2009), and the current through this sodium channel was: The activation gating particle was defined by: The inactivation gating particle was defined by:

The delayed rectifier potassium channel
The KDR model was adopted from (Dudman & Nolan, 2009), and the current through this potassium channel was: The activation gating particle was governed by the following equations:

The hyperpolarization-activated cyclic-nucleotide-gated channel
The HCN channel model was adopted from (Dickson et al., 2000;Fransen et al., 2004;Schmidt-Hieber & Hausser, 2013) and the current through this nonspecific cationic channel was: where ms HCN and mf HCN respectively defined the gating variables for the slow and fast components of the current through HCN channels, and HCN ms mf defined the ratio of the fast to slow HCN conductance values. The activation gating particles for the slow and fast HCN components were governed by the following equations: The persistent sodium channel The NaP model was adopted from (Magistretti & Alonso, 1999;Dickson et al., 2000;Fransen et al., 2004), and the current through this sodium channel was: The activation gating particle was defined by: The inactivation gating particle was defined by:

The transient A-type potassium channel
The KA model was adopted from (Dudman & Nolan, 2009), and the current through this potassium channel was: The activation gating particle was governed by the following equations: The inactivation gating particle was governed by the following equations:

The high voltage-activated calcium channel
The HVA calcium channel model was fitted with corresponding electrophysiological data (Bruehl & Wadman, 1999;Castelli & Magistretti, 2006). The current through this channel followed GHK conventions, with the default extracellular and cytosolic calcium concentrations set at 2 mM and 100 nM, respectively. The conductance evolution of this channel was modeled as follows: The activation and inactivation gating particles were governed by the following equations:

The low voltage-activated calcium channel
The LVA calcium channel model was fitted with corresponding electrophysiological data from (Bruehl & Wadman, 1999;Pastoll et al., 2012). The current through this channel followed GHK conventions, with the default extracellular and cytosolic calcium concentrations set at 2 mM and 100 nM, respectively. The conductance evolution of this channel was modeled as follows: The M-type potassium channel The KM model was adopted from (Shah et al., 2008), and the current through this potassium channel was: The activation gating particle was governed by the following equations:

Intrinsic measurements
We measured the resting membrane potential of the model neuron by allowing the model to settle at a steady-state potential when no current was injected for a period of 5 s. This was essential because there were several slow subthreshold conductances (especially the calcium-activated potassium channel) that contributed to resting membrane potential in our model. We set the passive membrane potential, in the absence of any subthreshold conductance, to be at -77 mV (Boehlen et al., 2013) and allowed the interactions among the several subthreshold conductances to set the steady-state resting membrane potential (V RMP ). After this initial 5 s period of the simulation, V RMP was computed as the mean of the membrane potential over a one-second interval (5 th to 6 th second; Fig. 1B). We also calculated the standard deviation (SD) of membrane potential over the same 1 s period to ensure that the RMP was measured after attainment of steady state. All intrinsic measurements reported below were always performed after this 6 s period (5 s for the transients to settle to steady-state and 1 s for RMP measurement).
To estimate sag ratio in the model, we injected a hyperpolarizing step current of 200 pA amplitude for 1s and recorded the voltage response. Sag ratio (Sag) was computed as the membrane potential deflection achieved at steady state (V SS ) during the current injection period divided by the peak deflection of the membrane potential (V peak ) within the period of current injection (Fig. 1C). In assessing supra-threshold excitability of the model neuron, we measured the number of action potentials (AP) elicited by the neuron in response to different depolarizing step current injections spanning 500 ms. We defined the number of APs fired for 100 and 400 pA current injections as N 100 (Fig. 1D) and N 400 (Fig. 1E), respectively.
Input resistance (R in ) was calculated from the steady-state voltage response (after 1 s of current injection) of the model neuron to subthreshold current pulses of amplitudes spanning -100 pA to 100 pA in steps of 20 pA. The steady state voltage response was plotted against the corresponding amplitude of injected current, and the slope of a linear fit to this plot was assigned as the input resistance of the model (Fig. 1F). Spike amplitude (V AP ) was computed from the first AP elicited during a 400-pA step current injection, and was defined as the difference between V RMP and the peak membrane potential achieved during the AP (Fig. 1G).
As entorhinal stellates reside within an oscillatory network, it was essential that excitability measures be computed in a frequency-dependent manner. To do this, we computed well-established impedance-based measurements from its amplitude and phase profiles (Hutcheon et al., 1996a, b;Hutcheon & Yarom, 2000;Haas & White, 2002;Erchova et al., 2004;Giocomo et al., 2007;Narayanan & Johnston, 2007. € where Re(Z(f)) and Im(Z(f)) were the real and imaginary parts of the impedance Z(f), respectively. The frequency at which |Z(f)| reached its maximum value was considered as the resonance frequency, f R , and resonance strength (Q R ) was defined as the ratio of |Z(f R )| to |Z(0.5)| (Fig. 1H). The impedance phase profile ϕ(f) was computed as: The total inductive area, Φ L , defined as the area under the inductive part of ϕ(f), was calculated based on the impedance phase profile : ∫ Sub-and peri-threshold membrane potential oscillations (MPO) were assessed in voltage responses of the model to depolarizing pulse current injections spanning 100-300 pA in steps of 10 pA, each lasting for 5 s (Fig. 1I). The last 3 s period of this 5-s period was transformed to frequency domain through the Fourier transform, and the frequency at which this spectral signal had maximum magnitude was defined as the MPO frequency (f MPO ). We defined the peri-threshold oscillation frequency (f osc ) as the f MPO of the subthreshold voltage response proximal to the spiking threshold of the model.

Multi-parametric Multi-objective Stochastic Search Algorithm
To generate an intrinsically heterogeneous population of LII SCs and to assess if the concomitant functional expression of all 10 intrinsic measurements manifested degeneracy in terms of the specific ion channel combinations that can elicit them, we employed a multi-parametric multi-objective stochastic search (MPMOSS) algorithm (Foster et al., 1993;Goldman et al., 2001;Prinz et al., 2003;Marder & Taylor, 2011;Rathour & Narayanan, 2012, 2014Anirudhan & Narayanan, 2015;Srikanth & Narayanan, 2015;Mukunda & Narayanan, 2017). This stochastic search was performed over 55 parameters (Table 1) and jointly validated against 10 sub-and supra-threshold measurements ( Fig. 1; V RMP , SD, Sag ratio, R in , f R , Q R , f osc , N 100 , N 400 , V AP ) towards matching electrophysiological recordings from LII SCs (Table 2). In executing the MPMOSS algorithm, we constructed a model neuron from specific values for each of the 55 parameters, each of which was independently and randomly picked from a uniform distribution whose bounds reflected the electrophysiological variability in that parameter (Table 1). For each such randomly chosen model, which ensured that we are not biasing our parametric ranges with any constraints, all 10 intrinsic measurements were computed and were compared against their respective electrophysiological bounds (Table 2). A model that satisfied all the 10 criteria for validation was declared valid. We repeated this procedure for 50,000 randomized picks of the 55 parameters, and validated these models against the 10 measurements. Intrinsic heterogeneity and degeneracy were then assessed in the resulting population of valid LII SC models and their parametric combinations.
This assessment was performed through the analysis of parametric distributions and pairwise correlations among valid model parameters.
To assess the impact of individual channels on each of the 10 intrinsic measurements within the valid model population, we employed the virtual knockout model (VKM) approach (Rathour & Narayanan, 2014;Anirudhan & Narayanan, 2015;Mukunda & Narayanan, 2017). In doing this, we first set the conductance value of each of the 9 active ion channels independently to zero for each of the valid models. We then computed all the 10 intrinsic measurements for each model, and assessed the sensitivity of each measurement to the different channels from the statistics of post-knockout change in the measurements across all valid models. When some of the channels were knocked out, certain valid models elicited spontaneous spiking or showed depolarization-induced block (when depolarizing current were injected). These VKMs were not included into the analysis for assessing the sensitivities, because this precluded computation of all 10 measurements from such models.

Spike triggered average and associated measurements
For estimation of STA, a zero-mean Gaussian white noise (GWN) with a standard deviation σ noise was injected into the neuron for 1000 s. σ noise was adjusted such that overall action potential firing rate was ~1 Hz in the model under consideration. This ensured that the spikes were isolated and aperiodic, thereby establishing statistical independence of the current samples used in arriving at the STA (Aguera y Arcas & Fairhall, 2003;Das & Narayanan, 2014. The STA was computed from the injected current for a period of 300 ms preceding the spike and averaged over all spikes across the time period of simulation, translating to around 1000 spikes for each STA computation. STA kernels were smoothed using a median filter spanning a 1 ms window for representation purposes and for computing quantitative measurements that were derived from the STA.
Quantitative metrics for spectral selectivity in the STA, for coincidence detection windows and intrinsic excitability were derived from the STA and its Fourier transform (Das & Narayanan, 2014. Specifically, the frequency at which the magnitude of the Fourier transform of the STA peaked was defined as the STA characteristic frequency (f STA ). STA selectivity strength (Q STA ) was defined as the ratio of |STA(f STA )| to |STA(0.5)|. The peak positive current in the STA kernel was defined as I STA peak , which constitutes a measure of excitability. To quantify the window for integration/coincidence detection, we defined the spike-proximal positive lobe (SPPL) as the temporal domain that was adjacent to the spike where the STA was positive (Das & Narayanan, 2015. The total coincidence detection window, CDW (T TCDW ) was computed as the temporal distance from the spike location (t=0 ms) to the first zero crossing in the STA. T TCDW constitutes the entire temporal expanse over which the inputs were positively weighted and hence covered the entire temporal spread of SPPL. To account for the specific shape of the STA in defining the coincidence detection window, we defined an effective CDW (T ECDW ), which was a STA-weighted measure of SPPL (Das & Narayanan, 2015: All simulations were performed using the NEURON programming environment (Carnevale & Hines, 2006) at 34° C, with a simulation step size of 25 µs. All data analyses and plotting were executed using custom-written software within the IGOR pro environment (Wavemetrics Inc.). All statistical analyses were performed using the R statistical package (R Core Team, 2013).

RESULTS
The principal goal of this study was to assess the impact of molecular-scale heterogeneities in channel properties on cellular-scale physiological signatures of LII stellate cells of the medial entorhinal cortex. We approached this by building an unbiased stochastic search-based conductance-based population of LII SCs that satisfied several of their unique electrophysiological signatures. Apart from providing an efficacious substrate for understanding the roles of channel parameters, intrinsic measurements and associated heterogeneities in entorhinal function, our goal in building these models was three fold. First, a heterogeneous population of LII SC models that satisfied several electrophysiological constraints would provide us the means to assess if there was significant degeneracy in the emergence of these measurements, or if there was a requirement on unique mappings between channel properties and physiological measurements. Second, such a heterogeneous population would allow us to establish the specific roles of different channels in mediating or regulating different physiological properties, and assess variability in such regulatory roles. Third, and importantly, as these models match with LII SC models in several ways, they provide an ideal foundation for making quantitative predictions about stellate cell physiology, which can be electrophysiologically tested.
Towards this, we first hand-tuned a conductance-based, biophysically and physiologically relevant base model that was characterized as a single compartmental cylindrical model with different passive and active properties (see Methods). The model was endowed with 9 different active ion channels (besides leak channels), and matched with 10 distinct electrophysiological measurements obtained from LII SCs (Fig. 1, Table   2). These electrophysiological measurements included the significant sag observed in response to pulse current injections (Fig. 1C), theta-frequency membrane potential resonance that exhibited strong spectral selectivity (Fig. 1H), and importantly the robust subthreshold membrane potential oscillations in the theta-frequency range at different depolarized potentials (Fig. 1I). The 55 active and passive parameters that governed LII SC models, and their respective values in the base model are listed in Table 1.

Diverse depolarization-dependent evolution of membrane potential oscillations in stochastically generated stellate cell models
We employed the base model parameters (Table 1) as the substrate for a multi-parametric stochastic search algorithm for models that would meet multiple objectives in terms of matching with the electrophysiological properties of LII SCs. The range of individual parameters over which this multi-parametric (55 parameters) multi-objective (10 measurements) stochastic search (MPMOSS) algorithm was executed is provided in Table 1. We generated a test population of 50,000 model cells by sampling these model parameters, and subjected these model cells to validation based on the physiologically observed range of sub-and supra-threshold measurements (Table 2). First, we found a subpopulation of these models where all measurements, except for the ability to express theta-frequency membrane potential oscillations, were within the specified bounds. We depolarized this subpopulation of models and asked if these models expressed robust subthreshold oscillations in their membrane potentials.
We found that the depolarization-dependent evolution of sub-and supra-threshold (regular spiking behavior) oscillations exhibited significant diversity across different models within this subpopulation (Fig. 2). In most models within this subpopulation, subthreshold oscillations and regular spiking with increasing current injections (Fig. 2G).

The stochastic search strategy yielded an intrinsically heterogeneous population of models that matched several electrophysiological signatures of stellate cells
Out of 50,000 models that were generated as part of the stochastic search strategy, only electrophysiological measurements (Table 2), including their ability to manifest robust peri-threshold theta frequency oscillations. We plotted all the 10 electrophysiological measurements in this model population to assess if they were clustered around their respective base model values (Fig. 1) or if they were distributed to span the range of valid model measurements (Table. 2). Whereas a clustered set of measurements would have implied a near-homogeneous population of models, a distributed pattern that spanned the range of respective electrophysiological counterparts would provide us with a heterogeneous model population that reflects experimental variability in the respective measurements. We found all measurements to spread over a large span, with most of them covering the entire min-max range of their respective bounds ( Fig. 3A; note that N 100 has not been plotted because it was required to be zero for model validity). We observed the emergence of sub-and supra-threshold membrane potential oscillations in these models when the average membrane voltage was between -59 mV to -45 mV (Fig.   3B). To distinguish between sub-and supra-threshold oscillations, we plotted the frequency of these membrane potential oscillations against their peak-to-peak amplitude ( Fig. 3C). Two clearly separable clusters were observed, with all subthreshold oscillations clustered at the low-amplitude range (< 25 mV), while the action potentials forming a cluster with amplitudes greater than ~60 mV. Importantly, these observations also demonstrate that the characteristics of membrane potential oscillations exhibited significant heterogeneity, thus matching the electrophysiological variability observed in LII SCs. Together, the biophysically and physiologically constrained MPMOSS algorithm yielded a population of LII SC models that manifested considerable intrinsic heterogeneities (Fig. 3) matching the ranges observed in corresponding electrophysiological measurements (Table 2).

The valid model population manifested cellular-scale degeneracy
What were the specific constraints on the 55 different parameters in yielding the valid model population that concomitantly matched several electrophysiological signatures of stellate cells? Were these parameters clustered around specific values thereby placing significant constraints on the different channels, their properties and their expression profiles? Would individual channels have to be maintained at specific expression levels for models to match all 10 electrophysiological signatures? To address these questions, we first randomly picked 5 of the 155 valid models that exhibited similar measurements ( Fig. 4A-H), and asked if the set of parameters governing these models were also similar.
We normalized each of the 55 parameters from these 5 models with reference to their respective min-max values, and found none of these models to follow any specific trend in their parametric values with most parameters spreading across the entire range that they were allowed to span (Fig. 4I). These observations provided the first line of evidence for the expression of degeneracy in stellate cell models, whereby models with very similar functional characteristics emerged from disparate parametric combinations. To confirm this across all valid models, we plotted the histograms of each parameter for all the 155 valid models, and found their spread to span the entire testing range for all parameters (Fig. 5A, bottom row).
Although the distributions of individual parameters span their respective ranges, was it essential that there are strong pair-wise constraints on different parameters towards obtaining valid models? To address this question and to test if parametric combinations resulting in the heterogeneous valid population were independent of each other, we plotted the pair-wise scatter plots across all 55 parameters for all the 155 valid models (Fig. 5A). We computed the Pearson's correlation coefficient for each pair-wise scatter plot (Fig. 5B), and found very weak correlations across different parameters ( Figure 5B- Together, these results demonstrated that the ability of a heterogeneous model population to concomitantly match several electrophysiological signatures of LII SCs was attainable even when individual ion channels were expressed with disparate densities and properties and when channels did not express strong pairwise correlations. These provided strong lines of evidence for the expression of cellular-scale degeneracy in LII SCs, whereby disparate combinations of channels with distinct parameters were robustly capable of eliciting similar functional characteristics.

Virtual knockout models: A many-to-many mapping between individual channels and physiological characteristics enabled the expression of degeneracy
A crucial requirement for the expression of such cellular-scale degeneracy is the ability of several channels to regulate different physiological characteristics (Drion et al., 2015;Rathour et al., 2016). In the absence of such capabilities, the system in essence will comprise of several one-to-one mappings between channels and physiological characteristics, thereby requiring the maintenance of individual channels at specific expression levels. In asking if there was a many-to-many mapping between channels and physiological properties, we employed virtual knockouts of individual channels on all models within the heterogeneous valid model population to assess the impact of their acute removal on physiology. We employed these virtual knockout models (VKMs) to assess the impact of individual channels on all the 10 physiological measurements by calculating the change observed in each measurement after setting individual conductance values to zero (Rathour & Narayanan, 2014;Anirudhan & Narayanan, 2015;Mukunda & Narayanan, 2017).
As expected, V RMP (Fig. 6A) was largely unaffected by the knockout of suprathreshold conductances (KDR, NaF and HVA) with all subthreshold channels showing differential and heterogeneous regulation of V RMP . Specifically, consistent with prior electrophysiological recordings (Dickson et al., 2000), HCN VKMs showed large and variable hyperpolarizing shifts in V RMP with reference to their respective base valid models. In addition, knockout of KM or SK channels resulted in variable depolarizing shifts to V RMP , but VKMs of NaP, KA and LVA channels did not significantly alter V RMP .
Although sag ratio was expectedly (Dickson et al., 2000) reliant on HCN channels (Fig.   6B), there were other channels, including NaP, LVA, SK and KM channels, which also significantly contributed to the specific value of sag ratio. Input resistance (Fig. 6C) was critically altered by HCN channel knockouts (Dickson et al., 2000), with other subthreshold channels including KM and SK also showing large impacts on R in .
Resonance strength (Fig. 6D) and resonance frequency (Fig. 6E) were dramatically reduced in HCN knockouts (Erchova et al., 2004;Boehlen et al., 2013), with NaP, KM and SK channel knockouts also showing significant changes in both measurements. Of all the assessed measurements, we found the frequency of perithreshold membrane potential oscillations to be the most sensitive measurement to channel knockouts, with most channels having a significant, yet variable, effect on f osc (Fig. 6F). We observed the most reliable and least variable effect on f osc with the removal of NaP, which consistently resulted in the loss of peri-threshold oscillations across all models (Fig. 6F). This is consistent with several studies that have demonstrated the importance of persistent sodium channels in the emergence of peri-threshold MPOs in LII SCs (Alonso & Llinas, 1989;Klink & Alonso, 1993;Dickson et al., 2000;Boehlen et al., 2013). Although NaP was the dominant channel that acted as the amplifying conductance (Hutcheon & Yarom, 2000) in the emergence of peri-threshold MPOs, we found heterogeneity across models in the specific resonating conductance that enabled these MPOs. Specifically, whereas in some models, the removal of HCN channels resulted in complete loss of MPOs, in other models the same result was achieved by the removal of KM channels. These observations suggested that the two conductances, HCN and KM, synergistically contributed as resonating conductances towards the emergence of perithreshold MPOs in stellate cells (Nolan et al., 2007;Boehlen et al., 2013). Although most of the other channels showed regulatory capabilities in terms of regulating peri-threshold f osc , unlike NaP, HCN and KM channels, their removal did not result in complete elimination of MPOs in most models. The critical dependence of f osc on KDR removal was just a reflection of high excitability of KDR VKMs, where the cells either spontaneously fired or abruptly switched to regular spiking with small current injections resulting in the complete absence of peri-threshold oscillations.
N 100 was significantly higher with the deletion of KM or SK channels, with little or no effect with deletion of other channels (Fig. 6G). As spiking activity is critically reliant on KDR and NaF channels, their removal significantly reduced N 400 (Fig. 6H). In addition, whereas the removal of the subthreshold regenerative conductance NaP reduced N 400 , knocking out the subthreshold restorative conductances, KM, SK and KA resulted in enhanced action potential firing that increased N 400 to variable degrees (Fig. 6H).
Although the two calcium channels mediate inward currents, their removal resulted in an increase (rather than a decrease) in N 400 because of the presence of SK channels.
Specifically, when either the HVA or the LVA channels was removed, the inward calcium current and cytosolic calcium concentration were lower, thereby resulting in lesser activation of SK channels and consequently leading to higher excitability ( Fig. 6G-H). Finally, AP amplitude ( Fig. 6I) was expectedly reliant on the presence of NaF channels, while KDR and HCN also had a regulatory role in setting the value of V AP . It should be noted that V AP is significantly dependent on V RMP , as V RMP determines the fraction of sodium channels that are inactivated and are thereby unavailable for activation. As the fraction of available sodium channels is higher with a hyperpolarized V RMP , V AP is higher when V RMP is hyperpolarized. As V RMP was significantly hyperpolarized when HCN channels were knocked out (Fig. 6A), this implied that V AP in HCN channel knockouts should be expected to be higher (Fig. 6I).
Together, analyses of all physiological measurements across valid models using VKMs of each of the 9 active ion channels demonstrated the clear lack of one-to-one mappings between channels and physiological characteristics. Although there was dominance of certain channels in their ability to alter specific measurements, there were several channels that were capable of regulating each measurement and each channel regulated several measurements. Additionally, the effect of virtually knocking out individual channels on all measurements was differential for different channels and measurements, and was variable for even a given channel-measurement combination. The electrophysiological support from LII SCs for several of our conclusions, including the regulatory role of specific channels and the differential and variable dependencies of measurements on channels, with reference to blockade of individual channels is strong (Alonso & Llinas, 1989;Klink & Alonso, 1993;Dickson et al., 2000;Erchova et al., 2004;Nolan et al., 2007;Garden et al., 2008;Pastoll et al., 2012;Boehlen et al., 2013).
These observations together point to a many-to-many configuration of the mapping between channel properties and cellular-scale physiological characteristics, a critical substrate for neurons to exhibit cellular scale degeneracy (Figs. 3-5).

Quantitative predictions: Spike initiation dynamics of stellate cells manifest thetafrequency spectral selectivity and gamma-band coincidence detection capabilities
As we now had a population of models that matched with LII SC physiology, we employed these models to make specific quantitative predictions about critical physiological characteristics of LII SC. Although it is well established that LII stellate cells exhibit spectral selectivity for subthreshold inputs (Haas & White, 2002;Erchova et al., 2004;Giocomo et al., 2007;Garden et al., 2008;Giocomo & Hasselmo, 2008, 2009, it is not known if such frequency selectivity translates to the spike initiation dynamics as well (Das & Narayanan, 2014Das et al., 2017). To assess this, we employed zero-mean Gaussian white noise (GWN) with standard deviation adjusted such that the over all firing rate is ~1 Hz (Fig. 7A). We computed the spike-triggered average (STA) as the mean of all the current stimuli (part of the GWN current, over a 300 ms period preceding each spike) that elicited a spike response in the model under consideration. We derived five distinct measurements of excitability, spectral selectivity and coincidence detection from the STA (Das & Narayanan, 2014Das et al., 2017), and repeated this for all the 155 valid models obtained from MPMOSS ( Fig. 3-5).
The STA computed from the 5 models shown in Fig. 4 showed characteristic class II/III excitability (Fig. 7B), marked by the distinct presence of negative lobes in these STA (Haas & White, 2002;Ermentrout et al., 2007;Haas et al., 2007;Ratte et al., 2013;Das & Narayanan, 2014Das et al., 2017). Employing quantitative metrics from the STA (Fig. 7B-C) for all the 155 valid models, we confirmed that these neurons were endowed with class II/III characteristics. Specifically, our analyses with the valid model population of LII SCs predict that the STA of these cells show theta-frequency spectral selectivity (Fig. 7D, f STA ) with strong selectivity strength (Fig. 7D, Q STA ). As class II/III excitability translates to coincidence detection capabilities in these neurons (Ratte et al., 2013;Das & Narayanan, 2014Das et al., 2017), we computed two distinct measures of coincidence detection window (CDW) from the STA. Whereas the total CDW (T TCDW ) considers the temporal span of the spike-proximal positive lobe of the STA, the effective CDW (T ECDW ) also accounts for the specific shape of the STA in arriving at the CDW (Das & Narayanan, 2014. We computed these CDW measures for all the 155 models, and quantitatively predict that the LII SCs are endowed with gamma-range (25-150 Hz translating to 6.6-40 ms) coincidence detection capabilities (Fig. 7D). Finally, although it is known that the impedance phase of LII SCs manifest a low-frequency inductive lead (Erchova et al., 2004), this inductive phase lead has not been systematically quantified. To do this, we computed the total inductive phase metric (Φ L ; Fig. 7E) developed in  for each of the 155 valid models, and quantitatively predict a prominent inductive phase lead in LII SCs ( Fig.   8F) with specific values for Φ L .

Pair-wise correlations among measurements
Finally, as our analyses (Fig. 6) demonstrated a many-to-many mapping between biophysical parameters and physiological measurements, we asked if there were significant correlations among the distinct measurements. Strong correlations across these measurements (which are reflective of distinct physiological characteristics) would imply that they could be mapped to a smaller set of "core" measurements with the other measurements relegated to redundant and correlated reflections of these core measurements. In addition, strong correlations across measurements would also imply that there are measurements that are strongly reliant on the expression and properties of one specific channel. To assess correlations across measurements, we plotted the pairwise scatter plots (Fig. 8A) spanning all 14 measurements (8 measurements from Fig. 1: Pearson's correlation coefficient for each scatter plot and analyzed the correlation coefficients across measurements (Fig. 8A-B). Although there were strong correlations across some measurements, a majority of these pair-wise correlations were weak (Fig.   8B). Measurements that showed strong pair-wise correlations were those that were known to be critically reliant on specific ion channels. For instance, sag ratio showed strong negative correlation with Φ L , Q STA and Q R , whereas f R manifested strong negative correlation with R in and T TCDW , where all these measurements are known to have a strong dependence on HCN channels (Fig. 6). However, even within this subset of measurements that were strongly dependent on one channel, we noted that only a subset of the pair-wise correlations were high. For instance, with reference to the specific examples discussed above, f R and Q R , both critically dependent on HCN channels and both derived from the impedance amplitude profile, do not show strong pair-wise correlation (Fig. 8A). Thus, although there were a small percentage of measurements that showed strong pair-wise correlations, a majority of their pair-wise correlations were weak, further emphasizing the absence of one-to-one relationships between channels and measurements ( Fig. 6).

DISCUSSION
The prime conclusion of this study is that LII stellate cells of the medial entorhinal cortex express cellular-scale degeneracy in the concomitant expression of several unique physiological characteristics of these neurons. We arrived at this conclusion through an unbiased stochastic search algorithm that spanned 55 parameters associated with biophysically constrained models of active and passive stellate cell components. We validated the outcomes of these stochastically generated models against 10 different physiological characteristics of stellate cells. This validation process provided us with a heterogeneous population of stellate cells that exhibited cellular-scale degeneracy with weak pair-wise correlations across parameters that governed these models. We employed these models to demonstrate the differential and variable dependencies of measurements on underlying channels, and also showed that the mapping between channels and measurements was many-to-many. Finally, we employed this electrophysiologically validated model population to make specific quantitative predictions that point to thetafrequency spectral selectivity and gamma-range coincidence detection capabilities in class II/III spike triggered average of LII SCs.

Correlations between electrophysiological signatures of distal dendrites in CA1 pyramidal neurons and LII MEC stellate cells: Instances of cellular-scale efficient encoding?
A cursory glance at the electrophysiological properties of distal dendrites of CA1 pyramidal neurons and LII MEC stellate cells presents significant correlations between some physiological characteristics of these two structures. Although superficial layers of the MEC project to the distal dendrites of CA1 pyramidal neurons, it is LIII, and not LII, principal neurons of the MEC that project to CA1 pyramidal neurons (Andersen et al., 2006). Despite this, there are several electrophysiological characteristics of LII MEC stellates that match with the distal dendrites of CA1 pyramidal neurons. Several of these similarities are strongly related to heavy expression profile of HCN channels in both these structures. Whereas the gradient in HCN channel density in CA1 pyramidal neurons implies a significantly sharp increase in HCN channels in the distal dendrites of CA1 pyramidal neurons (Magee, 1998;Lorincz et al., 2002), the heavy expression of HCN channels in the cell body of LII MEC stellates is also well established (Klink & Alonso, 1993;Dickson et al., 2000;Erchova et al., 2004;Giocomo et al., 2007;Nolan et al., 2007;Garden et al., 2008;Giocomo & Hasselmo, 2008, 2009Giocomo et al., 2011a;Boehlen et al., 2013).
As a consequence of this, these two structures are endowed with significant sag, similar input resistances and comparable theta-band suthreshold resonance frequencies (Erchova et al., 2004;Pastoll et al., 2012). In addition to these, our study predicts that the STA of LII SCs should be endowed with theta-frequency spectral selectivity and gammaband coincidence detection windows. Although it is known that the STA of LII SCs manifest class II/III excitability with coincidence detection characteristics (Haas & White, 2002;Haas et al., 2007), quantification of spectral selectivity in the STA or a systematic assessment of the specific frequency band of coincidence detection capabilities has not been assessed. Our study specifically predicts the coincidence detection window (Fig. 7D; T ECDW ) associated with the STA of LII SCs to be in the fastgamma frequency (60-120 Hz) range, with a high theta-range spectral selectivity in the STA ( Fig. 7D; f STA ). Interestingly, quantitative predictions for these measurements for the distal dendritic region of CA1 pyramidal neurons also fall within the same spectral bands (Das & Narayanan, 2015).
This confluence of STA measurements of distal dendrites in CA1 pyramidal neurons and of LII SC soma, especially that of the coincidence detection window falling within the fast gamma-frequency band, is striking from the perspective of gamma-band multiplexing that has been reported in the CA1 subregion (Colgin et al., 2009;Colgin & Moser, 2010;Fernandez-Ruiz et al., 2017). Specifically, gamma oscillations in the CA1 subregion have been shown to manifest stratification into distinct fast and slow frequency components, impinging respectively on distal and proximal dendritic regions. This spatially stratified frequency-division multiplexing allows for differential coupling of CA1 pyramidal neurons to afferent inputs from the MEC and CA3 through different gamma bands. Juxtaposed against this, and within the efficient coding framework (Barlow, 1961;Bell & Sejnowski, 1997;Simoncelli & Olshausen, 2001;Lewicki, 2002;Simoncelli, 2003) where the response filters in a single neuron should match the natural statistics of afferent network activity (Narayanan & Johnston, 2012), it may be argued that different dendritic locations should be equipped with filters (STA kernels) that are matched to the afferent inputs (different gamma frequencies) that are received by that specific location. As a specific instance of such location-dependent efficient encoding of afferent network statistics in hippocampal pyramidal neurons, it has been quantitatively postulated that the distal dendrites of CA1 pyramidal neurons are endowed with coincidence detection windows specific to fast-gamma frequencies, whereas those of the proximal dendrites matched with slow-gamma frequencies (Das & Narayanan, 2015).
In this context and given the several lines of evidence for the dominance of fast gamma oscillations in the superficial layers of MEC (Colgin et al., 2009;Colgin, 2016;Trimper et al., 2017), our predictions that the CDW for stellate neurons should be in fast gamma band point to a similar form of efficient encoding schema in the MEC.
Specifically, if the fast gamma oscillations are statistically the most prevalent in the superficial layers of MEC, it is imperative that neurons there are equipped with the machinery that is capable of detecting and processing inputs in this frequency range.
Additionally, from the efficient coding perspective, as neurons tune their response properties to efficiently represent the statistics of afferent inputs (Wiesel & Hubel, 1963b, a;Hirsch & Spinelli, 1970;Bell & Sejnowski, 1997;deCharms et al., 1998;Kilgard & Merzenich, 1998;Stemmler & Koch, 1999;Sharma et al., 2000;Kilgard et al., 2001;Simoncelli & Olshausen, 2001;Lewicki, 2002;Simoncelli, 2003;de Villers-Sidani et al., 2007), it is essential that the response properties of LII MEC neurons are tuned to the fast gamma frequency range. Thus, our prediction on a fast gamma-band CDW in the STA of LII EC cells suggests the possibility of efficient encoding spanning the hippocampal formation, whereby the neuronal properties in terms of their class of excitability and specific band of frequency where their coincidence detection windows lie match with the type of gamma-frequency band that is most prevalent in that subregion (or strata in case of the CA1). A direct test of this experimental prediction would be to measure the CDW of pyramidal neurons in the CA3, of different dendritic subregions of the CA1 and of LII MEC stellates and ask if these CDW match with the respective gamma-band inputs that are prevalent in these subregions.
Finally, encoding schema are state-dependent processes that are critically reliant on behavioral state and consequent changes in afferent activity, neuromodulatory tones and activity-dependent plasticity (Marder, 2012;Narayanan & Johnston, 2012;Bargmann & Marder, 2013;Ratte et al., 2013;Srikanth & Narayanan, 2015;Das et al., 2017;Deneve et al., 2017;Gallistel, 2017). Therefore, it is important that postulates on efficient codes and relationships between neuronal activity and afferent statistics are assessed in a manner that accounts for adaptability of coding within the neuron and across the network.
Such activity-dependent plasticity and neuromodulation of intrinsic properties, especially of signature characteristics such as the membrane potential oscillations, could be systematically assessed in entorhinal stellates. Specifically, as activity-dependent plasticity of several ion channels is well established across different brain regions (Magee & Johnston, 2005;Narayanan & Johnston, 2007;Sjostrom et al., 2008;Narayanan et al., 2010;Shah et al., 2010;Narayanan & Johnston, 2012), it would be important to ask if membrane potential oscillations, spectral selectivity characteristics (f R , Q R , f STA , Q STA , Φ L ) and coincidence detection windows are amenable to such activity-dependent plasticity that target different ion channels (Fig. 6).

Implications of cellular-scale degeneracy in LII stellate cell physiology
Degeneracy, the ability of disparate structural components to elicit similar function, is a ubiquitous biological phenomenon with strong links to robust physiology and evolution (Edelman & Gally, 2001;Price & Friston, 2002;Leonardo, 2005;Whitacre & Bender, 2010;Whitacre, 2010). Several studies spanning different systems have now demonstrated the expression of degeneracy, at different scales of analysis in neural systems as well (Marder & Prinz, 2002;Prinz et al., 2004;Marder & Goaillard, 2006;Marder, 2011;Marder & Taylor, 2011;O'Leary & Marder, 2014;O'Leary et al., 2014;Vogelstein et al., 2014;Drion et al., 2015). Within the hippocampal formation, recent studies have demonstrated the expression of degeneracy in single-neuron electrophysiology (Rathour & Narayanan, 2012;Srikanth & Narayanan, 2015;Mishra & Narayanan, 2017), intraneuronal functional maps (Rathour & Narayanan, 2014;Rathour et al., 2016), synaptic localization required for sharp-tuning of hippocampal place fields First, the many-to-many mapping between channels and physiological characteristics ( Fig. 6) and the consequent degeneracy in concomitantly achieving all signature electrophysiological characteristics implies that there is no explicit necessity for maintaining individual channels at specific levels or for maintaining paired expression between channel combinations . This provides several significant degrees of freedom to the protein localization and targeting machinery in achieving specific functions or in maintaining homeostasis in these functional characteristics. Second, this also implies that adaptability to external stimuli, in terms of achieving efficient codes of afferent stimuli or in encoding features of a novel stimulus structure, could be achieved through disparate combinations of plasticity in several constituent components. For instance, our results predict that the ability to achieve fast gamma-band coincidence detection capabilities could be achieved through distinct combinations of channel parameters (Figs. 4-5, Fig. 7). Experimental analyses of such degeneracy in achieving efficient matching of neuronal response characteristics with the statistics of oscillatory patterns, through the use of distinct pharmacological agents that target different channels, would demonstrate the ability of different channels to regulate such efficient encoding . Finally, degeneracy in the expression of excitability properties also implies that similar long-term plasticity profiles in these neurons could be achieved with disparate combinations of parameters. Specifically, several forms of neuronal plasticity are critically reliant on the amplitude and kinetics of cytosolic calcium entry, which in turn are dependent on neuronal excitability properties. As similar excitability profiles could be achieved with distinct combinations of constituent components, it stands to reason that similar plasticity profiles could be achieved through disparate parameter combinations (Narayanan & Johnston, 2010;Ashhad & Narayanan, 2013;Anirudhan & Narayanan, 2015).

Future directions: Electrophysiological and computational
A critical future direction for the study presented here is the incorporation of biological heterogeneities into entorhinal network models that assess grid cell formation. Most At the cellular scale, future studies could build heterogeneous models to account for the signature continuum of intrinsic physiological characteristics along the dorsoventral axis of the MEC, also accounting for specific channels that are known to change along this axis (Giocomo et al., 2007;Garden et al., 2008;Giocomo & Hasselmo, 2008, 2009Yoshida et al., 2011). Such studies will provide quantitative bases for exploring the expression of degeneracy in maintaining the dorsoventral gradients, and could be incorporated into network models for grid formation in assessing the relationship between grid-cell characteristic and neuronal intrinsic properties.
Additionally, the heterogeneous model population built in this study comprises a simple single-compartmental structure that did not account for dendritic properties or morphological heterogeneity of LII SCs. Electrophysiologically, the absence of systematic cell-attached recordings of specific channels and their properties in the stellate cells has been a significant impediment in building morphologically realistic models.
While future experimental studies could focus on recording channels and channel properties along the non-planar dendritic arbor of stellate cells, future computational studies could incorporate these channels into morphologically realistic models to assess the specific roles of dendritic channels and morphological heterogeneity in grid cell formation (Schmidt-Hieber et al., 2017). Finally, our study has specific quantitative predictions about the STA of these neurons and also presents an array of crossdependencies of measurements on different channel types. Future electrophysiological studies could systematically test these predictions, and assess efficient encoding in these structures apart from adding further evidence for the many-to-many mapping between channels and physiological characteristics. For instance, an important prediction from VKMs is on the critical role of SK channels in regulating several intrinsic measurements including membrane potential oscillations (Fig. 6). Although the expression of calciumdependent potassium channels is established in stellate cells (Khawaja et al., 2007;Pastoll et al., 2012), the specific role of these channels in regulating resonance, impedance phase, intrinsic excitability and membrane potential oscillations could be tested electrophysiologically using pharmacological blockers of SK channels.     (Table 1). Individual models are color coded across all plots and traces in A-I. Matrix consisting the pair-wise scatter plots between the 55 parameters (Table 1) for all 155 valid stellate cell models. Histograms in the bottom row depict the span of the corresponding parameter with reference to their respective min-max ranges (Table 1)