|Home | About | Journals | Submit | Contact Us | Français|
Conductance-based neuron models aid in understanding the role intrinsic and synaptic currents play in producing neuronal activity. Incorporating morphological detail into a model allows for additional analysis of non-homogeneous distributions of active and synaptic conductances, as well as spatial segregation of electrical events. We developed a morphologically detailed “Full Model” of a leech heart interneuron that replicates reasonably well intracellular recordings from these interneurons. However, it comprises hundreds of compartments, each increasing parameter space and simulation time. To reduce the number of compartments of the Full Model, while preserving conductance densities and distributions, its compartments were grouped into functional groups that each share identical conductance densities. Each functional group was sequentially reduced to one or two compartments, preserving surface area, conductance densities and its contribution to input resistance. As a result, the input resistance and membrane time constant were preserved. The axial resistances of several compartments were rescaled to match the amplitude of synaptic currents and low-threshold calcium currents and the shape of action potentials to those in the Full Model. This reduced model, with intrinsic conductances, matched the activity of the Full Model for a variety of simulated current-clamp and voltage-clamp data. Because surface area and conductance distribution of the functional groups of the Full Model were maintained, parameter changes introduced into the reduced model can be directly translated to the Full Model. Thus, our computationally efficient reduced morphology model can be used as a tool for exploring the parameter space of the Full Model and in network simulations.
Conductance-based models of neurons demonstrate how membrane conductances shape the electrical activity of a neuron, and how alterations in membrane conductances, such as by modulators or transmitters, can alter the activity of a neuron (for examples, see: DeSchutter and Bower 1994, Goldman et al. 2001, Golowasch et al. 1992, Traub et al. 1991). Additionally, conductance-based neuron models are used to explore how neurons contribute to the function of neural networks (for examples, see: Prinz et al. 2004, Zhang et al. 2003). In our work, the development of a conductance-based, single-compartment model of an oscillator interneuron of the leech heartbeat pattern generator has forged an understanding of how membrane conductances generate the bursting activity of these interneurons (Cymbalyuk et al. 2002; Hill et al. 2001; Nadim et al. 1995; Olsen et al. 1995). However, single-compartment models do not take into account that living neurons have spatial extent such that different regions of the membrane experience localized membrane potentials and non-homogeneous distributions of active and synaptic conductances. Adding morphological details to a model introduces non-uniform voltage changes and non-uniform distributions of active and synaptic conductances allowing us to explore the spatial segregation of electrical events such as synaptic potentials, plateau potentials and action potentials, and the propagation and attenuation of these signals. Building upon the insights we have gained from the single-compartment model, our goal was to create a morphology-based model of an oscillator heart interneuron to use for exploring effects of non-homogeneous distributions of active and synaptic conductances and for simulations of the heartbeat pattern generator network. Specifically, we wanted a model with detailed morphology with which we could study the role that conductance distributions play in producing the spiking and bursting activity of a heart interneuron in the heartbeat pattern generator.
The oscillator heart interneurons are spatially extended (Fig. 1A). From the soma branches the primary neurite, a focus of electrical activity and synaptic integration, from which issue many fine processes (secondary neurites) that participate in generating local electrical signals and serve as sites for both synaptic input and output. The primary neurite tapers into an axon that travels to the posterior adjacent ganglion. Our experimental data indicate that active conductance densities and synaptic inputs in the heart interneurons are non-homogeneous and give rise to spatially segregated electrical events. Photo-ablation studies have shown that synaptic inputs (and thus synaptic conductances) seem confined to the distal regions of the secondary neurites and action potentials are elicited only in the proximal portion of the axon (Ivanov and Calabrese 1998), indicating a higher fast voltage-gated Na conductance density. Not only are the electrical events within a heart interneuron affected by the spatial segregation of ionic currents and synapses but also so are network functions. For example, multicompartment models of coordinating fibers in the heartbeat timing network indicate that the extended morphology of the fibers is important for coordinating the activity of the network (Jezzini et al. 2004).
We created a Full Model of a single leech heart interneuron, incorporating the detailed morphology and likely active and synaptic conductance distributions of the interneurons. We tuned the conductance densities to match intracellular recordings from synaptically isolated interneurons. While the benefits of including spatial extent into a neuron model are obvious, each additional compartment increases the parameters that need to be tuned and the computation time of the simulation. A neuron model with detailed morphology such as our Full Model requires hundreds of compartments in order to avoid discretization errors associated with approximating a spatially continuous neuron by a discrete model. Consequently, a single model simulation becomes frustratingly slow, and a network simulation with tens or more of such neurons becomes nearly unfeasible (although not impossible with heavy computing resources, see: Traub et al. 1999, Traub et al. 2005). These disadvantages created a need for a simplified model with relevant morphological information, but far fewer compartments that would allow systematic parameter variation/searches and could be used in network simulations.
Most simplified-morphology models fall into one of three categories: (1) the compartment geometries and thus membrane-area specific properties are somewhat arbitrarily chosen (such as Traub et al. 1991, Pinsky and Rinzel 1994, and Davison et al. 2000), (2) the compartment geometries are chosen to preserve the axial resistance of a detailed morphology model (Bush and Sejnowski 1993), or (3) compartment geometries are chosen to preserve the surface area of a detailed morphology model, such as the ‘equivalent cylinder’ method developed by Rall (Burke 2000; Rall 1959; Rall 1964) and expanded upon by others (Agmon-Snir 1995; Burke 2000; Clements and Redman 1989; Zador et al. 1995). We wished to maintain similarity of the intrinsic and synaptic conductance distributions between the detailed and simplified models, and therefore we chose a method which preserved surface area. However, we sought a greater simplification than the equivalent cable techniques that reduce a branching structure into an unbranched cable. We developed our own method, similar, in part, to one developed by Destexhe (2001), to reduce the morphological details of the heart interneurons. We divided the neuron into functional regions corresponding to distinct anatomical and physiological role. We reduced all compartments within a functional region to one or two compartments, with a total reduction of ~900 to 9 compartments. While the surface area of each functional region was preserved, we set the axial resistance of each reduced compartment to preserve that functional region's contribution to the total neuron input resistance. Next we added active conductances and adjusted the axial resistances of several compartments to improve the similarities of the reduced model to the Full Model.
The result of our reduction method is a reduced model (the ‘Actively Reduced Model’) that reproduces not only the passive, but also active properties of the detailed morphology Full Model. The Actively Reduced Model, with only 9 compartments, runs over 160 times faster than the Full Model, with 871 compartments, while preserving the non-homogeneous active and synaptic conductance distribution as well as spatially segregated electrical events. The Actively Reduced Model, with its direct morphological mapping to the Full Model and its decreased simulation time, is an ideal tool for more expedient tuning and parameter exploration that can then be directly applied to the detailed morphology Full Model.
Leeches (Hirudo medicinalis) were supplied by Leeches USA (Westbury, NY) and Biopharm (distr. by Carolina Biological Supply Co., Burlington, NC) and maintained in artificial pond water at 15°C. Before dissection, animals were cold anesthetized (1-4°C). Ganglia from midbody segments 3 and 4 were removed, individually pinned in Petri dishes lined with silicone elastomer (SYLGARD, Dow Corning, Midland, MI), and covered with leech saline. The saline contained, in mM: 115 NaCl, 4 KCl, 1.8 CaCl2, 10 glucose, 10 HEPES buffer; pH is adjusted to 7.4 with NaOH. Immediately prior to the experiment, the glial sheath was removed with microscissors or scalpel, then the preparation was superfused continuously with saline (bath volume ~0.5 ml, flow rate ~2 ml/min). All experiments were performed at room temperature.
Borosilicate microelectrodes (1 mm o.d., 0.75 mm i.d.) were pulled to sharp tips, and filled with 4 M KAC, 20mM KCl solution, resulting in electrode resistances of 15 to 30 MΩ. To decrease capacitance, some electrodes were coated with silicone elastomer to within ~500 μm of the tip; all intracellular electrode tips were dipped in dimethyl polysiloxane. Heart interneurons were identified by their position in the ganglion and by their characteristic bursting pattern (Fig. 2). Neurons were impaled in leech saline, and only cells with input resistances greater than 60 MΩ, as measured by a -0.1 nA pulse, were accepted for experiments. At the end of each experiment, the bath potential was recorded, and only those experiments where it was within ± 5 mV of zero were accepted.
Current-clamp and voltage-clamp experiments were performed using an Axoclamp-2A or Axoclamp-2B amplifier (Axon Instruments, Union City, CA) operating in discontinuous current-clamp or discontinuous single electrode voltage-clamp mode with a sample rate of 2.5 – 2.8 kHz. The electrode potential was monitored to ensure that it settled during each sample cycle. Output bandwidth was 0.3 kHz. Voltage-clamp gain was 8 to 20 nA/mV. Clampex software (pClamp 9.0, Axon Instruments) generated current or voltage waveforms during current-clamp and voltage-clamp recordings respectively, digitized, and stored all electrophysiology recordings. Time constant fitting was performed with Clampfit (pClamp 9.0, Axon Instruments).
In voltage-clamp, during hyperpolarizing voltage ramp protocols, saline included 2 mM CsCl to block the hyperpolarization-activated cation current, Ih, and 0.5 mM bicuculline methiodide (Sigma, St. Louis, MO) to block synaptic current (Schmidt and Calabrese 1992). Membrane potential was clamped at -40 mV and linearly ramped to -85 mV in one second, then linearly ramped to -40 mV in one second. The voltage-clamp current corresponding to the descending and ascending voltage ramps was plotted versus voltage and the current between -70 and -80 mV, visually assessed to be the most linear portion of the current, was linearly fit. Input resistance, Rin, was calculated as the inverse of the slope of the fit; leak reversal potential, ELeak, was estimated as the y-intercept of the fit (zero current). Poor electrode penetration can yield misleadingly low Rin and less negative ELeak values, skewing calculated averages. Therefore, we did not use average values but rather chose a representative value from where measurements clustered.
Plotting the voltage-clamp current versus voltage and measuring the current offset between the hyperpolarizing (charging) and repolarizing (discharging) voltage ramps gave twice the capacitative current. Input capacitance was calculated by dividing capacitative current by the change in ramp voltage over time. Input capacitance was also measured by applying -0.1 nA current pulses and fitting the membrane time constant. Instead of calculating an average value, we chose a representative value from where measurements clustered; both protocols yielded similar results.
Spontaneous firing and spike frequency response to injected current was measured by synaptically isolating heart interneurons in saline containing 1 mM bicuculline (Cymbalyuk et al. 2002) and injecting current (sequential 12 s steps, in nA: -0.05, 0.05, -0.1, 0.1, -0.15, 0.15, -0.2, 0.2, with 12 s of 0 nA current between each). The spike count in the last 5 s of each step was divided by 5 s to calculate spike frequency. During the current pulse protocol, spontaneous spikes occurring during no current injection were selected as examples of spike shapes from isolated living heart interneurons.
The morphology of a heart interneuron was visualized by filling the cell with Neurobiotin, which was then developed with avidin, similar to methods described by Gilchrist et al. (1995) and Jaeger (2001). The geometry of the Full Model was reconstructed using Neurolucida software (MicroBrightfield, Inc.) (Fig. 1A). The traced cell was translated into isopotential compartments (871 total) loosely following an algorithm by Eichler West et al. (1997). In brief, the total electrotonic length was computed for each branch of neurite, and each branch was then divided into compartments such that the diameter of each compartment corresponded to the segment of neurite it represented, and the length of each compartment was less than 0.1 λ. This meshing procedure was accomplished using a version of the CVAPP utility by Robert Cannon (University of Edinburgh) that was modified in-house to mesh and export to the GENESIS neural simulation software (Bower and Beeman 1998).
The membrane potential Vm in each isopotential compartment is determined by:
where cm is the membrane capacitance (scaled by surface area), Iext is the externally injected current, Iaxial is the current from adjacent compartments, and n enumerates the membrane currents in the compartment. Simulations were performed using GENESIS neural simulation software (Bower and Beeman 1998) run with a Linux operating system on a 2.8 GHz Pentium III PC. The implicit Crank-Nicholson integration method was used with a time step of 0.1 ms.
We grouped the compartments into functional groups based on their anatomical position and physiological function. All compartments within a functional group were assigned identical intrinsic and synaptic conductance densities. These functional groups were: Soma, Neurite-soma (Neurite-to-Soma transition), Neurite, Neurite-axon1 and Neurite-axon2 (Neurite-to-Axon transition), Axon, Secondary Neurite, and Synaptic Compartments. Two Neurite-soma and eleven Neurite compartments constitute the primary neurite process branching directly from the single spherical Soma compartment (Fig. 1A). Distal to the primary neurite loop, the neurite process narrows into the Neurite-axon1, Neurite-axon2 and Axon regions, which were modeled with 10, 4, and 132 compartments, respectively. The secondary neurites, smaller processes branching from the primary neurite, were modeled with 663 Secondary Neurite compartments. Synaptic contacts were assigned based on electron microscopical observations (Tolbert and Calabrese 1985) as most likely occurring at the most distal portions of the secondary neurite and were represented by 48 Synaptic Compartments. The geometries of the compartments are given in Table 1. The conductance densities of the compartments are given in Table 2.
Axial resistivity of 250 Ωcm was chosen from literature on leech neurons (Fromherz and Müller 1994). Values for specific membrane resistance, Rm, leak reversal potential, ELeak, and specific membrane capacitance, Cm, were derived from current pulses and hyperpolarizing voltage ramps in heart interneurons (see Methods: Electrophysiology). In the model, Rm was hand fit to match the input resistance, Rin, measured as the amplitude of voltage change during a simulated current pulse; Cm was set to match the time constant of the voltage change.
Eight voltage-gated conductances have been identified in heart interneurons, all of which were implemented in the model. Simon et al. (1992) described three outward currents: a fast transient K current (IKA), and two delayed-rectifier-like K currents: one inactivating (IK1) and one persistent (IK2). Angstadt and Calabrese (1991, 1989) characterized a hyperpolarization-activated inward cation current (Ih) and two low-threshold Ca currents: one rapidly inactivating (ICaF) and one slowly inactivating (ICaS). Two Na currents have been identified: a low-threshold persistent Na current (IP) (Opdyke and Calabrese 1994) and a fast Na current involved in spiking. The inhibitory synaptic transmission between heart interneurons has both spike-mediated (Simon et al. 1994) and graded components (Angstadt and Calabrese 1991); only the spike-mediated component was modeled here. These currents are modeled as Hodgkin-Huxley type conductances (Hodgkin and Huxley 1952). The kinetics, voltage-dependencies and reversal potentials of these conductances were implemented in the model as published in Hill et al. (2001), except the voltage dependence of activation and inactivation, where applicable, of IK1, IK2, and IKA are shifted by +10 mV, and those of INa are shifted by -2 mV.
Conductance distributions were assigned to approximate physiological data and to limit parameter space for tuning the model. We restrict the Soma to only K conductances, the Axon compartments to only fast Na and K conductances, and the Secondary Neurite compartments to K and Ca conductances. The Synaptic Compartments' densities are fixed to be identical to those in the Secondary Neurites, except that they uniquely also contain synaptic conductance. The transitional compartments, Neurite-soma, Neurite-axon1, and Neurite-axon2, were designed to smooth the transition between compartments of abutting functional regions that have different conductance densities. Where we use one transition region, the Neurite-soma, the conductance densities are the average of those in the Soma and Neurite. Where we use 2 transitional regions, the Neurite-axon1 and Neurite-axon2, we use a 1/3, 2/3 graded average of the neighboring functional regions, such that (except for the fast Na conductance), the Neurite-axon1 conductance densities are the sum of 2/3 of those in the Neurite and 1/3 of those in the Axon; the Neurite-axon2 densities are 1/3 of those in the Neurite and 2/3 of those in the Axon. To set the spike initiation zone, fast Na conductance density was set as 1/2 of Axon density in the Neurite, 3/4 in Neurite-axon1, and 1 in Neurite-axon2 to insure action potentials initiate in the Neurite-axon2 compartments. Experimentally, spontaneous action potentials in heart interneurons are not observed when the area near the proximal axon/distal neurite is removed by photoablation, thus indicating the region near the Neurite-axon2 compartments is likely the action potential initiation zone in heart interneurons, as it is in the model (Ivanov and Calabrese 1998).
First we describe the construction and parameterization of a Full Model of a heart interneuron. Then we describe the construction of a Reduced Model by an iterative process and the evaluation of the similarity of the two models by a number of fit metrics.
Based on a Neurolucida image of a dye-filled heart interneuron, we created a model capturing its highly detailed morphology (Fig. 1A). To set passive parameters, we chose the specific axial resistance, Ra = 250 Ωcm, from the literature (Fromherz and Müller 1994), and set the specific membrane resistance, Rm = 1.6×104 Ωcm2, leak reversal potential, ELeak = -50 mV, and specific membrane capacitance, Cm = 2.2 μF/cm2, to match empirical data, namely a passive input resistance Rin = 105 MΩ, a leak reversal potential ELeak = -50 mV, and an input capacitance Cin = 360 pF, derived from current pulses and hyperpolarizing voltage ramps in heart interneurons (see Methods). As this model geometry captures the full branching pattern and spatial extent of the interneuron, we refer to it as the ‘Full Model’ to distinguish it from the reduced morphology models (Figs. 1B and C) we developed and will discuss in this text.
In addition to leak current (ILeak), which is accommodated in the present models by Rm and ELeak, eight voltage-gated conductances, as described in Methods, were included in the model. As mentioned above, there is likely a non-homogeneous distribution of membrane conductances in heart interneurons. Initiation of spikes at the distal neurite/proximal axon suggests a higher concentration of fast Na current in the neurite near the axon than in the regions near the soma. Synaptic conductances appear to be distributed at the branched distal ends of the secondary neurites (Tolbert and Calabrese 1985), and furthermore, calcium imaging studies suggest low-threshold Ca currents are denser in these distal regions (Ivanov and Calabrese 2000, 2003). These spatially distinct properties enabled us to postulate anatomically and functionally distinct regions of the neuron: the soma, primary neurite, secondary neurites, synaptic contacts, and the axon. To organize conductance distribution in the model neuron, and to reduce parameter space, we grouped all compartments in the model into functional or transitional groups, where all compartments within each functional/transitional group were assigned the same conductance densities. The functional groups correspond to the identified functional regions of the neuron: Soma, Neurite (corresponding to the middle part of the primary neurite), Secondary Neurite (corresponding to the secondary neurites), Synaptic Compartments and Axon. To roughly approximate linear changes in conductance distributions between functional regions, we also specified transitional groups: Neurite-soma (connecting the Soma to the Neurite), and the Neurite-axon1 and Neurite-axon2 (connecting the Neurite with the Axon) (Fig. 1), such that the conductance densities of each transitional group were assigned intermediate values of the densities of abutting functional groups. All the compartments within a functional/transitional group were assigned the same conductance densities.
We first set conductance densities so that simulated voltage-clamp protocols produced conductance amplitudes that matched conductance amplitudes from published voltage-clamp data, amplitudes from simulated voltage-clamp in the single-compartment model, or values that represented a compromise between the two (Table 3). As there are no voltage-clamp data for INa, we initially set the Axon conductance density of Na to that of the single-compartmental model of a leech heart interneuron (Hill et al. 2001). We then adjusted all conductance densities to reproduce voltage-clamp and current-clamp data from heart interneurons, focusing on five ‘fit metrics’, described below, for comparing the model neuron to living heart interneurons.
We chose fit metrics that tested both passive and active (slow and fast) properties of the model in different functional regions, that were relevant to the activity of the neuron in isolation and in networks, and that were experimentally measurable. Thus, we sought metrics that would constrain different combinations and distribution of conductances. The first fit metric was to evaluate current during somatic voltage-clamp ramps (over the range -40 to -90 mV) in leech physiological saline containing Cs+ and bicuculline. Bicuculline blocks synaptic transmission, effectively isolating the neurons (Schmidt and Calabrese 1992) and Cs+ blocks Ih (Angstadt and Calabrese 1989) to reveal a linear portion of the voltage ramp current that reflects the leak current, ILeak; these blocks were simulated by setting both h and synaptic conductances to 0. The non-linear regions of the voltage ramp current are largely determined by low threshold voltage-gated currents, such as IP, ICaS, ICaF, IK1, IK2, and IKA. The second fit metric was action potential shape as recorded in the soma, which is determined primarily by fast high-threshold currents, primarily INa and IK1, and by capacitance. Low threshold currents such as IP, ICaS, ICaF, and IKA shape the baseline membrane potential; the spike threshold and amplitude are largely dependent on INa, and the afterhyperpolarization is shaped by the potassium currents, IK1, IK2, and IKA. The third fit metric was the spike frequency response to constant injected current in the soma, both hyperpolarizing and depolarizing. Parameter variation in the single-compartment heart interneuron model indicates that IP, ICaS, IK2, and ILeak can all strongly affect the spike frequency (Hill et al. 2001). The fourth and fifth fit metrics were synaptic current amplitude and reversal potential as measured in somatic voltage-clamp. Both are substantially affected by axial resistance and the balance of voltage-gated currents in the secondary neurites that affect passive spread and space-clamp during somatic voltage-clamp. These 5 experimentally measurable fit metrics were chosen to provide comprehensive constraints on conductance densities for tuning the model to network-relevant activity. Because heart interneurons (soma diameter approx. 20 μm) degraded after 10 to 15 minutes of intracellular recording, and even faster during bathing in low Ca2+ saline or other pharmacological manipulations, we did not obtain data for all fit metrics from a single interneuron. Thus, we tuned our model to match typical responses collected from many neurons. Because changing each conductance affected several fit metrics, tuning the model required an iterative process of hand-tuning conductance densities and evaluating the match to the experimental data.
We measured the current during hyperpolarizing voltage ramps (first fit metric) in interneurons bathed in saline with bicuculline and Cs+. As described in Methods, passive properties were set to match those extracted from the linear portion (between -70 and -80 mV) of the voltage ramp current from interneurons; we tuned active conductances to match examples of the voltage ramp current, both quantitatively and by eye (Fig. 3). When using the potassium conductance kinetics as published in Hill et al. (2001), the current measured in the model exhibited an outward rectification during the negative-slope leg of the voltage ramp (not shown), opposite to the inward rectification seen in experimental data. Adjusting conductance densities did not significantly affect this concavity, so we shifted the activation curves of the potassium currents (IK1, IK2, IKA) by -10 mV, corresponding to those implemented in an earlier version of the single-compartment model of a leech heart interneuron (Nadim et al. 1995). This shift removed the outward rectification in the descending leg of voltage ramp but did not introduce the inward rectification observed experimentally (Fig. 3), which is apparently due to deactivation of IP (Olsen and Calabrese 1996). Additionally, heart interneurons show a pronounced inward current activating at the end of the depolarizing voltage ramp, due to low-threshold Ca currents (ICaS and ICaF) (Olsen and Calabrese 1996). Increasing the conductance of either low-threshold Ca currents or IP to ameliorate these flaws in the voltage ramp current worsened the match to other fit metrics. Therefore, we accepted the imperfections of fitting the voltage ramps as a compromise to the other fit metrics. The simulated voltage-clamp current of the Full Model, with the finally chosen parameter set, matched experimental voltage ramp currents (Fig. 3), particularly in holding current and slope in the linear region. Fitting the hyperpolarizing portion of the ramp current between -70 and -80 mV yields an input resistance of 128 MΩ for the Full Model and 114 MΩ for the displayed experimental example.
Examples of somatically recorded action potentials (second fit metric) were chosen from current-clamp experiments (example shown in Fig. 4). We adjusted conductance densities to match experimental data, namely the amplitude (quantitatively) and the time course (by eye) of the action potential and its afterhyperpolarization. In particular, to set action potential amplitude, we tuned the conductance of INa, and to match the threshold and time course, we shifted the activation and inactivation curves of INa by -2 mV. The action potential peak of the Full Model, using the finally chosen parameter set (-7 mV), was within observed variation of action potentials during tonic spiking in living heart interneurons (experimental example: -9 mV) (Fig. 4). The repolarization of the action potential was noticeably slower in all versions of the model than in the heart interneurons (Fig. 4). Fitting the action potential shape (amplitudes and time course) strongly influences the spike frequency response to injected current (third fit metric) by influencing the action potential's threshold characteristics. The model neuron spike frequency matched experimental data well for zero current injection (experimental average: 6.6 ± 1.2 Hz; Full Model: 6.2 Hz) (Fig. 5). The spike frequency response to injected current of the Full Model was slightly lower than that of heart interneurons for the most positive current injections and did not exhibit a smooth transition between silence and spiking. Despite this abruptness, the model can fire at low enough frequency to capture normal bursting activity, where minimum spike frequency is around 4 Hz (Hill et al.2001) and so was judged acceptable.
In heart interneurons, spike-mediated synaptic current (fourth fit metric) was measured with somatic voltage-clamp in normal physiological saline at a holding potential of -50 mV and was found to have an average peak of 190.1 ± 30.4 pA (Tobin and Calabrese 2005). The synaptic reversal potential in heart interneurons (fifth fit metric) has been measured most accurately for graded synaptic transmission, using somatic voltage-clamp in 0 Na+ saline while K currents were blocked, yielding a reversal potential of -62.5 mV (Angstadt and Calabrese 1991). We simulated these experimental protocols, and set the model synaptic conductance to produce a current peak amplitude measurement of 190.0 pA (Fig. 6). The synaptic reversal potential in the Synaptic Compartments was set to the empirically determined value of -62.5 mV. Passive cable properties and unblocked active membrane conductances (only ICaS and ICaF remain) could compromise space-clamp and therefore alter somatic measurement of reversal potential. We measured synaptic reversal potential in the model soma near -65 mV, sufficiently close to the set value of -62.5 mV (Fig. 6). The similarity between somatically-measured and set values indicates reasonably good somatic space-clamp of synaptic conductance in the model and that the model neuron is electrically compact, at least when voltage-clamped near -60 mV with most active conductances blocked.
The parameter set finally chosen for the Full Model (see Table 2) was the best of those tried in reproducing realistic activity of a single isolated heart interneuron using the fit metrics described above. Conductance amplitudes, as measured by simulated voltage-clamp, were within 14% of initially chosen values, except ICaS, which was increased by 42% (Table 3). These variations are within experimentally observed variability, including the increased ICaS conductance (Tobin and Calabrese 2005).
Through extensive hand tuning of conductances, we produced a satisfactory model of a single leech heart interneuron. Based on this morphologically-detailed Full Model, we wished to create a more computationally efficient reduced model for further parameter searches and for use in network simulations. The structure of our Full Model with functional and transitional compartment types, each with uniform conductance densities, lends itself to reducing compartment number while preserving the compartment types so that the reduced model can be mapped onto the Full Model.
In developing a reduced model, we first sought to decrease the simulation time by reducing the number of compartments. For a reduced model and the Full Model to share parameter settings (conductance densities and distributions), we reasoned that it would be necessary to preserve surface area of each functional/transitional region. Several reduction methods have been developed that preserve surface area and focus on reducing complex branching structures into a cable-like series of compartments (Burke 2000; Clements and Redman 1989; Rall 1959; Rall 1964). However, we desired an even more drastic reduction, representing the entire Secondary Neurite structure by only one or two compartments and the already cable-like axon by a single compartment. Destexhe (2001) has developed a method of reducing a complex branching structure to 3 compartments, while maintaining surface area, passive somatic voltage response to injected current, and somatodendritic attenuation of the somatic response. We began our reduction using similar methods by reducing all compartments within each functional/transitional region into one representative, reduced compartment with the same total length and surface area as the set of compartments it represents. As with Destexhe's method, the axial resistances of the reduced compartments were chosen to fit passive voltage responses to somatically injected current. Instead of fitting all axial resistances simultaneously, however, we chose a sequential process, combining each functional/transitional set of compartments into one compartment and scaling the axial resistance of that compartment before reducing the next functional/transitional set of compartments. In this fashion, we ensured that each axial resistance was appropriately set.
We started by reducing the 132 Axon compartments into a single, reduced Axon compartment with the same diameter and length as the average diameter and total summed length of the 132 full Axon compartments. These dimensions were chosen to match the surface area between the reduced compartment and the compartments it represented in the Full Model. Conceptually, this reduced compartment was equivalent to separating the entire membrane area of the Axon from the rest of the model by a large resistor, with resistance equal to the total summed axial resistance of the entire 485 μm of the Axon length. This substitution was reflected in an increase in total neuron input resistance from 105 MΩ to 114 MΩ (Fig. 7). To preserve the total neuron input resistance, we reduced the axial resistance (by hand) of the reduced Axon compartment until the total neuron input resistance, as measured by current pulses in the soma, matched the input resistance of the Full Model with the full Axon (105 MΩ). Any further reduction of axial resistance reduced the total neuron input resistance below 105 MΩ. The resulting axial resistance was 0.145 times the total summed axial resistance of the Full Model Axon. Scaling the axial resistance of a compartment is equivalent to scaling the length and diameter, with the constraint of constant surface area. Axial resistances of compartments are calculated using a global axial resistivity value, Ra, in units of Ωcm, scaled by the length, l, and diameter, d, of each compartment as follows:
axial resistance = 4lRa/(πd2)
To maintain the convention of using a global axial resistivity for all compartments in the model, we scaled the length and diameter of the Axon compartment, such that the new dimensions maintain surface area and produce the scaled value of axial resistance. The dimensions are given in Table 4 and shown in Fig. 1B.
Using the same method, we sequentially reduced all compartments of each functional/transitional region into a single reduced compartment (except for the Neurite which was divided into two equivalent compartments – see below), initially with the same length, average diameter, and total surface area of the represented compartments. The axial resistance of the reduced compartment was then scaled to preserve somatically measured input resistance, and the geometry of the compartment was adjusted appropriately. For example, the four Neurite-axon2 compartments were reduced to one. Model input resistance did not change measurably and thus no axial resistance scaling was required. In contrast, reducing the 10 Neurite-axon1 compartments and then the 2 Neurite-soma compartments required scaling the axial resistances by 0.22 and 0.953, respectively.
The highly branched secondary neurite structure did not permit the same sequential reduction of the Secondary Neurite and Synaptic Compartments as we used for the cable-like Neurite-axon and Axon sections. We first removed all 48 Synaptic Compartments and measured an input resistance of 108 MΩ. The 663 Secondary Neurite compartments were reduced to one compartment, whose axial resistance was scaled by 0.0004 to preserve the 108 MΩ input resistance corresponding to the Synaptic-Compartment-less model. Then, we added a reduced Synaptic Compartment, whose axial resistance was scaled to produce an input resistance of 105 MΩ corresponding to the Full Model. Finally, the 11 Neurite compartments were reduced to two equal compartments to allow the Secondary Neurite compartment to branch from between these two Neurite compartments. The axial resistance of the Neurite compartments was scaled by 1.28 to maintain input resistance of 105 MΩ.
All compartments of the Full Model were represented in the reduced model, maintaining the surface area of each functional/transitional region and somatically measured input resistance. In addition, the membrane time constant was preserved; τm of the Full Model was 32.4 ms whereas τm of the reduced model was 31.9 ms, only a 1.5% decrease. Because the compartment reduction and scaling was done to preserve passive responses, we refer to this model geometry as the ‘Passively Reduced Model’ (Fig. 1B).
We applied the active and synaptic conductances densities from the Full Model to the Passively Reduced Model and compared the Full and Passively Reduced Models in simulated voltage-clamp and current-clamp conditions. We compared the output of the models according to the fit metrics we used for tuning the Full Model: voltage-clamp current during hyperpolarizing voltage ramps, action potential shape, spike frequency response to injected current, and synaptic current amplitude and reversal potential.
The Passively Reduced Model replicated the Full Model well for hyperpolarizing voltage ramps (Fig. 3), spike frequency response to injected current (Fig. 5) and synaptic current reversal potential (Fig. 6). The synaptic current amplitude and the action potential shape, however, did not match. Synaptic current peak amplitude is 190.0 pA in the Full Model, but only 55.3 pA in the Passively Reduced Model (Fig. 6), and the action potentials had shallower after-hyperpolarization potentials than occur in the Full Model (Fig. 4). There could be several explanations why synaptic currents (IPSCs) and action potentials were particularly distorted by the morphological reduction. The axial resistances of the Passively Reduced Model were scaled to preserve input resistance measurements, but somatic input resistance was relatively insensitive to changes in the axial resistance of particularly distal compartments in the Passively Reduced Model. For example, when the Full Model Axon was collapsed to a single compartment, the total input resistance of the model became 114 MΩ. To restore the input resistance to 105 MΩ, a decrease of only 8%, the axial resistance of the Axon was decreased by 85%. As large variations in axial resistance produced only small variations on total input resistance, matching input resistance within approximately 1% may not have been a sufficiently strict constraint on axial resistances of particularly distal compartments. Additionally, because the action potentials are initiated in the Neurite-axon transitional compartments, and conduct through the Neurite and Neurite-soma to the soma with attenuation and altered time course due to the progressive decrease in INa density, their shape, as measured in the soma, may be particularly altered by the axial resistance settings. Moreover, synaptic events and action potentials may be particularly affected because of the severe morphological distortion of the compartments where these events originate. Additionally, the distal-to-proximal reduction sequence matched the current load of distal regions on proximal ones, but not of proximal regions on distal ones, where synaptic events and action potentials are initiated. We reasoned that the Passively Reduced Model would better replicate the Full Model active properties if the axial resistances of the Secondary Neurite, Synaptic Compartments, and Axon were rescaled with a focus on improving the synaptic currents (IPSCs) and action potentials.
To reduce attenuation of the synaptic current, we started by rescaling the axial resistances of the Secondary Neurite and Synaptic Compartments. Because much of the Ca current in the model is concentrated in the Secondary Neurite and Synaptic Compartments, we monitored whether rescaling the axial resistances of these compartments would affect the Ca current waveform, measured by simulated somatic voltage-clamp. The peak Ca current of the Passively Reduced model was slightly greater than the Full Model, but decreasing the axial resistances of the Secondary Neurite and Synaptic Compartments improved the match between the two models (Fig. 8). Rescaling the Secondary Neurite axial resistance by 0.8 and the Synaptic Compartment axial resistance by 0.2 produced the best match of Ca current waveforms (matched by eye) while increasing synaptic current amplitude to 185.8 pA, 2.2% less than that of the Full Model value (Fig. 6).
To produce action potentials more similar to the Full Model, we scaled the axial resistance of the Axon by 0.03. We arrived at this value as a compromise between fitting the trough of the afterhyperpolarization potential and the rise of the membrane potential prior to the spike. Further reduction of the Axon axial resistance brought the membrane potential 50 ms prior to the spike closer to that of the Full Model but also made the afterhyperpolarization potential too hyperpolarized. The resulting model morphology, with rescaled axial resistances of the Secondary Neurite, Synaptic Compartment and Axon, is referred to as the ‘Actively Reduced Model’ (Fig. 1C; Table 5). As the Passively Reduced Model was developed to have identical input resistance to the Full Model, rescaling the axial resistances to produce the Actively Reduced Model decreased input resistance to 99 MΩ, 5.7% less than the Full Model input resistance (Fig. 7, green line).
Finally, we compared the Actively Reduced Model to the Full Model for our five fit metrics. As mentioned above, the Actively Reduced Model has a lower input resistance than the Full Model, resulting in a steeper current waveform during hyperpolarizing voltage ramps in voltage-clamp (Fig. 3), but the differences are slight and satisfactorily within experimentally observed variability. As action potential shape was a determining factor in developing the Actively Reduced Model, the action potential was similar to those of the Full Model and the living heart interneurons except that all the model action potentials repolarized more slowly (Fig. 4). The spike frequency response to injected current of the Actively Reduced Model is similar to that of the Full Model, deviating at most by 6.3% at 0.2 nA current injection (Fig. 5). Although the synaptic current amplitude decreased slightly when the Axon compartment was rescaled, likely due to increasing the effect of the Axon as a current sink, the final current peak amplitude of the Actively Reduced Model, 182.3 pA, represents only a 4.0% decrease from that of the Full Model, 190 pA. In addition, the measured reversal potential of the synaptic current was similar in both the Full and Actively Reduced Models, near -65 mV (Fig. 6).
In all chosen fit metrics, the Actively Reduced Model satisfactorily matched the Full Model. Thus, we consider the Actively Reduced Model a useful representation of the Full Model for network simulations where we can sacrifice some morphological detail for computational efficiency. Moreover, it can be used in parameter searches for various activity regimes, such as endogenous bursting of which the oscillator heart interneurons are capable (Cymbalyuk et al. 2002). The results of these searches can be systematically mapped back onto the Full Model to help locate similar activity regimes. The computational efficiency gains are substantial: for example, our simulation computer required 29 cpu seconds to simulate 1 second of activity in the Full Model (871 compartments and 9 active conductances totaling 5772 simulated conductances), but only 0.18 cpu seconds for the Actively Reduced Model (9 compartments and 9 active conductances totaling 62 simulated conductances).
To explore how conductance distribution influences the electrical activity of heart interneurons, we developed a detailed morphology model, the Full Model (Fig. 1A). To reduce parameter space for tuning, we identified the functional regions of the interneuron, (Soma, Neurite, Axon, Secondary Neurite, and Synaptic Compartment), as well as several transitional regions (Neurite-soma, Neurite-axon1 and Neurite-axon2), and required homogeneous conductance densities in each of these regions. We tuned conductance densities to produce a realistic heart interneuron model with detailed morphology. To enhance computational-efficiency for parameter variations and network simulations, we simplified the morphological complexity; reducing nearly 900 compartments to 9. Our method, while similar to that developed by Destexhe (2001), reduces sections of the Full Model individually and sequentially, to ensure a unique match between the current load of each section of the Full Model and that of the corresponding reduced compartment. Additionally, our method readjusts the reduced model to match the active properties of the Full Model. As a result, we achieve two models whose passive and active properties are closely matched, but whose different computational efficiencies and levels of morphological detail make each uniquely suited for different uses.
We first assigned values to the active conductances based on available somatic voltage-clamp data (Table 3). Further tuning of the Full Model and the subsequent reduced models was directed by fitness criteria (fit metrics), largely based on passive properties, synaptic currents and spiking activity (see Figs. 3, ,4,4, ,5,5, ,6),6), and the choice of these metrics ultimately determines the utility of the model as experimental and theoretical tool. First and foremost, our choice of fitness metrics was dictated by measures that are experimentally feasible and replicable in our system. Morphological models definitely benefit greatly from recordings at multiple points within a neuron's extended structure (Keren et al. 2005), however, this technique is not currently feasible in heart interneurons, so we were constrained to metrics obtainable from somatic recordings. Fitness criteria were also chosen because they were relevant to the cell activity and tested the spatial distribution of conductances. For example, synaptic input and spike initiation are distributed distally within the heart interneuron structure and thus synaptic reversal potential and spike amplitude/shape and spike frequency response to injected current constrain both passive properties and active conductance distribution choices. These metrics led to substantial changes from the initial tuning to somatic voltage-clamp data (Table 3). Given the expected 2-4 fold range of active conductances, as observed in a variety of neuronal networks (Golowasch et al. 2002; Schulz et al. 2006; Swensen and Bean 2005) these final tuned values fall acceptably close to the exemplars provided from the literature (Table 3).
The tuning process and the final conductance distribution set chosen (Table 2) provide some insights into the electrotonic structure and spiking activity of heart interneurons. The neurons are relatively compact to synaptic inputs (somatically measured synaptic reversal potential differs by only a few mV from that set in the distal synaptic compartments). As might be expected by their role in regulating excitability, outward currents dominate in all compartments. The neurite is the focus of sub-threshold inward currents (Ih, IP, ICaS) that drive spiking activity, and it is poised to receive synaptic current from the secondary neurites and feed current to the spike initiation site in the Neurite-axon2. Spiking activity is well captured in the model with realistic somatically-recorded spike amplitudes and realistic spontaneous spike frequency. Spike initiation occurs in the Neurite-axon2 region, and both the peak of the action potential and the trough of the afterhyperpolarization potential attenuated greatly from the initiation zone to the somatic recording site (data not shown). Although the spike frequency response to injected current in the model captures the near linearity observed over a large range of currents, the transition to spiking is more abrupt than observed in the living neurons and may indicate that the model does not capture the proper dynamics of this transition (Rinzel and Ermentrout 1998). This discrepancy may not be a major concern, however because the model is capable of firing at low frequency (~3 Hz, Fig. 5A) and the normal bursting activity of heart interneurons (~4 Hz minimum intraburst spike frequency, ~12 Hz average intraburst spike frequency, (Hill et al. 2001)) may not require exact capture of the transition to spiking. Nevertheless, a re-evaluation of the kinetics of INa in the model may be warranted (Naundorf et al. 2006).
Our multi-step reduction method was facilitated by designing the Full Model using functional/transitional regions (Fig. 1A), which were constrained to have homogeneous conductance densities within each region. The functional regions correspond to areas we determined to be anatomically and physiologically distinct: soma, neurite, secondary neurite, synaptic contacts and axon. Within each region, active and synaptic conductance densities are uniform but differ from region to region, reflecting experimentally determined differences and presumed function. Such distinctions would be specific to the cell type and model goals. For example, Destexhe (2001), wishing to describe synaptic integration within a passive dendritic tree, divided a pyramidal cell into the three regions of soma/proximal dendrites, basal dendrites, and apical dendrites (including distal ends of basal dendrites). Non-homogeneous conductance densities are often distributed in an approximately linearly changing fashion (Bekkers 2000; Hoffman et al. 1997; Korngreen and Sakmann 2000; Magee 1999; Stuart and Hausser 1994). In our models, transition compartments, with intermediate conductance densities (see Methods) link functional regions with widely different active conductances. They roughly approximate linear changes in conductance densities and roughly match impedances to prevent reflection of propagating signals. Because the conductance densities in the transition regions are constrained by the conductance densities in abutting compartments, these compartments do not add to the parameter space of the model.
In the first reduction step, we collapsed the compartments of each functional/transitional region of the Full Model into a single reduced compartment (2 for the Neurite). We matched the current load of each reduced compartment to that of the represented functional/transitional regions in the Full Model by rescaling the axial resistance. We chose the initial length of each reduced compartment to be equal to the summed length of the compartments it represented. In his reduction procedure, Destexhe (2001) chose each reduced compartment to have the typical length of the portion of the dendrite it represented. The end results of these two approaches are similar, as the axial resistances, and thus the lengths of the compartments were rescaled to match passive properties. The methods differ, however, in that Destexhe's method (Destexhe 2001) fits the voltage decay of a somatic current pulse with distance from the soma, while our method relies on maintaining each functional region's contribution to the neuron's input resistance, as measured at the soma. Additionally, his method uses an algorithm to fit all axial resistances simultaneously, while we hand fit each axial resistance sequentially. Compartments were reduced sequentially to avoid the possibility of non-unique solutions, as the voltage spread to one compartment is affected by the axial resistances of other compartments in the model. Non-unique solutions may be particularly problematic in a model with more compartments, such as our 9 compartment model, as opposed to Destexhe's 3 compartment model (Destexhe 2001). By sequential reduction we ensure that each reduced compartment matches the current load of the compartments it represents before reducing the next set of compartments. Thus, as the current load is proportional to axial resistance, there is a unique value of axial resistance to match each set of compartments' current load.
The reduction order may influence the axial resistances chosen for the compartments. We chose a distal-to-proximal sequence, as we were assaying the current load on the soma and in this way the compartments being reduced were receiving current through non-reduced compartments. This does not necessarily preserve the current load in the reverse order (of proximal compartments on distal ones), but our subsequent readjustment, during formation of the Actively Reduced Model, to preserve attenuation of synaptic inputs and axon potential shape, likely improves the current load of proximal compartments on the distal compartments where these events are initiated, albeit at some cost to accurate representation of the current load on the soma as seen in the 5.7% decrease in Rin (Fig. 7)
While the simplicity of our method allows for easy hand-fitting of axial resistances, the passive reduction method could be automated by a program that sequentially calculates surface area of each functional group, replaces functional groups of compartments with a reduced compartment of equivalent surface area, and adjusts the axial resistance with an automated bisection method to preserve the neuron input resistance.
In the second reduction step, we applied the conductance densities from the Full Model to the Passively Reduced Models. We matched the active properties of the Passively Reduced Model to those of the Full Model by rescaling axial resistances, producing the Actively Reduced Model. This rescaling is in part necessary to correct the biases introduced into the model by the distal to proximal sequential nature of the passive reduction so that the effective load of proximal compartments on more distal ones is taken into account for events such as spikes and IPSPs that are generated in distal compartments. The resultant model is thus a compromise between one that faithfully reproduces events generated proximally (e.g. soma current injection) and those generated distally (spikes and IPSPs).
The electrotonic structure of a compartment can affect the activity in all other compartments by affecting its current load on those compartments. Thus, rescaling the axial resistances of more than one compartment in the Passively Reduced Model presents a circular process. For example, after rescaling the axial resistance of the Secondary Neurite and Synaptic Compartments of the Passively Reduced Model to match measured Ca currents (Fig. 8) and synaptic current amplitude (Fig. 6), we then had to rescale the axial resistance of the Axon to improve the match of the action potential shape to that of the Full Model. Rescaling the Axon decreased the amplitudes of the measured synaptic current (data not shown) and the Ca current waveform (Fig. 8). In our case, the changes in synaptic and Ca currents were small enough that we did not rescale any compartments, however, this potential problem highlights the need to recheck all metrics after all compartments are rescaled.
An automated process such as discussed above for the passive reduction could be extended to fit axial resistances to match fit metrics, starting with the axial resistance values obtained during passive reduction. Or, potentially, axial resistances could be chosen to best fit neuron input resistance and active properties simultaneously. Our emphasis here, which is unique to our approach, is to use sequential reduction to avoid non-unique solutions and, for neuron models with active structures such as axons and active dendrites, to include active properties in addition to passive properties in fitting the reduced model geometry.
While previous methods have been developed to simplify neuronal morphology, such as Rall's equivalent cable method (Agmon-Snir 1995; Burke 2000; Clements and Redman 1989; Rall 1959; Rall 1964; Zador et al. 1995), all such solutions, to our knowledge, are not applicable to structures with active conductances in which the electrotonic structure of the neuron is a function of voltage and time, not simply surface area and distance. The method we propose here, though necessarily empirical and approximate, is applicable to any neuronal structure with any set of conductances.
We have succeeded in tuning the Full Model to the intracellularly recorded voltage-clamp data and tonic spiking characteristics from synaptically isolated leech heart interneurons. Future applications of this model for analysis of normal cellular and network activity may require adjustment of some conductance densities that are likely critical for endogenous and half-center bursting, such as ICaS, Ih, ISyn, ELeak and Rm (Cymbalyuk et al. 2002; Hill et al. 2001) that are not as constrained during tonic spiking activity. Because the Actively Reduced Model is more computationally efficient, we will use it for exploring the parameter space where endogenous bursting and bursting in a half-center occurs.
Having two similar models with different levels of morphological detail allows us to choose the appropriate model based on the question we wish to pursue. The Full Model was constructed to investigate how morphology and intrinsic and synaptic conductance distributions affect cellular and network activity. However, it is a difficult tool to wield because it is computationally demanding. The Actively Reduced Model is most appropriate for studies where morphological detail can be sacrificed for computational efficiency, such as network simulations, for extensive tuning and exploring parameter space. Because the morphology and intrinsic and synaptic electrical properties of the two models are directly interchangeable, information gained from one can be applied to the other, providing us with a comprehensive tool for understanding the effect of intrinsic and synaptic conductance distributions on the activity of leech heart interneurons.
We thank Robert Cannon for providing java code for CVAPP. We thank Angela Wenning for critically reading the manuscript and John Rinzel for helpful discussions.
Grants: This work was supported by the National Institute of Health NRSA predoctoral fellowship NS-42983 to A-E. Tobin and by National Institute of Neurological Disorders and Stroke Grant NS-24072 to R. L. Calabrese.