Using a Semi-Automated Strategy to Develop Multi-Compartment Models That Predict Biophysical Properties of Interneuron-Specific 3 (IS3) Cells in Hippocampus

Abstract Determining how intrinsic cellular properties govern and modulate neuronal input–output processing is a critical endeavor for understanding microcircuit functions in the brain. However, lack of cellular specifics and nonlinear interactions prevent experiments alone from achieving this. Building and using cellular models is essential in these efforts. We focus on uncovering the intrinsic properties of mus musculus hippocampal type 3 interneuron-specific (IS3) cells, a cell type that makes GABAergic synapses onto specific interneuron types, but not pyramidal cells. While IS3 cell morphology and synaptic output have been examined, their voltage-gated ion channel profile and distribution remain unknown. We combined whole-cell patch-clamp recordings and two-photon dendritic calcium imaging to examine IS3 cell membrane and dendritic properties. Using these data as a target reference, we developed a semi-automated strategy to obtain multi-compartment models for a cell type with unknown intrinsic properties. Our approach is based on generating populations of models to capture determined features of the experimental data, each of which possesses unique combinations of channel types and conductance values. From these populations, we chose models that most closely resembled the experimental data. We used these models to examine the impact of specific ion channel combinations on spike generation. Our models predict that fast delayed rectifier currents should be present in soma and proximal dendrites, and this is confirmed using immunohistochemistry. Further, without A-type potassium currents in the dendrites, spike generation is facilitated at more distal synaptic input locations. Our models will help to determine the functional role of IS3 cells in hippocampal microcircuits.


Introduction
Uncovering the intrinsic cellular properties that control the input-output processing within brain microcircuits is highly challenging due to the large diversity of voltagegated ion channels (VGCs) expressed by individual neurons and their variable distribution across the cell structure. In particular, the spatial location of VGCs determines the integration of synaptic inputs by specific subcellular regions. Because of the difficulty in investigating subcellular properties with purely experimental means, multi-compartmental modeling has been used as a key tool to explore the contributions of different VGC types (Carnevale and Hines, 2006).
Traditionally, morphologically detailed multi-compartment models are developed using hand-tuning methods (De Schutter and Bower, 1994;Migliore et al., 1995Migliore et al., , 1999Saraga et al., 2003;Lawrence et al., 2006;Hu et al., 2010), where passive and active parameter values are adjusted manually until the model reproduces experimentally observed electrophysiology characteristics (e.g., spike frequency). Recognizing that similar outputs can be achieved with different sets of parameter values, along with inherent biological variability, model database approaches have been developed (Prinz et al., 2003;Günay et al., 2008;Hay et al., 2011;Sekulić et al., 2014). These approaches typically generate populations of models by varying VGC maximal conductances and extracting a subset of the models that best match the experimental data. Often, however, information on the types of VGCs and their biophysical details is not available. While having multi-compartment models to determine the roles of particular cell types in microcircuit dynamics is useful, how best to proceed in building multicompartment models with limited experimental data is unclear. Furthermore, some sort of cyclic approach between model and experiment is clearly warranted, as models require continuous updating as more experimental data be-comes available. In particular, Sekulić et al. (2014) designed a cyclic model database approach where the database design hinged on a particular question. In this way, they not only obtained populations of models representing the particular cell type for subsequent examination, but they also motivated specific experimental studies for future consideration. However, their starting point requires that a "reference" multi-compartment model with known densities and distributions of VGCs that can capture the experimental data already exists. This would typically have arisen from previous hand-tuning work. Given that not all cell types possess such multi-compartment model representations, and that hand tuning is not the most efficient method for investigating parameter spaces, developing other approaches would be useful.
Here, we devised a strategy that uses populations of models in conjunction with hand tuning to generate multicompartment models with characteristics that match features from electrophysiological recordings. We chose not to fully automate the process to allow flexibility in the model-building process, when no particular VGC data are available. At the same time, the partial automation through use of populations of models reduced the subjectivity and tediousness of pure hand tuning. We applied our strategy to obtain reference multi-compartment models for an identified cell type, the CA1 type 3 interneuron-specific (IS3) cell in the hippocampus. The models provided suggestions of types, densities and distributions of VGCs in IS3 cells.
The IS3 cell is a type of inhibitory interneuron that exclusively targets other interneurons and does not contact pyramidal cells (Acsády et al., 1996a,b;Gulyás et al., 1996). IS3 cells may receive excitatory inputs from the entorhinal cortex via the perforant path, from the CA3 area via Schaffer collaterals (SCs) and from CA1 local collaterals (Francavilla et al., 2015). Through photostimulation experiments and paired recordings, it has been shown that IS3 cells provide a major local source of inhibition to oriens-lacunosum moleculare (OLM) cells, and optimally control their firing at theta frequencies Tyan et al., 2014). OLM cells in turn are known to provide a gating function in the information flow from the perforant path to the hippocampus (Maccaferri and McBain, 1995;Leão et al., 2012). We used our developed multi-compartment models with different combinations of VGCs to investigate how the type and location of dendritic conductances in conjunction with synaptic input influenced the IS3 output.
Our strategy for building multi-compartment models can be considered for other cell types where no VGC information is available. The developed models can then provide guidance to experimentalists in narrowing down which biological details are both plausible and important in microcircuit function. Moreover, the models can be used as reference models for subsequent examination. In this study, our developed multi-compartment IS3 cell models with fast delayed rectifier potassium channels distributed in the soma (S) and proximal dendrites (D) can capture the experimental data.

Two-photon calcium imaging
Dendritic Ca 2ϩ imaging was performed using a TCS SP5 two-photon laser-scanning microscope (Leica Microsystems) based on a Ti-sapphire laser (Chameleon Ultra II, Coherent; Ͼ3 W, 140 fs pulses, 80 Hz repetition rate) tuned to 800 nm. A long-range water-immersion objective [40ϫ; numerical aperture (NA), 0.8] was used to collect photons in the epifluorescence mode with external photomultiplier tubes. Neurons were filled with Oregon Green ® 488 BAPTA-1 (OGB-1) via the patch electrode for 20 -30 min before imaging. Red fluorescence from Alexa Fluor 594 was used to locate dendrites of interest. To measure Ca 2ϩ signals, green fluorescence was collected during 400 Hz line scans across dendritic segments of 2-15 m. Fluorescence changes were quantified as changes in green fluorescence from the baseline (i.e., (F Ϫ F 0 )/F 0 ), referred to as ⌬F/F in Figure 1C-E.

Morphological reconstructions
For morphological reconstructions of IS3 cells, recorded neurons were filled with biocytin (Sigma-Aldrich) during whole-cell recordings. Slices with recorded cells were fixed with 4% paraformaldehyde (PFA), rinsed in TBS, and incubated with Streptavidin-Alexa Fluor 488 overnight at ϩ4°C, following which the slices were mounted. Confocal z-stacks were acquired using a Leica SP5 microscope (Leica Microsystems) with a 1 m step. Morphological reconstructions were performed using the Neurolucida 8.26.2 software (MBF Bioscience).

Neurochemical analysis
For neurochemical analysis, animals were perfused with 4% PFA and hippocampal sections (thickness, 30 m) were prepared using a vibratome (VT1000, Leica Microsystems). Sections were permeabilized with 0.1% Triton X-100 in TBS and incubated in blocking solution containing 20% normal goat serum for 1 h. After this step, sections were incubated with anti-GFP (chicken, 1:1000; Aves Labs) and either anti-KCNB1 (mouse, 1:1000; Sigma-Aldrich) or anti-KCNC1 (rabbit, 1:500; Sigma-Aldrich) primary antibodies at room temperature for 24 -48 h. For KCNC1, slices were incubated the following day with biotinylated anti-rabbit (goat, 1:200; Vector Laboratories) antibodies for 4 h at room temperature, then for 24 h with A/B reagents (1:100; Vectastain ABC Kit, Vector Laboratories) at room temperature. The following day, slices were incubated for 2 h with the following conjugated secondary antibodies: anti-chicken Alexa Fluor 488 (1:1000; Jackson ImmunoResearch); and either antistreptavidin Alexa Fluor 546 (1:200; Jackson ImmunoResearch) or Alexa Fluor 647 (1:250; Invitrogen). Slices were then rinsed and mounted on microscope slides. Confocal images of labeled sections were obtained using a Leica TCS SP5 imaging system equipped with a 488 nm argon, a 543 nm HeNe, and a 633 nm HeNe laser, with a 63ϫ (NA, 1.4) oil-immersion objective (Leica Microsystems). Based on previous observations (Tyan et al., 2014), GFPexpressing bipolar cells with a small size cell body (10 -15 m) and vertically oriented dendrites were considered as putative IS3 cells, and were used for analysis of the voltage-gated K ϩ channel subunit (Kv) expression.

Simulation and analysis software
All simulations were performed using the NEURON software environment (Carnevale and Hines, 2006). In some instances, these simulations were computed using the Neuroscience Gateway for high-performance computing (Sivagnanam et al., 2013). All simulations as well as experimental recordings were analyzed using a MATLAB toolbox called PANDORA (Günay et al., 2009) as well as customized MATLAB code. Models will be made publicly available on ModelDB (Hines et al., 2004) following publication of this article.

Model morphology and compartmentalization
Topological information provided by the morphological reconstruction of an IS3 cell was imported directly into NEURON. The morphological structure was then subdivided into multiple compartments until the point where the compartments were small enough to be considered isopotential. The topological model is shown in Figure 1A. Because inclusion of the full axon in the model requires more specifics regarding biophysical axonal properties (i.e., channels types and distributions) and more computational strain, we considered a reduced morphological version where the axon is cut down to a small remaining section (M2) relative to the full version with the axon intact (M1; Fig. 1A). In the M2 version, there were 221 compartments, compared with 653 in the full M1 version, and the surface area was reduced from 20,871 to 7297 m 2 (ϳ35% of the original surface area). While both morphologies were considered, we focused on the M2 version.

Model passive properties
The passive parameters of the model were adjusted such that resting membrane potential, input resistance, and membrane time constant were matched with their In cell 1, current step injections show hyperpolarization with minimal sag (Ϫ100 pA; red), passive response without spiking (ϩ10 pA; green), and depolarization with irregular spiking (ϩ50 pA; blue). In cell 2, current injections show hyperpolarization with minimal sag (Ϫ100 pA; red), depolarization with spiking (ϩ50 pA; blue), and depolarization block (ϩ500 pA; orange). Recordings were obtained in the presence of synaptic blockers (i.e., NBQX, AP5, and gabazine). C, Two-photon image of an IS3 cell filled with Alexa Fluor 594 and OGB-1. White lines indicate where backpropagating action potential-evoked Ca 2ϩ -transient (bAP-CaT) line scans were performed. The protocol consisted of three consecutive 2 ms, 800 pA somatic current step injections in conjunction with dendritic CaT recordings. D, Representative examples of dendritic bAP-CaTs evoked by three APs at the soma. E, Summary plot showing average changes (solid red) and individual changes (dotted gray) in the bAP-CaT amplitude at different distances from the soma. The dashed black line indicates the threshold level below which the calcium signal was indistinguishable from noise.
New Research appropriate experimentally observed values. The passive properties of the model were matched with raw data from the respective cell, and not with the experimental averages. The specific membrane capacitance and axial resistance parameters were restricted to appropriate experimentally estimated ranges. The average measurement of specific membrane capacitance in cortical pyramidal neurons, spinal cord neurons, and hippocampal neurons, including interneurons, is 0.9 F/cm 2 (Chitwood et al., 1999;Gentet et al., 2000), suggesting that the range for this parameter is fairly consistent across many neuron types. Values for axial resistance in individual neurons can range from 50 to 400 ⍀/cm (Dayan and Abbott, 2005;Holmes, 2010).
Furthermore, in the case of the M2 morphology, passive parameter values in the soma and dendrites were kept the same as the passive parameter values in the M1 morphology, but the passive parameter values in the remaining axon sections were readjusted in order to get appropriate resting membrane potential, input resistance, and membrane time constant values. This adjustment was performed to compensate for the axonal surface area difference and is a strategy that has been used previously (Lawrence et al., 2006). Passive properties for M1, M2, and experimental data are given in Table 1.

Electrotonic analysis in passive model
Using the passive model and the Electrotonic Analysis toolbox in NEURON, we looked at the natural log of voltage attenuation across the morphology of the model, with the soma as the electrode reference point. In this case, the natural log of voltage attenuation is considered the electrotonic distance (i.e., L ϭ log[A]), and attenuation (A) is measured as: A ϭ voltage upstream/voltage downstream (i.e., where voltage upstream is an applied 1 mV signal and voltage downstream is the downstream response to the 1 mV signal). This definition of electrotonic distance allows for a direct relationship to attenuation, regardless of the cellular morphology. For this analysis, both voltage spreading toward the soma (V in ) as well as away from the soma (V out ) were plotted.

Semi-automated strategy
To determine VGC types, densities, and distributions in IS3 cells, we used an approach that combined hand tuning with automation, allowing both automated simulations and examination of multiple features to simultaneously capture various modes of IS3 cell dynamics. The steps in this strategy (illustrated in Fig. 2) are outlined here, where each step is prefaced by whether it is an automated (Auto) or manual (Hand) step:  Fig. 1B). Experimental traces with these features were selected and they had CIPs of Ϫ100, ϩ20, ϩ50, and ϩ500 pA ( Table 2, first column). Additionally, an initial set of VGC types needs to be chosen and hand tuning used to find a single set of parameter values that approximately matches electrophysiological features. These parameter values are then used as a basis to generate a population of models. We chose three to four parameter values that encompassed each of the initial hand-tuned values (the maximal conductance values were used; see Table 6). 1.
[Auto] Generate model database: a population of models was generated in the NEURON software environment where each model possessed a unique combination of channel conductance values. The CIPs were applied to each model, and the voltage traces generated from each model were imported into MAT-LAB using PANDORA (Günay et al., 2009) and organized into databases. 2.
[Auto] Eliminate inappropriate models: the models that did not capture IS3 cell features at the given CIPs were eliminated. The criteria used were based on an examination of the experimental data (see Results and Table 2, fifth column). 3.
[Auto] Compute model measurements: for the remaining models, characteristic measurements for each feature were determined (Table 2, third column). Specifically, each CIP possessed its own set of characteristic measurements (e.g., spike amplitude mean at ϩ50 pA, potential sag at Ϫ100 pA). For Ϫ100 pA, 5 measurements were evaluated; for ϩ20 pA, 1 measurement was evaluated; for ϩ50 pA, 10 measurements were evaluated; and for ϩ500 pA, 2 measurements were evaluated. Although computation of most of these measurements already exists in the PANDORA toolbox, the following two customized measurements needed to be added: membrane potential (V m ) difference and membrane potential difference in the last 700 ms of the CIP. 4.
[Auto] Generate a quality metric for each model: each measurement for each model was compared with the same measurement from the selected experimental traces representing the signature feature. Dashed lines in Figure 3 histograms show the measurements obtained from the selected traces (Table 2, fourth column) used in comparison with the model measurements. The normalized Euclidean distance (Günay et al., 2009) between each model and selected experimental trace measurements was then computed as follows: For each CIP step (1st column) we note a signature feature that is observed experimentally (2nd column), characteristic measurements that are used to quantify the signature feature (3rd column), characteristic measurement values from the selected IS3 traces (4th column), as well as conditions where models would be rejected for clearly not capturing the experimental feature (5th column). CIP, Current Injection Protocol; V m , membrane potential. In this case, x i and y i represent the ith measurement ( Table 2, third column) of N total measurements for the simulation (x) and experiment (y; Fig. 3, dashed lines), respectively. Since d x,y was generated using only one experimental trace for each CIP, the experimental measurement SD i was set to 1 for each measure and, therefore, had no effect on the normalization. The goal here is to use a metric to sort the remaining models relative to the selected experimental traces and so obtain representative models. In this sense, d x,y serves as a quality metric. The normalized Euclidean distance was evaluated for each CIP, and then the four CIP normalized Euclidean distances were averaged together so that equal importance was given to each signature feature. However, since there are a different number of measurements for each CIP, each measurement did not have equal importance. We felt that this was reasonable to do at this stage, given that it is unknown which specific dynamic regimes are more or less important for IS3 cell function in hippocampal microcircuits. 5.
[Auto] Find representative models: the models were ranked using the quality metric from smallest ("best") normalized Euclidean distance to largest. This allowed us to obtain ranked, representative models for each database, albeit qualitative in nature. 6.
[Hand] Adjust parameter space: the model parameter space was adjusted to reduce the number of models that were eliminated in step 2 of the following iteration (i.e., in cycling back to step 1). This hand-tuning step was performed using clutter-based dimensional reordering (CBDR) plots, a method used previously by Taylor et al. (2006). CBDR plots allow a visualization of the parameter space such that it is easy to see when models in the database are eliminated. Further, these plots visualize the structure of the conductance parameter space by reordering the channel conductance parameters along the N-dimensional tensor in order to sort higher-order channels from lower-order channels (i.e., channels that have a greater or smaller impact, respectively, on normalized Euclidean distance). In other words, each pixel in these plots represents a different model with a unique set of VGC conductance values, whose value can be extrapolated from the scale bars on the x-and y-axes using the maximal conductance values found in Table 6. To minimize space in the CBDR plots where models had been eliminated ( Fig. 4B1-B3, black pixels) we adjusted the parameter values (in range and resolution) to use in generating the next database (step 1) focus-ing only on conductance parameter spaces that had generated good models. A cycling process of this is shown in Figure 5. In some instances, we avoided parameter sets for particular reasons (e.g., lowamplitude spikes caused by low transient sodium conductances in SDprox.1; Fig. 5, Cycle 2) even though they were not eliminated in step 2. Once it was no longer possible to increase the parameter space from the inspection of CBDR plots, we moved to step 7. 7.
[Hand] Adjust channel types: the model channel distributions and/or channel types were adjusted to improve the ability of the models to capture signature experimental features. That is, we analyzed different cases where each case possessed different distributions of channels or had added, altered, or removed channel types, and decided whether there were improvements ( Fig. 4C). Once manual adjustments were made, we returned to step 1.
The overall process is terminated when it is no longer possible improve on the quality metric and the resulting models capture the known experimental data. As these last steps are manual, this termination point is qualitatively determined.

Uniform channel distributions in soma and proximal dendrites
To incorporate dendritic VGC distributions consistent with the calcium imaging data, we specified channel conductance values that were uniform across the soma and proximal dendrites. We considered this as reasonable since there are no data to support particular distributions (e.g., exponential). A Boltzmann function was used to describe the relationship between the channel conductance and distance from the soma such that non-zero conductances were present only in the soma and proximal dendrites (so, overall, a non-uniform distribution), as follows: where f(p) is the distance-dependent conductance, p is the distance along the dendrite (in m), G is the desired conductance value in the proximal dendrite (in S/cm 2 ), and k is a value of 10 (in m Ϫ1 ) in order to allow a sharp decrease in conductance at around d, where a d of 55, 75, or 95 ensures that the decrease in conductance occurs at around 50, 70 or 90 m from the soma. Thus, this equation describes the relationship between channel conductance and distance from soma, and the parameters are continued hyperpolarization (A-E; CIPs include Ϫ100 and Ϫ90 pA for the "lack of sag feature"), passive depolarization (F; CIPs include 10 and 20 pA for the "no spiking feature"), depolarization with spiking (G-P; CIPs range from 10 to 140 pA with 10 pA intervals for the "normal spiking feature"), or depolarization block (Q, R; CIPs range from 450 to 700 pA with 50 pA intervals for the "depolarization block feature"). Note that the dashed lines indicate the measurements obtained from the selected IS3 cell experimental trace used to compute the distance metric for each model. The signature features, and characteristic measurements and their values are given in Table 2 (second, third, and fourth columns).
New Research tuned in such a way that the conductance is uniform at a value of G up to a distance of 50, 70, or 90 m from the soma, after which it drops sharply to 0 S/cm 2 .

Analysis of the input-output relationship in IS3 cell models
To understand how the intrinsic properties could affect spike generation in IS3 cells, models were analyzed regarding the required input for spike generation. This was done by defining a synaptic conductance that is activated due to a single presynaptic spike to a target point along the dendritic tree of the model. In this case, a point process is used to describe the synapse as a two-state kinetic scheme to produce synaptic current according to the following: where i is the synaptic current (in nA), G is the synaptic conductance (in S), v is the membrane potential (in mV), e is the reversal potential (set to 0 mV), weight is synaptic weight (in S), factor is a NEURON process used to The current injection protocol is Ϫ100, ϩ20, ϩ50, and ϩ500 pA. B1-B3, Parameter spaces for S.1, S.2, and SD (from left to right), as visualized using CBDR. Each pixel represents a single model and the distance metric of that model. Conductance axes are organized such that overall low-conductance models are in the bottom left quadrants and overall high-conductance models are in the top right quadrants. Black pixels represent models that are eliminated in step 2 and are assigned a distance value of 100 as a result. For S.1, S.2, and SD, respectively, 182 of 432, 59 of 81, and 42 of 81 models did not get rejected (i.e., the remaining colored pixels). C, Example experimental and model database (i.e., for S.1, S.2, and SD) voltage trace measurement histograms for spike half-width and first spike time during 50 pA stimulation. The additional flowchart above the figures shows example rationales for the changes that were made among the S.1, S.2, and SD distributions, including reasoning for why the S.2 models best captured IS3 experimental features. normalize the peak synaptic conductance to 1 (i.e., it ensures that the peak of G in response to a single presynaptic spike is equal to weight), t is time (ms), tau1 is the rise time (set to 0.2 ms), and tau2 is the decay time (set to 2 ms). Although rise time and decay time can vary across different cell types, the values used were chosen based on the measurements in other interneuron types in hippocampus (Tóth, 2010). A full examination of synaptic parameters and their effects on spike generation was beyond the scope of this study. The synaptic model was incrementally moved to different locations along the dendrite. At each location, the minimum weight value necessary to evoke an action potential at the soma was numerically determined by increasing the weight in increments until the somatic membrane potential surpassed Ϫ20 mV in response to a single presynaptic spike. After a somatic spike was recorded, the synapse was then moved to a further dendritic location, and the process was repeated until the entire dendritic arbor was analyzed. Importantly, increasing the weight of a single synapse would eventually lead to a saturation point where (v -e) ϭ 0 because v would become so depolarized that it was essentially equal to the synaptic reversal potential (e ϭ 0 mV) at the dendritic site of the synapse. At this saturation point, no additional current could be generated, regardless of whether weight was further increased (i.e., because i ϭ G‫(ء‬v -e) ϭ G‫0ء‬ ϭ 0 pA). To track this, we recorded the change in the membrane potential at the site of the synapse in the first 1 ms following the synaptic event.

Results
The morphological and membrane properties of IS3 cells have been reported Tyan et al., 2014); however, the diversity and spatial distribution of VGCs in these cells remain unknown. To simulate how specific types of VGCs located at distinct subcellular domains can affect the IS3 cell input-output properties, we developed a semi-automated strategy (Fig. 2), which generated models mimicking the voltage dynamics seen in IS3 cells. This required the following: (1) experimentally obtained electrophysiological signature features to capture; (2) a compartmentalized model morphology with appropriate passive properties; and (3) initial choices for VGCs.

Experiment: electrophysiological features of IS3 cells
To use our strategy to estimate IS3 cell VGC composition, we first determined the electrophysiological features that should be captured. For this, we examined voltage traces obtained in response to different current steps, a CIP, in the presence of synaptic blockers to explore IS3 intrinsic cell membrane properties. We note that previous work has shown that CIPs are a simple protocol that sufficiently exposes the dynamics of a cell (Druckmann et al., 2011). We observed the following.
First, during hyperpolarizing steps, IS3 cells exhibited a small membrane potential sag (Fig. 1B, red traces). With depolarizing steps, IS3 cells first exhibited a passive response ( Fig. 1B, green trace) and then began spiking (Fig.  1B, blue trace, cell 2). The data obtained here is consistent with a previously recorded rheobase in these cells of 42.8 Ϯ 8 pA (Tyan et al., 2014). Also, these cells exhibited irregular firing patterns, with occasional demonstrations of regular adaptive firing (Fig. 1B, blue traces). At higher depolarization steps, progressive spike amplitude adaptation was observed (Fig. 1B, blue trace, cell 2), with depolarization-dependent spike block at current steps exceeding 250 pA (Fig. 1B, cell 2, orange trace).
With this set of experimental data in mind, we made decisions on which "signature features" to capture in our models. Although there is variability in the electrophysiological data, we needed to make assumptions on how best to constrain the models without being overly restrictive or adding excessive detail that would unnecessarily add computational strain. Thus, we chose the following four signature features that occurred for specific CIP ranges: lack of sag during hyperpolarization; depolarization without spiking; depolarization with normal spiking; and depolarization block (Table 2, second column). Examples of all these signature features are shown in Figure  1B. As subthreshold membrane potential fluctuations in the theta range similar to those described in other types of hippocampal interneurons (Chapman and Lacaille, 1999;Morin et al., 2010) were observed in IS3 cells, the VGC types that could be involved were considered when making the initial choices of VGCs. Subthreshold fluctuations and irregular firing were included in the models by the generic addition of noise but were not part of the semiautomated strategy.
For each of the chosen features, we identified characteristic measurements (Table 2, third column). The experimental data are summarized in Figure 3, where each panel is a histogram of a particular measurement from all the experimental traces showing signature features. We selected experimental traces, and the particular CIP values and measurements for them are shown in Table 2 (first and third columns; Fig. 3, histograms, dashed lines). Trace selection was performed manually by going through the experimental data and identifying key electrophysiological regimes. While one might consider automating this continued SDprox.1, five cycles of steps 1 to 6 were required (Cycle 1 ϭ 47 remaining of 243 models; Cycle 2 ϭ 64 remaining of 243 models; Cycle 3 ϭ 22 remaining of 243 models; Cycle 4 ϭ 52 remaining of 243 models; Cycle 5 ϭ 70 remaining of 243 models). B, for SDprox.2, three cycles of steps 1 to 6 were required (Cycle 1 ϭ 24 remaining of 243 models; Cycle 2 ϭ 107 remaining of 243 models; Cycle 3 ϭ 146 remaining of 243 models). The CBDR plots show the quality of the parameter space of the model database during each cycle. Note that blue-boxed areas indicate parameter spaces of interest that were focused on. Red-boxed areas indicate parameter spaces that were purposely avoided. aspect, because of biological variability and not knowing the extent of IS3 cell dynamics that might be functionally important, it would be challenging to develop appropriate automated algorithms. Overall, we felt that this was a balanced strategy that would allow us to capture salient IS3 cell features. At the same time, it was necessary to choose experimental traces to be able to automate steps in the process (Fig. 2).

Experiment: action potential evoked dendritic calcium signals in IS3 cells
Experiments with dendritic patch-clamp recordings and two-photon Ca 2ϩ imaging of dendritic Ca 2ϩ elevations evoked by backpropagating action potentials (bAPs) have revealed that interneuron types that have a high density of VGCs in their dendrites (in particular, of sodium channels) exhibit non-decrementing bAP-evoked Ca 2ϩ signals (Martina et al., 2000;Topolnik et al., 2009;Hu et al., 2010;Camiré and Topolnik, 2014). As direct dendritic recordings from IS3 cells are technically challenging, we performed dendritic two-photon Ca 2ϩ imaging of bAPevoked Ca 2ϩ transients (bAP-CaTs) in these cells (Fig. 1C-E; Table 3). bAP-CaTs were evoked by the somatic injection of three consecutive current pulses (800 pA, 2 ms), and the resulting CaTs were recorded at different points along the dendrite. Our data showed that bAP-CaTs had the same amplitude up until 50 m from the soma (soma ⌬F/F, 0.49 Ϯ 0.08; 50 m ⌬F/F, 0.53 Ϯ 0.06; n ϭ 9, p Ͼ 0.05, signed rank test) and declined gradually at 150 m from the soma (⌬F/F, 0.22 Ϯ 0.08; n ϭ 7, p Ͻ 0.05 compared with somatic CaT, signed rank test). Although the data from individual cells were variable, all cells followed the same trend, exhibiting a significant decline in bAP-CaT amplitude between 50 and 150 m from the soma (50 m ⌬F/F, 0.53 Ϯ 0.06; 150 m ⌬F/F, 0.22 Ϯ 0.08; n ϭ 7, p Ͻ 0.05, signed rank test). These data are consistent with previous observations in different interneuron subtypes (Martina et al., 2000;Golding et al., 2001;Topolnik et al., 2009;Hu et al., 2010;Camiré and Topolnik, 2014) and indicate that specific types of VGCs must be present in IS3 dendrites, and shape the bAP propagation and Ca 2ϩ signal generation.

Model: starting choices for VGC types and distributions
We considered several different channel types with different distributions. At first, the following four types of channel currents were considered: transient sodium current (I Na,t ), slow delayed rectifier potassium current (I Kdrs ), fast delayed rectifier potassium current (I Kdrf ), and A-type potassium current (I Ka ), since they are known to be present in hippocampal interneurons (Saraga et al., 2003;Lawrence et al., 2006;Hu et al., 2010) and represent a reasonable, minimal set of channel types that should be able to replicate the observed firing patterns in IS3 cells. The channel models used were obtained from a previous OLM cell model developed by Lawrence et al. (2006) without any alterations. The model equations are given in Table 4 along with ModelDB (Hines et al., 2004) and ICGenealogy (http://icg.neurotheory.ox.ac.uk/) reference numbers.
Intrinsic, subthreshold activities in the theta rhythm and irregular firing were captured in previous interneuron models using I Ka , persistent sodium current (I Na,p ), and somatically injected white noise current (Morin et al., 2010;Sritharan and Skinner, 2012). Since IS3 cells exhibited similar electrophysiological features (Fig. 1B, cell 1, blue trace), somatic I Na,p was also included. The I Na,p model was obtained from Uebachs et al. (2010), and the steady-state activation equation parameter values were altered to match the parameters used by Morin et al. (2010) and Sritharan and Skinner (2012). Note that the steady-state mathematical structures for the I Na,p model are identical in all three of these articles (Morin et al., 2010;Uebachs et al., 2010;Sritharan and Skinner, 2012). Collectively, these channels were the first set of channels used (i.e., S.1; Table 5).
We now have (1) experimentally obtained electrophysiological signature features to capture, (2) have a compartmentalized model morphology with appropriate passive properties (see Materials and Methods), and (3) made initial choices for VGCs and can proceed with the semi-automated strategy (Fig. 2).

Model: determining types, densities, and distributions of VGCs using the semi-automated strategy
We investigated a total of 12 different scenarios (not all are shown). This was done in multiples of three: VGCs in the soma only; distributed uniformly in the soma and dendrites; and distributed uniformly in the soma and dendrites with A-type potassium channels restricted to the soma. The first three scenarios contained our initial set of channels, while the second and third triads of scenarios varied the VGC types (see below). The fourth triad of scenarios possessed L-type and T-type calcium channels (see Discussion). As most of these scenarios were not able to encompass the determined electrophysiological features, they are not described any further.
We describe and present three of these scenarios in detail (Table 5), as follows: the initial set of channels in the soma only (S.1); the set of channels with faster potassium kinetics in the soma only (S.2); and the set of channels with faster potassium kinetics in the soma and dendrites (SD). Focusing our results on these three scenarios allows us to show the clear progression from our initial assumptions to our improved models, and also allows us to demonstrate the impact of dendritic channels. Note that distribution labels containing an "S" means somatic channels, "D" means dendritic channels, and ".x" denotes different versions of similar distributions (e.g., different channel types). We also investigated channels in specific dendritic subregions. These models possessed the label "SDprox", meaning that they had channels in the soma and proximal dendrites.
For S.1 (Fig. 4 A1,B1), in visualizing the parameter space using the CBDR plots, we found that a fairly large parameter space yields models that are capable of eliciting the key features of IS3 cells (i.e., colored pixels represent models that were not eliminated in step 2). However, looking more closely at the histograms of mea-surements from the database of model populations, we noticed that they could fall completely outside those of the experimental histogram measurements. For example, the spike threshold was too low, the half-width was too large (Fig. 4C, left histogram), and the spike afterhyperpolarization was too large. In addition to this, the depolarization height in the membrane potential during depolarizing current injections was also too small (i.e., not just in the ϩ50 pA step but also in the ϩ20 pA and ϩ500 pA steps). To quantify this "depolarization height" feature, we included a measure that calculates the difference between the average membrane potential before the current injection and the average membrane potential during the current injection (Table 2, third column).
To illustrate the advantage of this semi-automated strategy over a purely hand-tuned strategy, we show the conductance ranges at the start and end (steps 0 and 7, respectively) of the cyclic approach for the different conductances shown in Table 6. Note that they vary in their range, and in their maximum and minimum values. Determining this by hand tuning only would be tedious, time consuming, and potentially without success, as changing Ϫ͑Vϩ33͒ 12 ͬϪ1 -Note that the Na ϩ reversal potential (E Na ) and the K ϩ reversal potential (E K ) are set to values of 50 and Ϫ77 mV, respectively, in all simulations. Note also that ModelDB reference numbers and ICGenealogy reference numbers are indicated in the first column. Comparison of channel kinetic traces are available from ICGenealogy.
any particular conductance and still capturing salient features may require shifting other conductances into other ranges.
In the absence of any particular knowledge of VGCs in IS3 cells or of particular IS3 cell dynamics of functional importance, we constrained our models with a few chosen signature features and a quality metric (see Materials and Methods). This meant that a number of different current injection steps from experiment were lumped together for a given feature (e.g., for the depolarization block feature, current steps of 450 -700 pA were used, giving ranging histograms; Fig. 3Q,R). In essence, we used the experimental histogram measurements (Fig. 3) as a guide. We felt that it made more sense to simply compare each measurement histogram independently when making choices on channel types, distributions, and conductance ranges, as our goal is to obtain starting reference models that would suggest the types, densities, and distributions of VGCs. In this sense, we considered biological variability in a manual fashion (i.e., steps 6 and 7).

Subsequent investigations (step 7)
In noting that particular measurements could not be captured with the initial choice of VGCs, we considered changes in the set of VGCs being used (Fig. 4C, flowchart). In particular, we had observed through hand tuning that by using delayed rectifier potassium channels with faster kinetics, the spike half-width and afterhyperpolarization were decreased, leading to an improvement in the ability of the models to resemble the experimental signature features. This, in turn, allowed for VGC conductance densities that increase the depolarization height in the membrane potential during depolarizing currents. In terms Note that in all cases, each channel has a uniform distribution, whether it is restricted to the soma or distributed across the soma and dendrites. of using faster delayed rectifier potassium kinetics, one option is to remove the slow delayed rectifier channels from the model, and the other is to use faster time constants in the fast delayed rectifier channel model. With this in mind, we applied the semi-automated strategy, but without a slow delayed rectifier potassium channel and with a more generic fast delayed rectifier channel containing faster delayed rectifier potassium channel kinetics (channel model was taken from Saraga et al., 2003). The equations are given in Table 4. We show the scenario with only somatic channels (i.e., S.2) and with channels uniformly distributed across the soma and dendrites (i.e., SD; Table 5; Fig. 4A2,A3,B2,B3). For both channel distributions, measurement distributions generated from the models in each database showed more overlap with the experimental measurement distributions when compared with S.1 (Fig. 4C). For these channel distributions, we obtained an increase in the depolarization height of the membrane potential during depolarizing currents. However, it was not as large as that seen in the experimental recordings. Specifically, during the ϩ50 pA step in the experimental trace, we see a depolarization height of 29.42 mV (Fig. 1B, cell 2) and an experimental trace distribution ranging from 15 to 34 mV (Fig. 3G). While the S.1 top model fell outside this range, with a depolarization height of 10.83 mV (Fig. 4A1), both S.2 and SD top models fell closer to the experimental range with depolarization heights of 21.37 and 14.55 mV, respectively (Fig. 4A2,A3). Other measurements such as the spike half-widths and the spike afterhyperpolarizations were also improved. Thus, we were able to distinguish that the changes made in S.2 and SD from S.1 led to measured improvements in the models.
In terms of distinguishing between S.2 and SD, we looked at several additional factors. For example, we looked at the first spike time distributions in each database (Fig. 4C, right). Experimental recordings indicate that during regular spiking regimes, the first spike time of the current injection is usually early on (i.e., between 0 and 40 ms). Although S.2 ranges from spikes starting as early as 40 ms to ϳ120 ms, SD mainly ranges from 50 to 350 ms. This indicates that spike onset is often too late in all the models, but that S.2 performs better. Another factor that distinguished S.2 from SD is the presence of spike amplitude decreases during the regular spiking regime (i.e., at ϩ50 pA). Whereas the models generated in SD have spike amplitudes that increase during regular spiking, S.2 has the opposite: spike amplitudes that start off large and then decrease over the course of the current injection. This is observed in the top model for S.2 in Figure 4A2, compared with the experimental recording (Fig. 1B, cell 2, blue trace). As mentioned previously, decreasing amplitudes is a common feature observed in the experimental data, and it is therefore an important feature to replicate in our models.
The determined conductance value ranges are listed in Table 6 (see final ranges). We note the following: in SD, distributing the channels uniformly across the dendrites requires an overall decrease in conductance values (i.e., with the exception of persistent sodium, which is already quite small). For the most part, the conductance value ranges are similar among S.1, S.2, and SD, with the exception of the potassium channels, which generally have higher conductance values in S.2 and SD, likely due to the removal of the slow delayed rectifier channel.

Delayed rectifier potassium channel subunit analysis
To investigate the predictions generated by the models, we sought to examine the delayed rectifier potassium channel subunit composition in IS3 cells. Specifically, we wanted to know whether fast or slow delayed rectifier potassium channels existed in IS3 cells.
For slow delayed rectifier potassium channel subunit composition, it is generally known that Kv2.1 or Kv2.2 combine with Kv5.1, Kv6.4, or Kv9.1 to Kv9.3 (Lien et al., 2002;Coetzee et al., 1999). Fast delayed rectifier potassium channels, on the other hand, are generally known to be composed of Kv3.1 and/or Kv3.2 subunits, of which Kv3.1 composition yields faster time constants (Hernández-Pineda et al., 1999;Lien et al., 2002). To get an idea of expression of these subunits in the stratum radiatum (SR) of the hippocampus (i.e., the layer where IS3 cell bodies and the large majority of its dendrites are located), we examined the in situ hybridization data from mouse brain slices using the Allen Mouse Brain Atlas (Lein et al., 2007; for full documentation, see http://help.brain-map.org/display/mousebrain) as a point of reference. For the slow delayed rectifier subunits, it appeared as though there is very little expression, if any, in the SR or stratum lacunosum-moleculare (SLM). For fast delayed rectifier subunits, it appeared as though Kv3.1 subunit expression was much more prominent in the SR than Kv3.2 subunits.
Using immunohistochemical analysis of the subunit localization in soma and dendrites of IS3 cells, we were able to examine whether Kv2.1 and Kv3.1 were expressed. We found that Kv3.1 is present in both the soma and proximal dendrites of putative IS3 cells in stratum pyramidale (SP) and SR (Fig. 6A,B). However, we found no evidence for Kv2.1 expression in IS3 cells residing in SR (Fig. 6D), and its sparse expression in IS3 cell bodies located within SP that, based on the punctate pattern around the cell body (Fig. 6C), could be attributed to axonal boutons contacting IS3 cells (Fig. 6C,D). These data indicate that Kv3.1, but not Kv2.1, is expressed in IS3 cells. Thus, the model predictions that delayed rectifier potassium channels with faster kinetics are present and slow delayed rectifier potassium currents are absent in IS3 cells were confirmed.

Model adjustments to consider spike propagation
So far, it seems as though S.2 models present a reasonable scenario for IS3 cell channel types and distributions. On the other hand, if we consider S.2 as the most reasonable models, this would indicate that channels in IS3 cells are restricted to the soma and that dendritic channel distributions are not a suitable aspect to capture normal IS3 cell electrophysiological activity. However, given that our immunohistochemical analyses indicate that potassium channels are present in the dendrites, this is not appropriate. Further, our models so far only considered uniform channel distributions in the whole dendritic tree. Thus, we turn to our Ca 2ϩ imaging data and consider the spatial profile of the spike amplitudes, as derived from the bAP-CaTs profile, in two of our models: the highest-ranking model from S.2 and the highestranking model from SD. In particular, we used a methodology similar to that used experimentally by injecting somatic current (800 pA, 2 ms) to induce bAPs. We then qualitatively compared the bAP amplitude decay with distance from soma to the calcium signal amplitude decay seen in experimental recordings. As can be observed in Figure 7A, although S.2 bAP amplitudes decay, bAPs are still seen at large distances from the soma (i.e., small ϳ40 mV amplitude spikes are still observed at 200 m from the soma).
At the 150 m location, these amplitudes are 40% of the amplitude at soma, which is comparable to the average calcium signal decay of 50% at 150 m (Fig. 1E). However, there is a lack of calcium signal decay in the first 50 m of dendrite, which is not comparable to the S.2 models. SD on the other hand does not show any spike amplitude decay (Fig. 7B), regardless of the distance from soma, which is expected because of the high density of VGCs distributed uniformly across the dendrites. Altogether, these results suggest the presence of VGCs in at least the first 50 m of dendrites.

Channels only in soma and proximal dendrites: SDprox.1 and SDprox.2
To better capture the experimental data of calcium signal decay, we adjusted our models to have dendritic channels only in the proximal dendrites (i.e., up to 50 -100 m from the soma; see Materials and Methods) and used the semi-automated strategy to estimate the required VGC conductance densities to be able to capture the electrophysiological signature features. In addition to looking at variations in maximal conductance values, we considered channels at different distances from soma (i.e., 50, 70, or 90 m; Fig. 7C2,D2). We also considered model variants where A-type potassium channels were confined to the soma (SDprox.2) or uniform in the soma and proximal dendrites (SDprox.1). Models in the databases for SDprox.1 (Fig. 7C2) and SDprox.2 (Fig. 7D2) can have channels at 50, 70, or 90 m (Fig. 7, denoted as "VGC Distance" in the CBDR plots) on top of having different VGC conductance densities.
In Figure 7, C1 and D1, we see the simulated voltage traces from top models from the SDprox.1 database and the SDprox.2 database. See Table 6 for conductance ranges, and Table 7 for top model conductance values and VGC distances along dendrites. From the traces for both of these models, we see that they appear to be very similar in quality to the S.2 models, which shows that a scenario with VGCs in the proximal dendrites (preferably at 50 and 70 m) is capable of capturing the electrophysiological features. Notably, from the measurement histograms we see that these databases, similar to the S.2 database, fall within many of the experimental measurement distributions. Also, using the same recording procedure (i.e., an 800 pA, 2 ms duration somatic current injection; Fig. 7C3,D3), we found that bAP amplitude decay is attenuated in the first 50 m of dendrites in these two scenarios, unlike the continuous spike amplitude decay seen in the S.2 model. With these developed models in hand, we next determined the location and amount of synaptic input necessary for spiking to occur.

Exploration of required input for spike generation in IS3 cell models
It was previously shown that unitary synaptic conductance at IS3-OLM cell synapses is small, but at the same time, the synchronous recruitment of many IS3 cells through optogenetic stimulation could control the OLM cell firing at theta frequencies (Tyan et al., 2014). To understand how incoming inputs can integrate to trigger spiking in IS3 cells, we used our developed multicompartment models of IS3 cells. As described above, these were S.2 models with VGCs only in the soma, SD models with VGCs uniformly distributed in the dendrites, and SDprox.1 and SDprox.2 models that had channels in the proximal dendrites. Only the SDprox.1 and SDprox.2 models were appropriate for capturing the salient features of IS3 cells and spike attenuation along the dendrites as  . VGC distributions in proximal dendrites. A, Action potential amplitude deterioration along a dendrite section (tree 1) in the S.2 top model following a somatic 800 pA current injection for 2 ms (G Na,t soma ϭ 0.25 S/cm 2 ; G Na,p soma ϭ 0.0001 S/cm 2 ; G Ka soma ϭ 0.15 S/cm 2 ; G Kdrf soma ϭ 1 S/cm 2 ). B, Action potential amplitude deterioration along a dendrite section (tree 1) in the SD top model following a somatic 800 pA current injection for 2 ms (G Na,t soma/dendrites ϭ 0.06 S/cm 2 ; G Na,p soma ϭ 0.0002 S/cm 2 ; G Ka soma/dendrites ϭ 0.1 S/cm 2 ; G Kdrf soma/dendrites ϭ 0.1 S/cm 2 ). C1, SDprox.1 top model with channels that are uniform from the soma until the first 70 m of the dendrites using the Boltzmann function. Conductance values are as follows: G Na,t soma/dendrites ϭ 0.07 S/cm 2 ; G Na,p soma ϭ 0.000075 S/cm 2 ; G Kdrf soma/dendrites ϭ 0.25 S/cm 2 ; G Ka soma/dendrites ϭ 0.07 S/cm 2 . Note that the well as the immunohistochemical analyses, indicating the presence of fast delayed rectifier potassium channels in the proximal dendrites. To explore the impact of different spatial locations on the input attenuation along the dendritic tree of our models, we first performed an electrotonic analysis.

Electrotonic analysis of passive model
We used the Electrotonic Analysis tools in NEURON to compute the electrotonic distance throughout the entire dendritic arbor of the passive model (see Materials and Methods). Figure 8A shows the morphology with labeled trees, while the number of branch points, surface area, most distal length, and summed length from all branches are given in Table 8. Tree 1, for example, has less branching and longer sections, and has a larger surface area than the other dendritic trees, 2A and 2B.
In the axon initial segment and trees 1, 2A, and 2B, we found that electrotonic distance increases with distance from soma (Fig. 8B,C), which indicates signal attenuation throughout the model. We also observed differences be-tween trees in the magnitude of electrotonic distance. For example, an input into the dendrite of tree 1 versus tree 2A at 300 m is attenuated by an electrotonic distance of 0.9 versus 1.3 (Fig. 8B). Similarly, a somatic input attenuates by an electrotonic distance of 0.2 versus 0.4 when measured at 300 m along tree 1 versus tree 2A (Fig. 8C).
This electrotonic distance appears to be smaller in dendritic trees with less branching and larger surface areas (Table 8, tree 1). For example, while tree 2A showed more attenuation (i.e., larger electrotonic distances) than tree 1, we also noted that tree 1 had less branching (seven branching points) and a larger surface area (2463.79 m 2 ; determined by section diameters and lengths) than tree 2A (eight branching points and 2385.98 m 2 ; Table 8). For similar reasons, both tree 1 and tree 2A showed less attenuation than tree 2B (nine branching points and 1982.60 m 2 ; Fig. 8, Table 8). Although a greater level of attenuation was observed in the remaining axon segment, this is likely because of the adjusted passive parameters that were used to compensate for the removal of the continued current injection protocol is Ϫ100, ϩ20, ϩ 50, and ϩ500 pA. C2, Parameter spaces for SDprox.1, as visualized using CBDR. Note that 70 of 243 models did not get rejected. C3, Action potential amplitude deterioration along a dendrite section (tree 1) in the above SDprox.1 top model following a somatic 800 pA current injection for 2 ms. D1, SDprox.2 top model with channels that are uniform from the soma until the first 70 m of the dendrites using the Boltzmann function. Conductance values are as follows: G Na,t soma/dendrites ϭ 0.055 S/cm 2 ; G Na,p soma ϭ 0.00015 S/cm 2 ; G Kdrf soma/dendrites ϭ 0.295 S/cm 2 ; G Ka soma ϭ 0.07 S/cm 2 . Note that the current injection protocol is Ϫ100, ϩ20, ϩ 50, and ϩ500 pA. D2, Parameter spaces for SDprox.2, as visualized using CBDR. Note that 146 of 243 models did not get rejected. D3, Action potential amplitude deterioration along a dendrite section in the above SDprox.2 top model following a somatic 800 pA current injection for 2 ms. axon. Collectively, these results suggest that synapses located on dendritic trees with higher branching and smaller surface areas (e.g., trees 2A and 2B relative to tree 1) are less likely to elicit spikes at large distances from the soma because of a greater level of attenuation (Fig. 8).

Synaptic inputs
We simulated single presynaptic spikes that give rise to synaptic inputs at different locations along the dendritic arbors (see Materials and Methods) of the S.2 top model, the SD top model, the SDprox.1 top model, and the SDprox.2 top model. These simulations allowed predictions of the minimum synaptic weights (threshold weight) necessary to evoke spikes at the somata of IS3 cells (Fig. 9).
For the S.2, SDprox.1, and SDprox.2 models, there were distinct dendritic points where threshold weight rose sharply ( Fig. 9; i.e., the point marked by the bottom red dashed line at 0.5 S). These points correspond to the level of synaptic current at which the membrane potential approached the reversal potential (i.e., 0 mV). In other words, they correspond to when the increase in membrane potential following an excitatory synaptic event was ϳ70 mV (Fig. 9A1-A3). We identified this "saturation" by plotting the maximal change in membrane potential (i.e., maximum potential Ϫ minimum potential) at the synapse location within the first 1 ms of the presynaptic spike ( Fig.  9A1-A3, S.2 top model). When the change in membrane potential had a magnitude equal to the absolute value of the resting membrane potential (ϳ70 mV), the synapse was not generating any additional current, regardless of the magnitude of the weight. This occurred at synaptic weights beyond ϳ0.5 S, as indicated by the red dashed line.
Furthermore, we found that for each main dendritic tree (i.e., trees 1, 2A, and 2B), the point at which the threshold weight increased sharply (red dashed line) occurs at a different distance from the soma, as shown in the plot of Figure 9B for the S.2 model. In these plots, each vertical blue line marks a main dendritic tree (Fig. 8A, subset morphologies; Table 8, morphological analysis). Synaptic inputs to tree 1, with its smaller number of branching points and longer sections than other trees, generated somatic spikes up until 350 -400 m. Synaptic inputs to tree 2A could elicit somatic spikes up until 200 -250 m. Tree 2B had the largest amount of branching and the smallest surface area, and elicited spikes up until 150 -200 m. These observations make sense following our electrotonic analyses in which trees 2A and 2B had more attenuation relative to tree 1.
In Figure 9C1-C3, we show differences among S.2, SDprox.1, and SDprox.2, noting that SDprox.1 and Figure 8. Electrotonic analysis of the M2 morphology. A, A guide of the M2 morphology designating all of the morphological subsections, as well as showing both the spatial scale and diameter. B, Electrotonic distance along the dendritic arbor for voltage flowing into the soma. C, Electrotonic distance along the dendritic arbor for voltage flowing away from the soma. Note that the electrotonic distance is equal to the log value of attenuation, Figure 8. continued where attenuation is measured as voltage upstream/voltage downstream. More specifically, voltage upstream is an applied 1 mV signal, and voltage downstream is the downstream response to the 1 mV signal. In this sense, electrotonic distances Ͼ1 would imply a 10-fold attenuation in the signal.
SDprox.2 do not have passive dendrites like S.2. In SDprox.1, some of the dendritic trees show a slight decrease in the distance from soma at which the threshold weight increases past 0.5 S, suggesting that the dendritic channel distribution and densities (i.e., I Na,t , I Kdrf , and I Ka ) do not appreciably decrease the required input for spike generation when compared with a passive dendritic scenario (i.e., S.2). However, in SDprox.2 ( Fig. 9C1-C3), the location at which synaptic input could generate spikes was shifted forward by ϳ10 m in all of the branches, indicating that, by not including A-type potassium channels in the dendrites, spike generation is facilitated.
Finally, in the SD top model, we found that the required input (threshold weight) for spike generation was minimal (i.e., Ͻ Ͻ0.5 S) throughout the dendritic arbor and even decreased with distance from the soma (Fig. 9D, left y-axis). This is likely because synaptic saturation is never reached in the SD model (Fig. 9D, right y-axis). This suggests that uniform VGCs across the dendritic arbor could amplify distal inputs over proximal ones. However, since the threshold weight is fairly small across synaptic locations, it seems unlikely that any significant preference for distal inputs over proximal inputs would be observed.

Discussion
We have presented an efficient, semi-automated strategy to estimate the active VGC properties of identified single neurons that uses electrophysiological data and dendritic calcium imaging to constrain various possibilities. Using this approach, we generated databases of models that mimic the various features seen in IS3 cell experimental recordings. From these models, we obtained estimates of the VGC types, densities, and distributions that may be found in IS3 cells. Specifically, our models predict that fast potassium channel kinetics play a large role in generating appropriate electrophysiological activity, and their presence was confirmed with immunohistochemical analyses. Our models also predict that VGCs in the proximal dendrites best facilitate spike generation when A-type potassium channels are restricted to the soma. We note that, although the approach still involves a fair amount of hand tuning, the automated cycling brings forth a structured aspect that ensures that the included VGCs have conductance values that are balanced relative to each other so as to produce the observed electrophysiological features. In summary, this approach is a much more efficient approach (relative to hand tuning) that takes advantage of automated model databases and visualization tools, successfully predicting the presence of particular VGCs in IS3 cells.
Notably, in all cases where starting values were recorded, the starting hand-tuned parameters and ranges always differed from the final top model conductance values and ranges (Tables 6, 7). The particular VGC conductance differences could be small or large, but the improvement in model quality between parameter refinement cycles is apparent when looking at the CBDR plots between cycles (Fig. 5). Specifically, these plots highlight how the different channel-type choices and conductance values have clear impacts on model database quality. As such, the predicted relative conductance values are meaningful and are not simply a particular hand-tuned choice.

Proximal dendritic distributions of VGCs on IS3 cells
When analyzing the minimal synaptic weights necessary to elicit a spike along the dendritic arbor of our models, we found that these weights are minimized when the models possess uniform distributions of channels along the dendrites. Interestingly, for the models with VGCs in the proximal dendrites, the synaptic weights necessary to elicit a somatic spike are minimized when A-type potassium channels are restricted to the soma. This makes sense since the presence of additional outward potassium currents in the dendrite would effectively require larger, excitatory inputs to generate a spike. Furthermore, previous experimental (Magee et al., 1998;Magee, 2000;Cai et al., 2004;Yang et al., 2015) and modeling (Tigerholm et al., 2013) work shows that A-type potassium channels contribute to sublinear summation of dendritic spikes, thus dampening dendritic excitability during dendritic integration in pyramidal cells.
In terms of how well these different model types replicate IS3 cell features, we found that uniform distributions of channels along the dendrites (i.e., SD models) do not acceptably reproduce the electrophysiological features of IS3 cells, whereas passive dendritic cases (i.e., S.2 models) do. However, S.2 models have inappropriate spike amplitudes and propagation when considering dendritic Ca 2ϩ imaging data from IS3 cells. For these reasons, we suggest that the models with VGCs in proximal dendrites (i.e., SDprox) offer the most likely scenario for IS3 cell active properties, particularly when VGCs extend up to 70 m from the soma, since this is seen in the top models from the SDprox cases. This type of dendritic ion channel distribution would fall somewhere between the ion chan- ‫ء‬Note that to compute the surface area, the effective area was computed for each compartment via a trapezoidal integration across the compartment length using the NEURON area(x) function.

Figure 9.
Threshold weight with distance from soma in top models. A1-A3, Left y-axes, Synaptic weight threshold (in S) necessary to evoke a somatic spike in response to a single presynaptic spike applied incrementally at different points along the dendritic arbor for the S.2 top model. Right y-axes, S.2 top model local changes in membrane potential (i.e., maximum potential Ϫ minimum potential in the first 1 ms following the presynaptic spike) at the site of the synapse, where synaptic current saturation is reached once the local membrane potential reaches the reversal potential (i.e., an increase of ϳ70 mV). Note that tree 1 is plotted in A1, tree 2A is plotted nel distributions found in basket cells and OLM cells. Particularly, in previous basket cell modeling work, low bAP amplitudes in dendrites have been associated with a high potassium-to-sodium channel ratio (Hu et al., 2010). Conversely, in previous OLM cell models, the ratio of dendritic sodium channels to potassium channels was not low, leading to high bAP amplitudes across most of the dendrites (Martina et al., 2000;Saraga et al., 2003). Similar results for both basket and OLM cells were obtained in previous dendritic Ca 2ϩ imaging studies investigating bAP-CaTs amplitudes (Topolnik et al., 2009;Camiré and Topolnik, 2014). Since the bAP-CaTs observed for IS3 cells (Fig. 1E) fall somewhere in between what is seen in basket and OLM cells, it might follow that IS3 dendritic ion channel densities fall somewhere in between as well, as is indicated by our models. Additionally, it is known that in pyramidal cells, A-type potassium channel conductance values are much larger in distal apical dendrites than in proximal dendrites and soma (Magee et al., 1998;Golding et al., 2001). Although we did not investigate fully nonuniform dendritic ion channel distributions in our models, we explored two possible scenarios for A-type potassium channels: A-type potassium channel in soma and proximal dendrites (SDprox.1) versus that restricted to soma (SDprox.2). In these models, we explain a reduction in bAP amplitude by only having ion channels uniformly spread across a portion of the dendrites. In reality, this could be due to an increase in A-type potassium with distance from soma similar to what is seen in pyramidal cells (Golding et al., 2001). On the other hand, it could be similar to basket cells in that, past a certain distance from soma, the ratio of potassium to sodium channels becomes larger, resulting in a decay of bAP amplitude (Hu et al., 2010). Of course, the question of whether or not A-type potassium channels are present in IS3 dendrites remains to be determined.
Proximal dendritic VGC distributions on IS3 cells could serve several functions. For example, having large densities of VGCs in the proximal dendrites of IS3 cells may facilitate CA3 input gating, spike backpropagation, and induction of Hebbian forms of long-term potentiation at SC synapses (Golding et al., 2002;Topolnik, 2012). Accordingly, IS3 cells, through increased recruitment in response to CA3 input, may be responsible for OLM cell silencing during sharp wave-associated ripples (SWRs; Katona et al., 2014). In other words, VGCs in the proximal dendrites of IS3 cells might facilitate the transient SWR recruitment of IS3 cells. In terms of input onto the passive distal dendrites (i.e., in SLM) of IS3 cells in these models, we predict that either stronger input (e.g., simultaneous spatially distributed presynaptic spikes or high-frequency presynaptic input) or adjusted synaptic parameters (e.g., rise time or decay time) would be necessary to recruit IS3 cells to fire a somatic spike, since we found that input from a single presynaptic spike was inefficient in doing so.
This may be similar to CA1 pyramidal cells, which exhibit denser innervation through the perforant path relative to SC inputs (Kajiwara et al., 2008). In addition, in CA1 pyramidal cells, perforant path inputs evoked AMPAmediated EPSCs with longer rise and decay times than SC inputs (Otmakhova et al., 2002). Although the amplitudes of AMPA-mediated EPSCs from both inputs are not statistically different, the longer perforant path-evoked time course suggests a higher likelihood of synaptic summation. This also raises the possibility that the activation of distal inputs might trigger local cooperative mechanisms that would enhance the influence of SLM synapses on somatic spiking. Perforant path inputs on distal dendrites of CA1 pyramidal cells can initiate local Ca 2ϩ spikes through NMDAR-and voltage-gated calcium channelmediated Ca 2ϩ influx (Golding et al., 2002). It is conceivable that a functionally similar mechanism could exist in IS3 cells, albeit one that would likely involve a different combination of Ca 2ϩ sources, given the heterogeneity of local Ca 2ϩ mechanisms in CA1 interneurons (Camiré and Topolnik, 2012). Furthermore, the input from the entorhinal cortex generates theta rhythmic sinks coupled with theta rhythmic sources in the CA1 SLM and SLM/SR border (Kamondi et al., 1998). Since IS3 cells exhibit intrinsic theta oscillations and are capable of synchronizing OLM cells at theta frequencies (Tyan et al., 2014), they are well positioned to contribute to the generation of these SLM/SR rhythmic patterns.

Other model database approaches: similarities and differences
In comparison with fully automated approaches using other brute-force techniques (Prinz et al., 2003;Günay et al., 2008;Sekulić et al., 2014), evolutionary algorithm multi-objective optimization techniques (Druckmann et al., 2007(Druckmann et al., , 2013Hay et al., 2011) or control theoretical approaches (Brookings et al., 2014), the semi-automated continued in A2, and tree 2B is plotted in A3. The plots show that the location along the dendrites where there is an exponential increase in threshold weight (i.e., dashed red line at around 0.5 S) for each main dendritic tree approximately co-locates with the location along the dendrites where synaptic current saturation occurs. B, S.2 top model synaptic threshold weights (in S) in all dendritic trees. Plot shows that there are differences in the location along the dendrites where there is an exponential increase (i.e., dashed red line) in the threshold weight, depending on the dendritic tree of interest. C1-C3, Relative to passive dendrites (S.2 top model), active dendrites with A-type potassium channels in the first 70 m of dendrites (SDprox.1 top model) show a decrease in the distance from the soma at which the exponential increase in threshold weight is observed. On the other hand, active dendrites in the first 70 m of dendrites with A-type potassium channels restricted to the soma (SDprox.2 top model) shows an increase in the distance from the soma at which the exponential increase in threshold weight is observed. This is shown in all main trees (C1-C3). D, Note that the SD top model does not reach input saturation, and the amount of input to elicit a somatic spike does not increase exponentially in any of the dendritic trees. strategy has a somewhat different goal. Rather than a focus only on conductance values of the different VGCs in the model populations, it aims to also suggest what VGCs need to be present to match given electrophysiological features in cell types whose intrinsic properties lack any characterization. As such, we felt that retaining elements of hand tuning as well as lumping experimental data around chosen features (Fig. 3) were reasonable at this stage. Moving forward, our subsequently developed models with experimental confirmation of VGC types could then be used as reference models in fully automated approaches with more sophisticated algorithmic techniques.
Brookings et al. (2014) offer a control theoretical approach that focuses on temporal alignment with experimental data and is not feature based, unlike the approach used here. Hay et al. (2011) possessed experimental backpropagating action potential-evoked calcium spike recordings from proximal and distal dendrites in neocortical layer 5 pyramidal cells to constrain their model explorations. While this is feasible to obtain for pyramidal cells, dendritic recordings in many interneuron types are much more challenging, due to the smaller diameters of interneuron dendrites and the heterogeneity of interneuron types. For this reason, we relied on dendritic Ca 2ϩ imaging data, which can be indicative of action potential propagation based on previous works (CA1 pyramidal cells: Spruston et al., 1995;Golding et al., 2001;CA1 OLM cells: Martina et al., 2000;Topolnik et al., 2009;CA1 basket cells: Hu et al., 2010;Camiré and Topolnik, 2014).
As well, the approach used by Hay et al. (2011) relies on a base knowledge of layer 5 pyramidal cell channel types and conductance ranges, which has not yet been obtained for IS3 cells. Similarly, to generate a large database of OLM cell models necessary to explore conductance densities and channel distributions, Sekulić et al. (2014) based their models on experimental data for CA1 OLM cell channels and conductance ranges obtained from previous work before using a high-performance computing cluster. In this sense, our approach is useful for focusing on intrinsic properties that are unknown and are of potential theoretical and functional importance.
Although it is unclear whether or not the IS3 features chosen for the models to capture are sufficient, in previous work from Druckmann et al. (2011), it was found that simple current steps were appropriate stimuli to reveal the dynamics of the cell. Here we chose representative IS3 features from the experimental data that we wanted to capture, and particular current step sizes were chosen to encompass each of the four chosen features (Table 2). Our CIP choices were based on the experimental data and features in hand, and with these choices we aimed to obtain representative models. As we were not aiming for optimal models per se, we did not focus on examining different CIPs to see what might be ideal, but rather to make choices that were encompassing of the data. Also, our approach did not focus on automating the capturing of the variability in the experimental data (Druckmann et al., 2007;Brookings et al., 2014), but rather used non-automated, flexible hand tuning to ensure that model feature measurements fell within the distributions found experimentally.

Limitations
The major morphological limitation is whether or not axonal branches are important in this scenario (M1 vs M2 morphology). From one perspective, including this feature is computationally expensive and requires additional assumptions regarding axonal biophysical properties. From another perspective, excluding it can limit the uses of the model in future projects. Also it seems clear that a large proportion of IS3 surface area is occupied by axonal arborization (i.e., 65% in this case), which means that there is a fairly large loss in surface area once axonal branches are removed. Importantly though, once the passive parameters are optimized in both morphologies, the inclusion of somatic channels has similar effects (data not shown) on the action potential shape and timing measurements in both the M1 and M2 morphologies. Therefore, we can assume that adjusting the passive properties in the remaining axon segments is enough to counteract the effects of surface area loss without largely altering the spiking properties of the cell.
Since our models are minimalistic in regard to only having components that are needed to reasonably capture experimental output, several channels that are likely to be present in IS3 cells were not included. This is namely L-type calcium (I CaL ) and hyperpolarization-activated cyclic nucleotide-gated (HCN; I h ) channels. For one, I h seems to be present due to observations of channelspecific effects on the experimental voltage recordings (i.e., hyperpolarization sag). Second, our preliminary Ca 2ϩ imaging data show that I CaL makes a small contribution to bAP-CaTs in proximal dendrites. Additionally, a previous study (Vinet and Sík, 2006) has also indicated that calretinin-positive cells in the CA1 area of the hippocampus express L-type calcium subunits in small proportion, as well as T-type, N-type, R-type, and P-type calcium channels in larger proportions. Although we examined simulations with some of these channels using our approach as well as with only hand tuning, we did not observe any marked improvements in the measurements of the model, in comparison with the experimental data. In other words, model quality is usually maximized when conductance values for these channels are very small, regardless of the channel distribution or assortment of channel types. It is to be noted that these channel models were based on a previous OLM cell model (Lawrence et al., 2006). The involvement of these channels should, however, be investigated in larger database simulations (Sekulić et al., 2014) since they may likely play subtler roles in governing IS3 electrophysiology and dendritic Ca 2ϩ dynamics. For example, although the inclusion of a small HCN channel conductance in the model would drastically improve the hyperpolarization regime measurements, too high of a conductance would detrimentally affect all of the depolarization regime measurements. Finding an appropriate conductance parameter space for HCN channels might be possible using higher-resolution database searches with inclusion of calcium VGC types.
With more detailed parameter searches, it may also be possible to ameliorate the measurements in which the SDprox.1 and SDprox.2 databases had little overlap with the experimental data (e.g., spike rate, interspike interval, spike threshold).
Since our models include only minimal types of VGCs required, it might be premature to make direct comparisons with conductance values reported for other cell types. Nevertheless, when considering the estimates of channel conductance using a fast dynamic-clamp technique in CA1 oriens/alveus interneurons (Lien and Jonas, 2003), it was found that the addition of perisomatic Kv3 conductance [maximum conductance (G max ) ϭ 140 nS] with an increased deactivation rate could generate irregular firing patterns with dampened action potential amplitudes, similar to what is seen in IS3 cells. In our approach, however, we found that optimal conductance ranges for S.2 fast delayed rectifier potassium channels were between 0.25 and 1 S/cm 2 . Assuming a spherical soma with a radius of 10 m and a minimal conductance of 0.25 S/cm 2 , we obtain the following conductance estimate: G max ϭ (conductance) ϫ (Somatic Surace Area) ϭ (0.25 S/cm 2 ) ϫ (4r 2 ) ϭ (0.25 ϫ 10 6 S/cm 2 ) ϫ (4(0.001 cm) 2 ) ϭ S ϭ 3142 nS.
Although this is only an approximation of G max (i.e., at least 22 times larger than that in the study by Lien and Jonas, 2003), it points to a likely consequence of having a very minimal set of active properties (i.e., two types of inward currents and two types of outward currents in this case). We expect that with the addition of other channel types, the fast delayed rectifier potassium conductance values of the model will undergo a "balanced decrease" to be able to continue capturing IS3 cell features.
It is also worth mentioning that many assumptions regarding bAP amplitudes were derived from the bAP-CaT amplitudes, while the observation of bAP-CaTs implies the activation of calcium channels along with sodium and potassium channels, and, thus, may depend on the spatial distribution of Ca 2ϩ channels. However, for the purposes of this analysis, we have assumed that bAP-CaTs reflect the bAP spatial profile in at least the proximal dendrites. Notably, previous studies have shown that dendritic Ca 2ϩ imaging can be highly predictive of AP propagation in dendrites in other hippocampal cell types (CA1 pyramidal cells: Spruston et al., 1995;Golding et al., 2001;CA1 OLM Cells: Martina et al., 2000;Topolnik et al., 2009;CA1 basket cells: Hu et al., 2010;Camiré and Topolnik, 2014).
Finally, in our IS3 cell model, we incorporated a relatively simple model representation of stochastic gating (i.e., using Gaussian white noise), compared with more realistic models of stochastic gating (Fox, 1997;Dorval, 2006;Goldwyn et al., 2011). Despite this, since we know from previous work that Gaussian white noise, in combination with specific VGC types, is sufficient to elicit both irregular firing (Stiefel et al., 2013) as well as subthreshold spectral properties (Morin et al., 2010;Yoshida et al., 2011;Sritharan and Skinner, 2012), this minimal representation seems reasonable. Also, in the absence of particular biological details of stochastic gating in IS3 cells, a more detailed representation is not warranted.

Concluding remarks and future studies
As mentioned, it is possible to use our developed models as base reference models for larger-scale model database approaches to overcome some of the limitations when searching for appropriate parameter values (Sekulić et al., 2014) in the absence of detailed and larger sets of experimental recordings.
The developed models can also be used to consider functional contributions of IS3 cells by providing theoretical predictions on their recruitment and microcircuit interactions, keeping in mind the above limitations.
Although IS3 cells have been shown to exhibit both irregular and regular adaptive firing, the conditions required for either of these are unknown. It will therefore be necessary for experimentalists to investigate the functional ratios of excitatory and inhibitory inputs that control IS3 cell firing.
Furthermore, by modeling different layer-specific synaptic inputs to the developed IS3 model, we can predict what types of inputs can drive firing patterns similar to those observed in IS3 cells during electrophysiological recordings. Ultimately, we aim to use these multicompartment models to help understand the functional contribution of these cells to network oscillations such as theta rhythms.