A model of electrophysiological heterogeneity in periglomerular cells

Olfactory bulb (OB) periglomerular (PG) cells are heterogeneous with respect to several features, including morphology, connectivity, patterns of protein expression, and electrophysiological properties. However, these features rarely correlate with one another, suggesting that the differentiating properties of PG cells may arise from multiple independent adaptive variables rather than representing discrete cell classes. We use computational modeling to assess this hypothesis with respect to electrophysiological properties. Specifically, we show that the heterogeneous electrophysiological properties demonstrated in PG cell recordings can be explained solely by differences in the relative expression levels of ion channel species in the cell, without recourse to modifying channel kinetic properties themselves. This PG cell model can therefore be used as the basis for diverse cellular and network-level analyses of OB computations. Moreover, this simple basis for heterogeneity contributes to an emerging hypothesis that glomerular-layer interneurons may be better described as a single population expressing distributions of partially independent, potentially plastic properties, rather than as a set of discrete cell classes.


INTRODUCTION
Olfactory bulb (OB) periglomerular (PG) cells are heterogeneous with respect to several features, including morphology, synaptic connectivity, patterns of protein expression, transmitter identity, transcription factor expression, and electrophysiological properties Kosaka, 2005, 2012;Allen et al., 2007;Parrish-Aungst et al., 2007;Kiyokage et al., 2010). Interestingly, while subclasses of PG cells can be defined with respect to each of these features, the subclasses based on different features generally do not correlate-for example, knowing about a given marker expressed in a PG cell generally does not enable conclusions about its morphological or synaptic properties, its neurotransmitter, or its expression profile with respect to other molecular markers. Occasionally, however, correlations between specific features are observed. For example, while calbindin is expressed in a variety of OB juxtaglomerular interneurons including PG cells and superficial short-axon (sSA) cells (Kosaka and Kosaka, 2010), calbindin-positive PG cells in rats have been reported to receive no direct input from the olfactory nerve , indicating that calbindin is not expressed in the morphologically defined olfactory nerve-driven subtype of PG cell (PGo; Shao et al., 2009). As another example, only morphologically bipolar PG cells appear to express nicotinic cholinergic receptors in mice (Castillo et al., 1999). This pattern of mixed feature correlations suggests that PG cell heterogeneity may arise from multiple distinct and interacting processes rather than representing a discrete hierarchy of defined cell classes.
The electrophysiological properties of PG cells also are heterogeneous and defy easy categorization into discrete classes (Puopolo and Belluzzi, 1998;McQuiston and Katz, 2001;Shao et al., 2009). For example, PG cells have been grouped into subtypes based upon different patterns of potassium conductance expression (Puopolo and Belluzzi, 1998) as well as on whether or not they exhibit bursting properties based on a low-threshold calcium spike, as opposed to one of a range of non-bursting responses to somatic inputs (McQuiston and Katz, 2001). However, these physiological subtypes do not correlate with neuronal morphology, again reflecting the pattern of mixed feature correlations underlying PG cell heterogeneity. A further possibility is that the properties of any given cell may not be fixed, but plastic, enabling interneuronal response properties to be regulated in response to changing conditions or neuromodulatory state. Toward this end, we sought to build a biophysically-based reduced compartmental model of the PG cell in which the diverse electrophysiological response properties of PG cells that have been observed experimentally by McQuiston and Katz (2001) each could be evoked by varying the expression levels (relative conductances) of different ion channels, without modification of their kinetic properties or relative spatial distributions within the neuron. Specifically, McQuiston and Katz (2001) recorded from 21 presumptive PG neurons (in addition to recordings from other juxtaglomerular cells) and found them to exhibit diverse response profiles, from simple non-accommodating or accommodating spike trains to calcium T current-dependent low-threshold spikes (LTSs) crowned by single or bursts of sodium action potentials, along with some rarer responses. We here present a core PG cell model capable of replicating the diverse responses evoked under current clamp by McQuiston and Katz (2001) and suitable for extension pursuant to more specific simulation goals. Specifically, we replicate the various PG cell response profiles presented by these authors within a single computational model, modestly varying the conductance ratios of membrane channel classes in order to produce the diversity of responses exhibited by PG cells.

MODEL ARCHITECTURE
The PG cell model was implemented in the NEURON simulation environment (Hines and Carnevale, 1997) and was shaped by four principal objectives (Rubin and Cleland, 2006). First, it should represent cellular properties with an accuracy and precision adequate for its intended purposes. Second, it must be reduced in complexity for computational efficiency, such that it can be usefully incorporated into moderately-sized network simulations. Third, it should be reasonably compartmentalized, so that elements of the model can be upgraded or redesigned as required without rendering the revised model incompatible with its predecessors. Fourth, unconstrained variables should be minimized, as they can improve apparent fits to data while reducing models' predictive value. Consequently, as this model is intended as a baseline model for simulations incorporating PG cell heterogeneity, six sections were implemented: a soma, two equivalent dendrites, a gemmule (spine) shaft and gemmule body attached to one of these dendrites, and an axon stub (Figure 1). The separate gemmule sections are relevant to olfactory nerve-driven PGo cell models in which feed-forward olfactory nerve input onto the PGo gemmule results in local neurotransmitter release onto the dendrites of other OB neurons, a property thought to underlie odor similarity-dependent computations in the OB (Cleland and Sethupathy, 2006;Cleland, 2010). The axon compartment, appropriately extended, is relevant to the modeling of PG cell interglomerular projections, and potentially to the construction of an sSA neuron model. The multiple dendrites are relevant to modeling studies of dendrodendritic interactions across PG cells. In the present modeling study, we emphasized replication of current-clamp responses measured from PG cell somata in an OB slice preparation (McQuiston and Katz, 2001).

MODEL MORPHOLOGY AND MECHANISMS
The geometric parameters of the cell were derived from Pinching and Powell (1971). Specifically, the soma was 8 um in length and diameter; all other model sections were 1 um in diameter. The gemmule and shaft compartments were 1 um in all dimensions. The axon stub was modeled as 50 um long, and each dendrite as 20 um long. Passive membrane capacitance was 1.2 uF/cm 2 . Additional passive parameters are listed in Table 1. The modeling of cable properties for appropriate compartmentalization  Table 3); channel types marked as absent in a given section (-) were always absent. The nicotinic cholinergic receptor channel (G nic ), an ohmic cation channel, was expressed only in the gemmule compartment (see text for details). was governed by NEURON's d_lambda rule; accordingly, the axon comprised three compartments whereas all other sections were modeled as single isopotential compartments. Increasing compartmentalization by factors of 3 or 9 produced identical results. We utilized hyperpolarization-activated cation current (I H ) kinetics derived from PG cell recordings (Cadetti and Belluzzi, 2001). For other channel types, we adapted existing kinetic models as follows: a fast sodium spike current (I Na ) and delayed Frontiers in Computational Neuroscience www.frontiersin.org April 2013 | Volume 7 | Article 49 | 2 rectifier potassium current (I K ) (Destexhe et al., 1998), an A-type potassium current (Migliore et al., 1995), a low voltage-activated T-type calcium current [I Ca(T) ] (Destexhe et al., 1996), an Ltype high voltage-activated calcium current [I Ca(L) ] (Carlin et al., 2000;Markaki et al., 2005), a calcium-dependent non-specific cation current (I CAN ), a calcium-dependent potassium current [I K(Ca) ] (Destexhe et al., 1994), and a current noise source (David et al., 2008). The noise source constituted white noise (bandwidth 0-4 kHz, standard deviation 50 fA) colored by convolution with a single exponential function (τ = 5 ms). Perimembrane calcium concentration after calcium influx (within 0.1 um of the membrane) was modeled as decaying exponentially to a basal level via diffusion into the cell volume, using a mechanism based on that of Destexhe and colleagues (1998). Simulation kinetics assumed a temperature of 23 • C, consistent with the temperature at which the relevant slice recordings were performed (McQuiston and Katz, 2001). Membrane current kinetics were implemented in NMODL (Hines and Carnevale, 2000) and were identical across all model sections for each of the channel types implemented in the model ( Table 2). However, the levels of expression of each current mechanism (i.e., their maximum conductances) were varied to evoke different electrophysiological response properties in model simulations ( Table 3). The current I x through each channel was determined at every time step by the Hodgkin-Huxley formalism (Hodgkin and Huxley, 1952):
Frontiers in Computational Neuroscience www.frontiersin.org April 2013 | Volume 7 | Article 49 | 3 Table 3 | Active channel conductance densities (S/cm 2 ). Channel conductance densities were uniform across all model sections in which they were present at all (Figure 1B), except for the T-type low threshold calcium conductance, which was always expressed 5.667 times more densely in the dendritic and spine compartments than in the soma. Irregular spiking responses ( Figure 2D) also required the addition of a current noise mechanism (see Materials and Methods).

Figure PG response type
where, for each channel type x, G max defines the maximum conductance density of that channel in the designated compartment (Table 3), m and h are the activation and inactivation gating coefficients (of order a and b, respectively), V m is the membrane potential within that compartment, E x is the net reversal potential of the permeant ion(s), and is the temperature coefficient ( Table 2). Each gating coefficient was calculated from a differential function of its steady state (∞) and time constant (τ), both of which were functions of voltage. The gating coefficients of the I CAN and I K(Ca) currents also were functions of the internal perimembrane calcium concentration, which was governed by calcium influx and by the constitutive diffusion rate of calcium out of the 0.1 um deep perimembrane region.

COMPUTATIONAL METHODS AND PARAMETERIZATION
All simulations were run in NEURON version 7.2 (http://www. neuron.yale.edu) Carnevale, 1997, 2001;Carnevale and Hines, 2006). Timesteps were dynamically determined by NEURON's CVODE algorithm. The change in voltage in each compartment was calculated as the sum of all ionic currents, injected currents, and currents flowing from the neighboring compartments: where C m was the membrane capacitance and I x denotes particular channel-specific currents except as follows: I leak represents the leak current, I nic an ohmic cation current associated with the nicotinic cholinergic receptor channel, I adj denotes currents from adjacent connected model sections, and I inj denotes experimental current injection into the soma compartment. In the simulations presented here, membrane potential recordings were made from the soma in order to replicate experimental methods (McQuiston and Katz, 2001). Final parameterization was performed as described in Rubin and Cleland (2006). Parameters were additionally constrained so that the same parameter set produced both the depolarizationinduced and hyperpolarization-induced responses recorded experimentally from the same neuron, consistent with the corresponding experimental figures (McQuiston and Katz, 2001).

RESPONSE PROPERTIES COMPRISING SIMPLE SPIKE TRAINS
The expression levels of different membrane currents first were adjusted to replicate the experimentally observed range of simple (non-bursting) electrophysiological properties depicted in Figures 1A-D of McQuiston and Katz (2001). (The fifth simple response profile described in that work ( Figure 1E in McQuiston and Katz, 2001) was never observed in PG cells ( Table 2 in McQuiston and Katz, 2001) and hence was not modeled). Under the first parameter set (Table 3), depolarizing current evoked a non-accommodating series of action potentials (Figure 2A), whereas hyperpolarizing current evoked a modest H-current dependent sag and generated an action potential upon release (anode break response; Figure 2A). Higher levels of depolarizing activation evoked higher rates of non-accommodating action potentials, whereas stronger inhibition led to anode break responses evoking multiple action potentials. Other ratios of the same membrane conductances replicated the remaining simple responses of different PG neurons; specifically, trains of accommodating spikes upon depolarization with a prominent postinhibitory rebound after hyperpolarization (Figure 2B), single spikes upon depolarization without significant sag or rebound after hyperpolarization (Figure 2C), or irregular spiking patterns upon depolarization with a prominent sag and rebound-evoked spike burst after hyperpolarization ( Figure 2D). All responses were generated by altering the relative expression levels of five membrane mechanisms (Table 3), with two exceptions. First, evoking a non-accommodating train of spikes (Figure 2A) required a relatively depolarized leak reversal potential (−55 mV, rather than the −70 mV used for all other figures; Table 1). Alternatively, activation of a nicotinic receptor-induced current-modeled as a 5 mS/cm 2 ohmic cation conductance expressed in the single dendritic spine-changed accommodating and single-spiking neurons to instead produce non-accommodating spike train responses. In a network context, nicotinic neuromodulation of PG cells enhances their feedforward inhibition of mitral cells, contributing to the sharpening of odor representations within OB (Li and Cleland, 2013). Restriction of nicotinic receptor channels to the spine compartment reflected the distribution of cholinergic synaptic inputs to PG neurons observed via choline acetyltransferase immunohistochemistry (Kasa et al., 1995), though similar results were obtained using a uniform expression of nicotinic receptor channels. Second, the irregular spiking patterns of Figure 2D, replicating a relatively unusual form of activity among recorded neurons (McQuiston and Katz, 2001), were produced with aid of a current noise model (David et al., 2008; see Materials and Methods). Such activity may reflect irregular synaptic inputs from unusually active neurons, possibly owing to damage during slice preparation.

RESPONSE PROPERTIES INCLUDING CALCIUM T CURRENT-DEPENDENT LOW-THRESHOLD SPIKES
Additional distinct combinations of membrane conductance densities (Table 3) evoked a family of responses based on calcium T current-dependent LTSs and directly corresponding to those experimentally evoked from PG cells in Figures 2A-C of McQuiston and Katz (2001) and also observed by Zhou and colleagues (2006). Depolarizing current injected into one such model neuron produced an LTS crowned with a single fast sodium spike, whereas injection of hyperpolarizing current produced an anode break response upon release leading to an LTS with similar properties (Figure 3A). A different ratio of currents generated an LTS response upon depolarization that produced a burst of decrementing spikes, and a similar LTS response upon release from hyperpolarization ( Figure 3B). Both LTS responses resulted in substantially higher transient intracellular calcium concentrations in the model compared with non-LTS responses. Experimentally, these LTS-based responses were observed in PG cells twice as frequently as were the simple non-LTS spike trains (McQuiston and Katz, 2001). The repeated bursting observed in some JG neurons ( Figure 5E of McQuiston and Katz, 2001) is a property of external tufted cells (Hayar et al., 2004)-a discrete, glutamatergic class of JG neuron-and was not incorporated in the present PG cell model. The effects of diverse JG membrane currents on LTS properties have been analyzed in detail by Masurkar and Chen (2011a,b).

RESPONSE PROPERTIES INCORPORATING PERSISTENT PLATEAU POTENTIALS
Finally, the model was able to generate a persistent plateau potential upon both depolarization and release from hyperpolarization ( Figure 3C). This response was very rarely observed experimentally and therefore is poorly defined physiologically (McQuiston and Katz, 2001). In the model, this response profile required the inclusion of additional conductances beyond those that were capable of generating simple spiking responses and calcium LTSs (see next section; Table 3). The persistent plateau was not LTSdependent; reducing the T-type calcium current to zero did not affect plateau duration (Figure 3D), whereas modestly reducing the potassium delayed rectifier current ( Table 3) substantially extended the plateau (Figure 3E).

CONDUCTANCE TIMECOURSES
We mapped the timecourses of the most important membrane conductances in the model to examine their interplay. The singlespike LTS, whether evoked by direct depolarization (Figure 4Ai) or by an anode break response (Figure 4Aii), depended strongly on a low-threshold T-type calcium conductance. Its timecourse was primarily determined by that of the Ca (T) current itself, with a short post-LTS hyperpolarization punctuated by a modest deinactivation of the A-type potassium current. A roughly similar dynamical pattern of currents underlaid the LTS responses that included a decrementing burst of spikes, although the absolute dynamical levels of these conductances were reduced by roughly a factor of ten ( Figure 4B; note vertical scale). Specifically, whereas the maximal conductances of voltage-dependent membrane currents in PG neurons evoking this response ranged from 25 to 110% of those in PG neurons evoking single-spike LTSs (Table 3), the dynamical conductances were on the order of 10-20% of those in the latter group. The overall dynamical dependence on the interactions among these different membrane conductances was, however, similar for both forms of the LTS-based response.
In contrast to the simple and LTS-based responses described above, the persistent plateau potential responses replicated in Figures 3C-E required the inclusion of additional conductances. Lacking constraining pharmacology from their uncommon observation in source data, we modeled them using an established plateau potential combination of a voltage-dependent calcium current, a calcium-dependent cation current, and a calcium-dependent potassium current (Table 3). These currents are present in JG neurons, though they are not important determinants of the more common LTS timecourses (Masurkar and Chen, 2011a,b) and hence were omitted from simple and LTS simulations to reduce the dimensionality of the model (Table 3). Calcium entering the cell in response to direct depolarization (Figure 4Ci) or in response to H current-dependent depolarization following release from a hyperpolarizing pulse (Figure 4Cii) activated both the cation conductance G CAN and the potassium conductance G K(Ca) . The large cation conductance generated and maintained the plateau, declining along with the gradual decline in perimembrane calcium concentration until it could no longer be sustained. The calcium-dependent potassium conductance was primarily responsible for the termination of the burst of spikes following depolarizing input, but otherwise also followed the declining perimembrane calcium concentration. Unlike the much more common LTS responses, the role of the calcium T current in these persistent plateau potentials was negligible. Due to the absence of constraining pharmacological data regarding this response profile, it was not investigated further in the model.

PHARMACOLOGY OF LOW-THRESHOLD SPIKING RESPONSES
The ionic basis of LTS responses in PG cells was experimentally determined using pharmacological manipulations. Specifically, blockade of sodium influx by choline replacement of external sodium or by intracellular infusion of the fast Na channel antagonist QX314 did not eliminate the LTS from PG cells (McQuiston and Katz, 2001). Similarly, in the model (Figure 5), reduction of the fast sodium conductance to zero prevented action potential generation but preserved the LTS response to both depolarizing current injection and release from hyperpolarization ( Figure 5B).
In contrast, reducing the T-type calcium current to zero (as if superfusing the tissue with nickel ions) completely blocked the LTS responses ( Figure 5C), as demonstrated experimentally (McQuiston and Katz, 2001). The ionic basis for LTS responses in the model thus reflects those ascertained experimentally in PG neurons.

DISCUSSION
Network models incorporating specific biophysical elements of interest can produce useful predictions about the impact of these elements even if the implementation of other cellular properties, such as the cable properties of extensive neuritic arbors, is generic. Moreover, for purposes of understanding network properties, computational efficiency is critical, and it has been clearly established that morphologically reduced (oligocompartmental) models are capable of replicating the global response patterns of complex neurons, even when single-compartment models fall short (Davison et al., 2000;Rubin and Cleland, 2006;Li and Cleland, 2013). Within these constraints, the present model replicates the diverse physiological and pharmacological responses of PG neurons as observed experimentally. Pursuant to particular applications, the model architecture can be morphologically elaborated or further reduced with manageable effects on whole-cell performance.

INTERNEURONAL DIVERSITY IN OLFACTORY BULB: DISTRIBUTIONS OR DISCRETE CLASSES?
We demonstrate computationally that each of the diverse simple and LTS-based responses observed experimentally in PG cells by McQuiston and Katz (2001) can be generated by relatively modest changes in the relative levels of expression of five membrane conductances without altering their kinetic properties or relative spatial distributions across the cell (Figures 2A-C, Figures 3A,B; Table 3; this summary excludes the noise model of Figure 2D and the rare persistent plateau potential responses of Figures 3C-E.).
While the prevalence of each of these "common" physiological response profiles in vivo is not known, this simple basis for heterogeneity contributes to the emerging hypothesis that glomerularlayer GABAergic and dopaminergic interneurons may be better described as a single population expressing distributions of partially independent properties, rather than as a series of discrete cell classes . PG cells have been subdivided into categories on the basis of morphology (bipolar vs. monopolar, axonless or with axon), extending dendrites to either contact the axons of primary olfactory sensory neurons (PGo, or "type 1") or not (PGe, or "type 2"), marker expression (e.g., calbindin, calretinin, parvalbumin, nitric oxide synthase, nicotinic cholinergic receptor), and presumptive neurotransmitter (GABA, dopamine, both) (Kosaka and Kosaka, 2005Shao et al., 2009). While there are occasional, possibly important correlations between certain feature sets (e.g., see Kosaka et al., 2001), generally these properties do not predict one another, appearing to be partially independent and forming a complex pattern of mixed feature correlations. Indeed, even sSA neurons, while exhibiting broad axonal projection patterns quite dissimilar from those of classically described PG cells (Aungst et al., 2003;Cleland, 2010), may not constitute a genuinely separate class of neurons. Like PG cells, sSA neurons are thought to utilize GABA and/or dopamine as neurotransmitters [though they were initially identified as glutamatergic; (Aungst et al., 2003;Brill et al., 2009) relatively proximate glomeruli, not unlike that nominal subclass of PG neurons that projects axons to one or a few neighboring glomeruli. The sSA cells that project long-distance axons to distant glomeruli are relatively uncommon (Aungst et al., 2003)albeit essential to the small-world network computation proposed to mediate global feedback normalization (Cleland et al., 2007;Cleland, 2010)-suggesting that these sSA cells may constitute one tail of a distribution of morphological properties exhibited by a single heterogeneous population of interneurons, with the axonless form of "PG" interneurons constituting an opposing tail. Indeed, very recent work (Kiyokage et al., 2010;Marbach and Albeanu, 2011;McGann, in press) has begun to employ a revised nomenclature relabeling some traditionally-construed PG cells as sSA cells. Specifically, GAD67-expressing juxtaglomerular cells-which often coexpress tyrosine hydroxylase and include traditionally recognized sSA cells as well as many traditionallyconstrued PG cells-are being labeled as sSA (or simply SA) neurons in this nomenclature whereas GAD65-expressing PG cells retain the PG label. Both groups are heterogeneous with respect to other features, though GAD65-expressing PG cells do not exhibit the most broadly projecting of the multiglomerular morphologies. Whatever its merits, this revised nomenclature highlights the arbitrariness of discrete classification schemes for PG/sSA neurons. Moreover, with continuing study, increasing numbers of OB interneurons have been observed that violate the defining features of any of these nominal classes. For example, some nominal sSA neurons feature extrabulbar axonal projections, and some DA/GABAergic juxtaglomerular neurons project to the contralateral OB Kosaka, 2008, 2011), a profile generally associated with the glutamatergic external tufted cells (Schoenfeld et al., 1985). If this single heterogeneous interneuron population hypothesis is correct, then the fundamental experimental question shifts from "what are each of these cell types for?" to "what factors cause or require the expression of particular markers or properties within neurons of this [single] population?" That is, the observed pattern of mixed feature correlations suggests that PG cell heterogeneity may arise from multiple distinct and interacting processes rather than representing a discrete hierarchy of defined cell classes, and that certain nominal class markers may signal relatively transient neuronal states rather than durable identities. For example, while the calcium binding protein markers expressed heterogeneously among PG neurons have not typically been associated with particular functions, recent work has implicated calbindin-D28k, particularly, in some forms of synaptic plasticity (Schmidt, 2012;Westerink et al., 2012). Consequently, the presence of this marker in a subset of PG cells might identify those PG cells that have been recently synaptically modified. More broadly, the properties of any given PG cell may not be fixed, but plastic and separately regulated, thereby enabling interneuronal response properties to adapt in response to changing conditions or neuromodulatory state. The concept of discrete, static classes of OB interneuron may no longer serve.