HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) All Versions of this Article: 96/4/2089    most recent 00025.2006v1 Alert me when this article is cited Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via ISI Web of Science (2) Citing Articles via Google Scholar Google Scholar Articles by Tobin, A.-E. Articles by Calabrese, R. L. Search for Related Content PubMed PubMed Citation Articles by Tobin, A.-E. Articles by Calabrese, R. L.

Endogenous and Half-Center Bursting in Morphologically Inspired Models of Leech Heart Interneurons

Department of Biology, Emory University, Atlanta, Georgia

Submitted 10 January 2006; accepted in final form 1 June 2006

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES

 
Based on a detailed morphology "Full Model" of a leech heart interneuron, we previously developed a computationally efficient, morphologically inspired "Reduced Model" to expedite tuning the model to produce endogenous bursting and alternating bursting when configured as a half-center oscillator (paired with reciprocally inhibitory synapses). To find conductance density distributions that produce endogenous bursting, we implemented a genetic algorithm automated parameter search. With multiple searches, we found eight parameter sets that produced endogenous bursting in the Reduced Model. When these parameter sets were applied to the Full Model, all produced endogenous bursting, although when the simulation time was extended from 80 to 300 s, only four parameter sets produced sustained bursting in the Reduced Models. All parameter sets produced alternating half-center bursting in the Reduced and Full Models throughout the entire 300 s. When conductance amplitudes were systematically varied for each of the four sustained burster sets, the effects on bursting activity differed, both for the same parameter set in the Reduced and Full Models and for different parameter sets with the same level of morphological detail. This implies that morphological detail can affect burst activity and that these parameter sets may represent different mechanisms for burst generation and/or regulation. We also tested the models with parameter variations that correspond to experimental manipulations. We conclude that, whereas similar output can be achieved with multiple different parameter sets, perturbations such as conductance variations can highlight differences. Additionally, this work demonstrates both the utility and limitations of using simplified models to represent more morphologically accurate models.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES

 
To understand network function in nervous systems, the neurons composing the networks must be identified and characterized at the cellular and subcellular levels. To this end, the most progress has been made in invertebrate motor pattern generating networks. Their few and individually identifiable neurons enable identification of the component neurons in these networks and characterization of their synaptic connections and membrane currents (for reviews, see De Schutter et al. 2005; Harris-Warrick 2002; Marder and Calabrese 1996; Marder et al. 2005). The plethora of data on intrinsic and synaptic conductances and their kinetics from detailed electrophysiological studies has necessitated modeling to elucidate how intrinsic and synaptic conductances interact to shape neuronal and network electrical activity (e.g., Hill et al. 2001; Jezzini et al. 2004; Marder et al. 2005; Nadim and Manor 2000; Nadim et al. 1998). Although such models have been very successful in promoting analysis and understanding of neuron and network activity, realistic morphological details of the component neurons have not been included. This omission is explained partly by the computational costs of detailed compartmental models of neurons and partly by a lack of detailed information about how conductances are distributed within a neuron. Yet structure does matter for neuronal and thus network function. For example, Traub and colleagues (1991) produced an eight-conductance, 19-compartment model of a CA3 hippocampal neuron to analyze several neural activities such as generation of dendritic calcium spikes, depolarizing afterpotentials following spikes, and rhythmic bursting with injected current that simply were not observable in single-compartment models containing the same conductances. This model was subsequently reduced to two compartments, maintaining most of the model activity by segregating fast and slow conductances in two electrically coupled compartments for computational efficiency and dynamical analysis (Pinsky and Rinzel 1994). These studies demonstrate how conductances and their distributions may influence neural activity and provide an inspiration for our exploration of effects of neuronal structure on activity in the heartbeat pattern generating network of the leech.

The network driving the rhythmic heartbeat of the leech is a pattern generating network where not only the neurons and their connections (Thompson and Stent 1976b,c) but also the rhythmic motor output (Thompson and Stent 1976a; Wenning et al. 2004) have been well characterized. Within this central pattern generating network, consisting of seven bilateral pairs of heart interneurons, the interneuron pairs in ganglia 3 and 4 are each linked with reciprocally inhibitory synaptic connections, forming two half-center oscillators that produce alternating bursts of action potentials and pace activity in the network (Peterson 1983). The alternating bursting of each of these oscillators is maintained when their ganglia are isolated (Masino and Calabrese 2002b; Peterson 1983). Extensive electrophysiological characterization of the ionic conductances in the oscillator interneurons has led to the development of a conductance-based single-compartment model of the oscillator interneurons (Hill et al. 2001). An earlier generation of this model (Nadim et al. 1995) guided experiments investigating which conductances are targeted by the endogenous leech peptide FMRFamide to accelerate the rhythmic bursting of the oscillator heart interneurons (Nadim and Calabrese 1997).

Building on our single-compartment model (Hill et al. 2001), we developed a morphologically detailed "Full Model" (871 compartments) of an oscillator heart interneuron, described in the companion paper (Tobin et al. 2006), that permits localized membrane potential changes (such as synaptic inputs) and nonhomogeneous conductance distributions. This hand-tuned model successfully captures the tonic spiking activity of pharmacologically isolated heart interneurons recorded with sharp microelectrodes, but does not express endogenous bursting activity of which the living neurons are capable (Cymbalyuk et al. 2002).

For computational efficiency in network simulations and in parameter searches we reduced the compartment number (nine compartments), preserving passive and active properties and the functional regions. This Reduced Model maps in a systematic way onto the Full Model so that conductance densities chosen for the Reduced Model can be applied to the Full Model (Tobin et al. 2006). Here we use this model for criterion-directed retuning that can be directly applicable to the Full Model. Because of large parameter spaces of active-membrane multicompartmental models, parameter optimization algorithms can be extremely useful for finding parameter sets producing desirable model activity. Vanier and Bower (1999) demonstrated that genetic algorithms are among the most effective parameter-search methods for large parameter spaces [23 parameters in Vanier and Bower (1999)]. This algorithm has been used by others to parameterize compartmental neural models (Hobbs and Hooper 2003; Keren et al. 2005).

Our primary goal in the research described here was to recover parameter sets (conductance density distributions) from our parameter searches that capture endogenous bursting and half-center bursting in our Reduced and Full Model neurons, so that we might study the effects of neuronal morphology and conductance distribution on endogenous and network bursting activity. We implemented an automated genetic algorithm parameter optimization routine to find multiple parameter sets that produced endogenous bursting in a single-neuron Reduced Model. We recovered four sets that produced desired activity and requisite stability, and that were remarkably different, suggesting that different conductance distributions can lead to similar neuronal and network activity, as has been found in crustacean model neurons and networks (Prinz et al. 2003, 2004). The Full Models, both single-cell and half-center, generated from these parameter sets, were similar in activity to their corresponding Reduced Models but more stable to parameter variation, suggesting that neuronal morphology can stabilize endogenous and network activity. We then sought to determine whether, for each of these parameter sets, the models responded appropriately to parameter changes associated with known neuromodulatory effects. Our findings indicate that challenging models with such modulatory influences can be useful in distinguishing different parameterizations of morphological models.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES

 
Model parameters

The morphology, compartment functional/transitional groups, and passive parameter settings of the Full Model and Actively Reduced Model, heretofore called the "Reduced Model," are described in the companion paper Tobin et al. (2006). A specific axial resistivity of 250 cm and specific membrane capacitance of 2.2 µF/cm², chosen as described in Tobin et al. (2006), were used for all simulations described here. Leak reversal potential E_Leak and specific membrane resistance R_m were varied with the parameter optimization routine, as described below. We modeled nine active membrane conductances, as previously identified in leech heart interneurons using voltage-clamp procedures (I_Na was identified by current clamp). There are five inward currents: a fast Na current (I_Na), a persistent Na current (I_P), a rapidly inactivating low-threshold Ca current (I_CaF), a slowly inactivating low-threshold Ca current (I_CaS), and a hyperpolarization-activated cation current (I_h). There are three outward currents: a fast, transient K current (I_KA), a delayed-rectifier–like K current (I_K1), and a persistent K current (I_K2). Additionally, we included a spike-mediated synaptic current (I_Syn), with reversal potential –62.5 mV, as described in Hill et al. (2001). In brief, this postsynaptic conductance is triggered by presynaptic spikes and is modeled as the sum of a rising (2 ms) and a falling (11 ms) exponential. The peak of the synaptic conductance is modulated by presynaptic voltage, in accordance with experimental findings. The kinetics and reversal potentials of all the currents are as described in Hill et al. (2001), except for modifications noted in the companion paper (Tobin et al. 2006).

Conductance distributions were assigned to best approximate physiological data and to limit the parameter space for tuning the model. We restrict the Soma to only I_K1, I_K2, and I_KA and the Axon to only I_Na, I_K1, I_K2, and I_KA. In contrast to the Secondary Neurite conductance distributions implemented in the models described in Tobin et al. (2006), the Secondary Neurite compartments of these models have I_h and I_P in addition to I_K1, I_K2, I_KA, I_CaF, and I_CaS. Densities of the Synaptic Compartments are fixed to be identical to those in the Secondary Neurite, except that the synaptic conductance is present only in the Synaptic Compartments and was fixed to the conductance density of 100 mS/cm². The transitional compartments, Neurite–soma, Neurite–axon1, and Neurite–axon2, were designed to smooth transitions between compartments of adjacent functional regions that have different conductance densities. They roughly approximate linear changes in conductance densities that are often distributed in a linear fashion (Bekkers 2000; Hoffman et al. 1997; Korngreen and Sakmann 2000; Magee 1999; Stuart and Hausser 1994) and roughly match impedances to prevent reflection of propagating signals. The Neurite–soma conductance densities are an average of those in the Soma and Neurite. The Neurite–axon1 and Neurite–axon2 have the same density of K conductances as that of the Neurite. For I_h and I_P, the Neurite–axon1 conductance densities are the sum of two thirds of those in the Neurite; the Neurite–axon2 densities are one third of those in the Neurite (the Axon has neither I_h nor I_P). To set the spike initiation zone, fast Na density is set as one half of Axon density in the Neurite, three quarters in Neurite–axon1, and one in Neurite–axon2 to ensure action potentials initiate in the Neurite–axon2 compartments. Experimentally, spontaneous action potentials in heart interneurons are not elicited 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).

Parameter optimization

We used a genetic algorithm parameter optimization routine implemented in MATLAB (The MathWorks, Natick, MA) by Houck and colleagues (1995). In brief, the genetic algorithm works by creating an initial population of parameter sets where the value of each parameter in a set is randomly chosen (within defined bounds). Each parameter set is rated (assigned a fitness value) by comparing the activity of the model with the parameter set against the desired activity (fitness criteria). In an evolution-inspired fashion, new parameters sets are created by "breeding" the highest-rated parameter sets. The new offspring are then rated, and the best-rated sets are bred to produce new offspring. The breeding continues throughout either a specified number of generations or until the fitness value passes a threshold. The specifics of the options used in calling the genetic algorithm MATLAB command are given in the Supplemental Methods, "Genetic Algorithm Command Options." Because of our large number of parameters, we chose a large initial population size of 10,000 unique parameter sets. Instead of specifying a convergence criterion based on goodness of fit, we specified the optimization routine to produce the best parameter set (evaluated by the fitness criteria described below) found after 200 generations.

When the genetic algorithm converged on a parameter set, to test whether the parameter set could be improved within local variations, we used a simplex search to find the local error minimum (described in LeMasson and Maex 2001). We chose the Nelder–Mead simplex (direct search) method, implemented in MATLAB as the standard function "fmins." Beginning with the parameter set optimized by the genetic algorithm, the simplex function adjusted parameters and ran the GENESIS simulations with the adjusted parameter set to produce new output from the model. This activity was evaluated with the same fitness criteria, described below, as that used for the genetic algorithm parameter optimization, although the final fitness value was multiplied by –1 for the simplex search. This simplex parameter optimization is an iterative process of adjusting parameters and evaluating how the parameter adjustments improve the fitness (bring the fitness value close to zero). We used all default options, as implemented in MATLAB.

Choice of boundaries for conductance densities

The upper and lower bounds for each conductance density were determined based on the conductance densities chosen while previously tuning the Full Model to synaptically isolated tonic spiking activity (Tobin et al. 2006). As described in the companion paper (Tobin et al. 2006), initial conductance densities were chosen to match voltage-clamp amplitudes from published voltage-clamp data, from simulated voltage-clamp in the single-compartment model, or values that represented a compromise between the two. After extensive conductance tuning to match synaptically isolated activity, the resulting parameter set produced conductance amplitudes, as measured by simulated voltage clamp, that were within 14% of initially chosen values, except I_CaS, which was increased by 42%. Because these conductance densities were close to those initially chosen to match voltage-clamp data, except for I_CaS, we used these initial conductance densities as the "default" values to determine the genetic algorithm boundaries. However, because tuning to synaptically isolated activity increased I_CaS densities in the Neurite compartments, we chose for I_CaS a default value that was a compromise between the initial and tuned values. Although these conductance densities were tuned for the Full Model, when they were applied to the Reduced Model the resulting conductance amplitudes, as measured by simulated voltage clamp, did not appreciably change (amplitudes were reduced between 0 and 15%; data not shown). Thus we used these default values to set the boundaries from which the genetic algorithm could choose conductance density values for tuning the Reduced Model. The conductance boundaries for the genetic algorithm were chosen as between about one third of the default value and about threefold the default value (Supplemental Table 1). However, we decreased the lower bound of I_h and I_P in the Secondary Neurite to 0 to allow for no conductance there, in accordance with the Full Model tuned parameter set (Tobin et al. 2006). To constrain action potential height and initiation in the Neurite–axon2 compartment, and to limit parameter space, we fixed I_Na to the default value (3.7 mS/cm²), as tuned in the Full Model (Tobin et al. 2006).

Evaluation of the fitness of a parameter set was based on whether it could produce bursting activity, as defined below, and how closely the burst characteristics matched ideal values, as defined below. Low fitness values were assigned to those parameter sets that most closely matched the desired characteristics. For each parameter set, model activity was simulated with an initial 1-s, –0.15-nA somatic current injection followed by 79 s of run time. Soma voltage was written to an ASCII file and the portion between 20 and 80 s was analyzed as follows. Spike detection and burst discrimination were performed as described in Masino and Calabrese (2002a). In summary, spikes were detected when the voltage crossed a threshold of –30 mV. To prevent double counts of a single spike, a 10-ms refractory period was used (Masino and Calabrese 2002a). Bursts were defined as groups of three or more spikes with 700 ms between each group of spikes. Burst period was calculated as the time between the median spikes of consecutive bursts. Burst duration is defined as the time from the first to the last spike of each burst, and duty cycle is the ratio of burst duration and period for each burst. All parameter sets that did not produce bursting activity were classified as "tonic spiking" or "no spiking," and were assigned the highest (worst) fitness value of 400.

For each parameter set for which bursting was detected, the mean period, SD of period, mean duty cycle, and mean maximum intraburst spike frequency were measured. The criteria for assigning fitness values to parameter sets were developed through successive trials of genetic algorithm runs to find criteria that successfully specified desirable bursting activity. Several types of undesirable bursting activity, described below, were assigned intermediate fitness values of 100, so that they were more likely to contribute to successive generations than the nonbursting parameter sets, but less than those with more desirable bursting activity. Irregular bursting was defined as bursting for which the SD of bursting was >10% of the mean period and was assigned the intermediate fitness value of 100. Nonsustainable bursting was defined as bursting for which there was more than 1.5-fold the duration of the mean period between the last spike of the last burst and the end of the trace (seen as bursting that falls into steady-state membrane potential before the end of the trace). Nonsustainable bursting was assigned a fitness value of 100. If the duty cycle was >80 or <40%, or if the time between bursts was <2.0 s, the parameter set was assigned a fitness value of 100.

Regular bursting that did not fit the above criteria was evaluated based on how closely the period and spike frequency matched the desired values of 8.0 s for period and 13.0 Hz for maximum spike frequency (Supplemental Methods, "Evaluation Criteria for Genetic Algorithm and Simplex Search"). These desired values were chosen to approximate the measured burst characteristics of endogenous bursters. As measured by Cymbalyuk et al. (2002), the period of endogenously bursting heart interneurons ranged from 4.8 to 8.6 s (mean: 6.5 ± 1.0 s) and mean spike frequency ranged from 6.6 to 14.7 Hz (mean: 9.3 ± 2.2 Hz). Our measurements of endogenous bursting in 0.5 mM bicuculline indicated burst periods from 5.0 to 8.4 s, with a mean of 6.6 ± 1.2 s; mean spike frequencies from 5.1 to 11.2 Hz, with a mean of 7.2 ± 2.4 Hz; and maximum spike frequencies from 7.7 to 16.6 Hz, with a mean of 11.3 ± 3.6 Hz (Tobin and Calabrese 2005). This variability in the rhythmic output of the heart interneurons indicates that no one parameter set could describe the range of oscillator heart interneurons, and thus we are seeking multiple parameter sets rather than only one ideal set.

The fitness value was calculated as a sum of the "period fitness" (a weighted difference between the desired period and actual period) and the "spike frequency fitness" (a weighted difference between the desired spike frequency and model spike frequency). For periods <3.0 or >14.0 s, "period fitness" was calculated as 16-fold the absolute value of the difference between the desired period and the period; if the period was between 3.0 and 14.0 s, "period fitness" was eightfold the absolute value of the difference between the desired period and the period. This weighting system was designed to more heavily penalize bursting with periods outside the range normally measured in heart interneurons. Likewise, if the maximum spike frequency was <4.0 Hz (much lower than experimentally measured), an overall fitness value of 100 was assigned. If the maximum spike frequency was >20.0 or <8.0 Hz (outside experimentally measured ranges), "spike frequency fitness" was twice the absolute value of the difference between the desired maximum spike frequency and the maximum spike frequency; if the spike frequency was between 8.0 and 20.0 Hz, the "spike frequency fitness" was the absolute value of the difference between the desired maximum spike frequency and the maximum spike frequency (MATLAB programming code provided in Supplemental Methods¹ ). In addition to penalizing values of burst characteristics that were not within the experimentally observed range, this weighting system also weighed differences in period more heavily than differences of spike frequency, as it was found in trial runs that spike frequency was consistently more closely matched than period.

Sustained endogenous and half-center bursting

We defined sustained bursting as bursting that continued throughout 300 s of simulation time. After parameter sets were chosen for their ability to sustain bursting for 300 s, all further simulations discussed in this text were run to 300 s and burst characteristics were averaged from 240 to 300 s of simulation time. When models were configured as half-centers, one cell received an initial 1-s, –0.15-nA somatic current injection, identical to that in the single-cell models. This asymmetry in starting conditions can cause asymmetric bursting between the two model neurons of the half-center. Unless otherwise noted, activity of models in a half-center was assumed to be symmetric (burst characteristics of noninjected cell were within 10% of those of the injected cell) and the values for the injected cell were quoted for comparison to the single-cell simulations with identical simulated current injection.

Mimicking microelectrode penetration

We modeled the electrode-induced soma leak conductance in a way similar to that described in Cymbalyuk et al. (2002), by introducing into the soma compartment a non-voltage-sensitive leak conductance with reversal potential . The original soma leak conductance and reversal potential values were replaced by the new soma leak conductance g_Leak and reversal potential E_Leak, calculated as (Cymbalyuk et al. 2002; their Eq. 2)

We increased the electrode-induced soma leak conductance until the endogenous bursters spiked tonically, and thus defined a critical value for .

Modeling Na/K pump inhibition

In previous studies exploring the effect of the neuropeptide myomodulin on membrane properties, we measured an inhibition of the Na/K pump, indicated by a myomodulin-induced average shift of 0.058 nA in the voltage-clamp holding current during hyperpolarizing voltage ramps in 0 Ca²⁺, 1.8 mM Mn²⁺ saline with 2 mM Cs⁺ (average holding current in myomodulin: 0.0423 ± 0.007 nA; control: 0.10 ± 0.008 nA; mean ± SE) (Tobin and Calabrese 2005). To simulate the effect of inhibiting the Na/K pump, we injected a constant current ("pump current") into each compartment, scaled by the compartment membrane area. The amplitude of simulated pump current was set for each model by simulating hyperpolarizing voltage-clamp ramps without and with pump current and choosing pump current amplitudes to shift the holding current by 0.058 nA, as empirically observed during myomodulin application. The effect of this pump current on burst characteristics was evaluated by adding the pump current to all compartments during simulated current clamp, with an initial 1-s, –0.15-nA somatic current injection, as in all other current-clamp simulations. The typical mechanisms by which voltage clamp and current clamp are simulated in GENESIS overwrite any additional current injections into the clamped compartment. Thus all simulations using pump current do not include pump current in the soma. Because nonsomatic pump current must be increased to compensate for the missing current in the soma, pump current amplitudes are likely overestimated.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES

 
Parameterization of the Reduced Model using a genetic algorithm

The genetic algorithm converged on a parameter set called "Parameter Set 1" (Supplemental Table 2) that produced an endogenously bursting Single-Cell Reduced Model with a period of 4.6 ± 0.1 s, a mean spike frequency of 9.2 ± 0.5 Hz, and a maximum spike frequency of 12.9 ± 1.2 Hz (Table 1). When we configured the Reduced Model into a half-center oscillator with reciprocally inhibitory synapses, the two model cells burst in alternation, with a period of 6.5 ± 0.0 s, mean spike frequency of 17.4 ± 0.3 Hz, and maximum spike frequency of 26.5 ± 0.4 Hz. This result corresponds to an increase in period and spike frequency compared with that in the Single-Cell Model. In heart interneurons bursting in a half-center configuration, we measured periods of 6.1 to 13.2 s (mean 10.2 ± 2.4 s) and mean spike frequencies from 6.8 to 15.3 Hz (mean 10.6 ± 1.5 Hz) (Tobin and Calabrese 2005), and from this same data set we here report maximum spike frequencies of 11.3 to 22.1 Hz (mean 15.8 ± 3.4 Hz). These models' burst characteristics for endogenous bursting are within experimental observations and, for half-center bursting, the period is similar to experimental observations and the spike frequencies are slightly higher.

To test whether the Reduced Model and Full Model are sufficiently similar that the former can be exploited with parameter optimization methods to parameterize the latter, we applied Parameter Set 1 to the Full Model and found that the Full Model did indeed burst endogenously. The burst period in the Full Model, 3.2 ± 0.2 s, was much less than that in the Reduced Model, whereas the mean (10.4 ± 3.2 Hz) and maximum (13.6 ± 6.0 Hz) spike frequencies were greater. When the Full Model was configured into a half-center oscillator, the burst period increased to 6.5 ± 0.03 s, the mean spike frequency to 20.4 ± 0.3 Hz, and the maximum spike frequency to 32.1 ± 0.4 Hz. The differences in period and spike frequency between the Full and Reduced Models demonstrate that the models are not identical but apparently similar enough to share/exchange parameter sets and achieve qualitatively similar activity patterns: endogenous and alternating bursting.

Multiple-parameter sets produce endogenous bursting

Modeling studies have shown that different conductance sets can produce similar bursting activity in models of rhythmically active neurons (Prinz et al. 2003). Moreover, experimental studies indicate that neurons may use a variety of conductance sets to produce the same activity (Golowasch et al. 1999, 2002). Therefore we asked whether a different parameter set could produce endogenous bursting in our Reduced Model. We ran the parameter optimization routine seven more times and assessed the resulting model parameter sets (Supplemental Table 2). From the eight total genetic algorithm runs, the resulting parameter sets all produced endogenous bursting with periods between 4.1 and 6.1 s, mean spike frequencies between 9.2 and 14.0 Hz, and maximum spike frequencies between 12.9 and 19.7 Hz (Table 1). As measured by Cymbalyuk et al. (2002), periods of endogenously bursting heart interneurons ranged from 4.8 to 8.6 s (mean: 6.5 ± 1.0 s) and mean spike frequencies ranged from 6.6 to 14.7 Hz (mean: 9.3 ± 2.2 Hz). Our measurements of endogenous bursting in 0.5 mM bicuculline indicate burst periods from 5.0 to 8.4 (mean: 6.6 ± 1.2 s), mean spike frequencies from 5.1 to 11.2 Hz (mean: 7.2 ± 2.4 Hz), and maximum spike frequencies from 7.7 to 16.6 Hz (mean: 11.3 ± 3.6 Hz) (Tobin and Calabrese 2005). Although a few model parameter sets produced maximum spike frequencies higher than those observed in the living system, overall the burst characteristics produced by these parameter sets are biologically realistic.

All parameter sets that produced endogenous bursting in the Reduced Model produced alternating bursting when the Reduced Model was configured as a half-center oscillator. When applied to the Full Model, all parameter sets produced endogenous and half-center bursting. In both the Reduced and Full Models, with all parameter sets, the period, mean, and maximum spike frequencies of alternating bursting were higher than those of endogenous bursting, in accordance with experimental data (Tobin and Calabrese 2005). In all models, half-center coupling increased the period by an average of 67.4 ± 21.1% and mean spike frequency by 74.5 ± 61.3% (as calculated from activity 20 to 80 s after the start of the simulation).

Sustaining bursting

Multiple parameter sets produced endogenous and half-center bursting in our models. However, to enable the genetic algorithm to converge within a reasonable time, we had limited the model simulation time to 80 s. To single out which of these parameter sets could produce sustained bursting, we ran the Reduced and the Full Models in single and half-center configurations for 300 s with each parameter set. Of the Single-Cell Reduced Models, only those with Parameter Sets 2, 3, 4, and 8 maintained endogenous bursting (Fig. 1), whereas those with Parameter Sets 1, 5, 6, and 7 developed a nonspiking steady resting potential within 300 s. All Half-Center Reduced Models maintained bursting (for examples, see Fig. 1) as did all Single-Cell Full Models and Half-Center Full Models (for examples, see Fig. 2). Because Parameter Sets 2, 3, 4, and 8 were the only ones to maintain endogenous and half-center bursting in both the Reduced and Full Models, we chose these sets as the most promising for realistically representing bursting in oscillator heart interneurons. That only these four parameter sets sustained endogenous bursting in the Full and Reduced Models, whereas all models sustained half-center bursting, supports findings from a single-compartment model of heart interneurons that the half-center configuration stabilizes bursting (Cymbalyuk et al. 2002). In all subsequent studies here only models with these four "robust" parameter sets are considered.

View larger version (26K): [in this window] [in a new window]   FIG. 1. Soma membrane potential is recorded for the Reduced Model with the parameter sets that sustained bursting in all models: Reduced 2 (A), Reduced 3 (B), Reduced 4 (C), and Reduced 8 (D) Models. For A–D: 1, endogenous bursting in Single-Cell Models, shown from 280 to 300 s from start of simulation; 2, alternating bursting in Half-Center Models, shown from 280 to 300 s from start of simulation. For 1 and 2: period and mean spike frequency are calculated for activity from 240 to 300 s after start of simulation.

View larger version (26K): [in this window] [in a new window]   FIG. 2. Soma membrane potential is recorded for the Full Model with the parameter sets that sustained bursting in all models: Full 2 (A), Full 3 (B), Full 4 (C), and Full 8 (D) Models. For A–D: 1, endogenous bursting in Single-Cell Models, shown from 280 to 300 s from start of simulation; 2, alternating bursting in Half-Center Models, from 280 to 300 s from start of simulation. For 1 and 2: period and mean spike frequency are calculated for activity from 240 to 300 s after start of simulation.

Comparing the conductance densities for all parameter sets chosen by the genetic algorithm (including parameter sets where bursting was sustained as well as those where it was not sustained), some trends are noticeable (Fig. 3). For example, although most conductance densities varied throughout the allowable range, neuritic I_KA values were clustered within the top half of the allowable range, suggesting that low neuritic I_KA may be unfavorable for bursting with our desired characteristics. I_K1 densities tended toward a bimodal distribution in the primary and secondary neurites, most noticeably in the axon. Of the inward currents, I_P values clustered in the lower half of the allowable range in both the primary and secondary neurites and, in the secondary neurite, I_h was relatively low, whereas I_CaS clustered in the upper half of the allowable range. Because the four parameter sets that produce sustained bursting likely represent only a few of the potentially large variety of bursters, we can highlight trends that appear, but cannot determine whether these trends represent rules of which distributions can or cannot produce bursting. Whereas a variety of axonal I_K2 conductance densities were observed among all bursters, the sustained bursters clustered near the top of the range. Similarly whereas I_KA and I_K1 spanned the range of allowable values in the secondary neurite, sustained bursters clustered at low values. The average I_h of sustained bursters was lower than that of the nonsustained bursters, as was the membrane resistance.

View larger version (51K): [in this window] [in a new window]   FIG. 3. For all parameter sets chosen by the parameter optimization routine for producing endogenous bursting in the Reduced Model, conductance densities are plotted, with black crosses, by compartment type. ELeak and Rm, which had the same value across all compartments, are plotted under the heading "homogeneous." Those parameter sets that exhibited sustained bursting in all models (Single-Cell, Half-Center, Reduced, and Full) are shown as black crosses highlighted by colored circles; legend displays parameter set number. Gray bars indicate the boundaries used by the parameter optimization routine (see Supplemental Table 1).

View larger version (42K): [in this window] [in a new window]   FIG. 4. Conductance density distributions of parameter sets 2, 3, 4, and 8 are compared. Icon represents compartment types as labeled. For each conductance, densities are normalized to the maximum density occurring in any of the 8 parameter sets from the genetic algorithm (regardless of compartment type). Thus densities are scaled with respect to only those values that are associated with endogenous bursting to enable comparison of densities within a model and across models.

To determine whether the bursting observed in both the single-cell and half-center configurations of the Reduced and Full Models was sensitive to small parameter variations, which would suggest "small hypervolume" in parameter space pinpointed by the genetic algorithm, we did a sensitivity analysis (Calabrese et al. 2001). Sequentially, we varied by ±5% the values of R_m and E_Leak and the density of each of the seven conductances, in the compartments where they were present, and we measured whether the Reduced and Full Models with Parameter Sets 2, 3, 4, and 8 could still produce endogenous and alternating bursting for 300 s.

When the conductance of I_K1 was varied, endogenous bursting was sustained for ±5% changes in Single-Cell Reduced 2, 3, and 8 and was not sustained in Single-Cell Reduced 4. For all other conductances, the Single-Cell Reduced 2, 3, 4, and 8 models were unable to sustain endogenous bursting when outward currents (I_K2, I_KA) were increased, when I_Leak was increased (modeled as a decrease in R_m), or when E_Leak was hyperpolarized. Similarly, decreasing inward currents (I_h, I_P, I_CaF, and I_CaS) resulted in nonspiking traces. The Single-Cell Reduced 8 differed from this trend only in that it was able to sustain bursting for both ±5% changes in I_CaF. All Half-Center Reduced Models were able to sustain bursting for all variations except that when E_Leak was hyperpolarized, the Half-Center Reduced 2, 3, and 4 stopped spiking. When the conductance variations were applied to the Full Model, all variations elicited endogenous bursting except for changes in E_Leak, for which hyperpolarizing caused no spiking in Single-Cell Full 2, 3, and 8, and depolarizing caused tonic spiking in Single-Cell Full 8. However, all Half-Center Full Models sustained bursting for all variations.

In general, the Single-Cell Reduced Models could not sustain endogenous bursting for hyperpolarizing changes: when outward currents were increased, E_Leak hyperpolarized or inward currents decreased. Half-Center Reduced Models and Single-Cell Full Models sustained bursting for all variations except E_Leak. Half-center bursting in the Full Models was robust to all changes. In experiments, endogenous bursting in heart interneurons is not always regular and simulations with the single-compartment model indicate it is very sensitive to leak parameters, whereas half-center bursting, in experiments and simulations, is far more regular (Cymbalyuk et al. 2002). Both the Reduced and Full Models corroborate these experimental and model results. Furthermore, these results indicate that bursting is more robust in the Full Model than in the Reduced Model, even though these parameter sets were tuned to the Reduced Model.

Assessing changes in burst characteristics with parameter variations

To assess changes in burst characteristics with parameter variations, we varied parameters by ±20% (to accentuate the small changes seen in 5% variations) in both Reduced and Full Models and calculated the percentage changes in average period and mean spike frequency. Full data sets are presented in Supplemental Tables 3–6. Because these burst characteristics have some variability within a simulated run, only those changes that are larger than the percentage SD of the burst characteristics measured from a simulation with the original parameter set are considered relevant. Because this criterion is less stringent than a statistical analysis (ANOVA with post hoc testing) for each variation, some statistically insignificant changes may be considered here.

Increasing the amplitude of the parameter variation to 20% led to a cessation of bursting for some models and some parameters, especially increases in outward currents (I_K1, I_K2, I_KA) and decreases of burst-supporting inward currents (I_P, I_CaF, I_CaS) or I_h. Single-Cell Reduced 4, on the other hand, unable to burst with a 5% decrease in I_K1, exhibited bursting with a 20% decrease in I_K1 (Fig. 5 B2). Otherwise, the responses of the models to 20% variations followed patterns similar to those with 5% variations (Figs. 5–8). In describing the responses of the models to the parameter variations here, we limit ourselves to models where bursting was maintained and focus on the general trends, mentioning exceptions only where they appear pertinent.

View larger version (30K): [in this window] [in a new window]   FIG. 5. Changes in burst characteristics arising from 20% conductance variations are measured in Single-Cell Reduced Models. Each conductance density, including Rm and ELeak, was varied by +20% (A1 and A2) and –20% (B1 and B2) in all compartments where the conductance was present, and changes in period (A1 and B1) and spike frequency (A2 and B2) were measured. Open shapes correspond to changes that were deemed nonrelevant, defined as the changes being within the SD of burst characteristics from the original parameter set. Where conductance variations stopped bursting, activity is indicated as "No spiking" or "Tonic spiking."

View larger version (31K): [in this window] [in a new window]   FIG. 8. Changes in burst characteristics arising from 20% conductance variations are measured in Half-Center Full Models. Each conductance density, including Rm and ELeak, was varied by +20% (A1 and A2) and –20% (B1 and B2) in all compartments where the conductance was present, and changes in period (A1 and B1) and spike frequency (A2 and B2) were measured. Open shapes correspond to changes that were deemed nonrelevant, defined as the changes being within the SD of burst characteristics from the original parameter set. Where conductance variations stopped bursting, activity is indicated as "No spiking" or "Tonic spiking."

The Half-Center Reduced Models generally exhibited qualitatively similar changes to one another and to their corresponding Single-Cell Models for all parameter variations (Fig. 6). A notable exception was for variation in I_h, where now increasing I_h decreased period and decreasing I_h increased period according to expectation from the single-compartment model (Hill et al. 2001).

View larger version (31K): [in this window] [in a new window]   FIG. 6. Changes in burst characteristics arising from 20% conductance variations are measured in Half-Center Reduced Models. Each conductance density, including Rm and ELeak, was varied by +20% (A1 and A2) and –20% (B1 and B2) in all compartments where the conductance was present, and changes in period (A1 and B1) and spike frequency (A2 and B2) were measured. Open shapes correspond to changes that were deemed nonrelevant, defined as the changes being within the SD of burst characteristics from the original parameter set. Where conductance variations stopped bursting, activity is indicated as "No spiking" or "Tonic spiking."

View larger version (32K): [in this window] [in a new window]   FIG. 7. Changes in burst characteristics arising from 20% conductance variations are measured in Single-Cell Full Models. Each conductance density, including Rm and ELeak, was varied by +20% (A1 and A2) and –20% (B1 and B2) in all compartments where the conductance was present, and changes in period (A1 and B1) and spike frequency (A2 and B2) were measured. Open shapes correspond to changes that were deemed nonrelevant, defined as the changes being within the SD of burst characteristics from the original parameter set. Where conductance variations stopped bursting, activity is indicated as "No spiking" or "Tonic spiking."

The effect of microelectrode penetration in Reduced and Full Models: endogenous bursting versus tonic spiking

Synaptically isolated oscillator heart interneurons can burst endogenously when recorded extracellularly, but when impaled by a sharp microelectrode, they spike tonically (Cymbalyuk et al. 2002). We tested whether the Reduced and Full Models would reproduce this phenomenon. We introduced a nonspecific (E_Leak = 0 mV) "electrode-induced" leak conductance , calculated as in Cymbalyuk et al. (2002) (see METHODS), into the soma compartment of each neuron model and increased the electrode-induced leak until the endogenous bursters spiked tonically. All Reduced and Full Models had a critical value of this electrode-induced leak for which the endogenous bursting became tonic spiking [Reduced 3 example: Fig. 9A; for comparison to biological data, see Fig. 2B in the companion paper Tobin et al. (2006)]. Normally connected oscillator heart interneurons (i.e., in half-center configuration) burst in alternation, even when both cells are penetrated with microelectrodes. Likewise, all Half-Center Models, with either Reduced or Full morphology, also continued bursting in alternation with the critical value of the somatic electrode-induced leak conductance [Reduced 3 example: Fig. 9B; for comparison to biological data, see Fig. 2A in the companion paper Tobin et al. (2006)].

View larger version (27K): [in this window] [in a new window]   FIG. 9. Soma voltage is recorded for the Reduced 3 Model in the (A) Single-Cell and (B) Half-Center configurations while mimicking a microelectrode-induced soma leak.

The effects of blocking I_h with Cs⁺ in the Half-Center Reduced and Full Models

Only a few types of conductance variations are useful for comparing the burst regulation of the models to experimental data; for example, blocking Na⁺ currents stops bursting in interneurons. We chose to test the effect of Cs⁺, which blocks I_h in heart interneurons (Angstadt and Calabrese 1989). Adding 2 mM Cs⁺ to the saline increases period (Masino and Calabrese 2002b; Tobin and Calabrese 2005) and decreases spike frequency of extracellularly recorded heart interneurons (Tobin and Calabrese 2005). To test whether our Half-Center Reduced and Full Models could replicate these effects, we mimicked the biophysical effects of Cs⁺ by setting h conductances to 0. As indicated by the parameter variations, decreasing I_h by 20% increased period and spike frequency in most Half-Center Reduced Models (Fig. 6 and Supplemental Table 4) and Half-Center Full Models (Fig. 8 and Supplemental Table 6). Completely blocking I_h, however, stopped bursting and spiking in all Half-Center Reduced and Full Models. Either our models do not fully represent the response of heart interneurons to the Cs⁺-induced I_h block or there are other membrane properties that Cs⁺ affects, for which we have not accounted. For example, Cs⁺ may enhance the Na/K pump in leech neurons (Catarsi and Brunelli 1991; Skou 1960, 1965). Potential nonspecific effects of Cs⁺ are indicated by the nonmonotonic nature of the dose–response relation of Cs⁺ on period; period reduction peaks at 2 mM and increasing concentrations decrease the effect (Masino and Calabrese 2002b). Moreover,intracellularly recorded heart interneurons fire tonically in Cs⁺-containing saline (Angstadt and Calabrese 1989), indicating that leak and I_h interact in a way that the models might not yet capture. The drug ZD7288, shown to block I_h selectively in vertebrate preparations (Gasparini and DiFrancesco 1997), does not block I_h in leech heart interneurons (Tobin and Calabrese 2005), precluding a more specific block at present.

The effects of myomodulin in the Single-Cell and Half-Center Reduced and Full Models

Myomodulin decreases period and increases spike frequency when bath applied at 1 µM to normally connected and synaptically isolated heart interneurons (Tobin and Calabrese 2005). We demonstrated that myomodulin enhances the maximum conductance of I_h by 30% and likely inhibits the Na/K pump (Tobin and Calabrese 2005). In the Full and Reduced Single-Cell and Half-Center Models, we simulated the effects of myomodulin on membrane properties to test whether they could account for the observed changes in burst period and spike frequency of endogenous and half-center bursting in heart interneurons.

To mimic the 30% increase of I_h caused by myomodulin, we increased I_h conductance uniformly until we measured a 30% increase in simulated somatic voltage clamp (for Parameter Sets 2, 3, 4, and 8, I_h conductance densities were increased by 30, 35, 50, and 30%, respectively, in Reduced Models and by 38, 33, 48, and 30%, respectively, in the Full Models). The Na/K pump inhibition, measured experimentally as a inward current shift in the linear portion of voltage ramp currents (Tobin and Calabrese 2005), was mimicked by current injections of 45, 54, 50, and 40 A/cm² in the Reduced 2, 3, 4, and 8 Models, respectively, and by 46, 55, 52, and 41 A/cm² in the Full 2, 3, 4, and 8 Models, respectively. In experiments, values for the inhibited pump current were measured only in Cs⁺ saline in which the pump may already be enhanced by the Cs⁺ (Catarsi and Brunelli 1991; Skou 1960, 1965). Thus we used these values as only an estimate of the effect myomodulin may have in normal saline, and for Single-Cell Models where the full amplitude of pump current caused tonic spiking (Reduced 8, Full 3, Full 8), we use half the amplitude of the inhibited pump current for both the Single-Cell and Half-Center Models. For Reduced 3, half-amplitude pump current was used only when both I_h and pump current were combined. For each of the modulatory effects, separately and combined, we measured burst characteristics for the Single-Cell and Half-Center Full and Reduced Models (Fig. 10 and Supplemental Table 7).

View larger version (23K): [in this window] [in a new window]   FIG. 10. Effects of the neuromodulator myomodulin are mimicked in the models by increasing Ih conductance by 30% and inhibiting the Na/K pump (simulated by current injection; see RESULTS). Each effect was applied separately and combined to qualitatively assess the individual and combined effects of the biophysical changes on period and spike frequency for (A) Single-Cell Reduced Models, (B) Single-Cell Full Models, (C) Half-Center Reduced Models, and (D) Half-Center Full Models. Where the full amplitude of pump current caused tonic spiking in Single-Cell Models (Reduced 8, Full 3, and Full 8), we use half the amplitude of the pump current for both the Single-Cell and Half-Center Models. Half-amplitude pump current was used for Reduced 3 only when both Ih and pump current were combined. B: Single-Cell Full 8 exhibited tonic spiking even when the pump current was reduced to half-amplitude. D: *Half-Center Full 4 exhibited asymmetric bursting such that the spike frequency of one model neuron decreased by 27.1%, whereas the spike frequency of the other model neuron increased by 17.4%. Periods of each model neuron decreased by 17.4 and 17.5%, respectively. Heart interneuron response to myomodulin is plotted with open triangles with SD bars. Response of synaptically isolated (bicuculline) heart interneurons to myomodulin is plotted with Single-Cell Models under the heading "both." Response of synaptically intact heart interneurons to myomodulin is plotted with Half-Center Models: under the heading "Ipump" when Ih is blocked by Cs+ and under the heading "both" when Ih and Ipump are intact (normal saline).

In the Single-Cell Full Models (Fig. 10B), increasing I_h, adding pump current, and both effects combined decreased period and spike frequency (except for Full 8, where even half the pump current caused tonic spiking). In each case, the effects were sublinearly additive. Whereas the decrease in period matches that measured experimentally, the decrease in spike frequency does not (Tobin and Calabrese 2005) (Fig. 10B). In the Half-Center Full Models (Fig. 10D), increasing I_h, adding pump current, and both effects combined decreased period and increased spike frequency, as seen experimentally with myomodulin (Fig. 10B) (Tobin and Calabrese 2005) (except for Half-Center Full 2, where increasing I_h slightly decreased spike frequency, and Half-Center Full 4, where the I_h and pump effects combined increased spike frequency in one model neuron and decreased it in the other). These effects were sublinearly additive as well. In general, the Single-Cell and Half-Center Full Models indicate that each of the biophysical effects of myomodulin, separately and combined, should decrease spike frequency in endogenously bursting cells, increase spike frequency of normally connected cells, and decrease period in either configuration.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES

 
Using an automated parameter optimization routine, we recovered four parameter sets (conductance density distributions) that capture stable endogenous bursting and half-center bursting in our Reduced and corresponding Full Model of a leech oscillator heart interneuron. We then subjected the models with these parameter sets to parameter variations and tests corresponding to experimental measures: microelectrode-induced soma leak (Cymbalyuk et al. 2002), Cs⁺ application, and application of the neuromodulator myomodulin (Tobin and Calabrese 2005).

Choice of optimization algorithm and convergence criteria for model tuning

Our Full Model of a heart interneuron is too computationally intensive for effective automated parameter optimization, but the design of a corresponding Reduced Model with single compartments that map onto functional regions of the Full Model (Tobin et al. 2006) allowed us to use the Reduced Model as a proxy for the Full Model in parameter searches. Because of the high dimensionality of our Reduced Model (22 parameters), we chose a genetic algorithm for automated parameter optimization; this method has been shown to be among the most effective parameter search methods for large parameter spaces (Vanier and Bower 1999). We designed the fitness criteria to select for endogenous bursting with burst characteristics (period, spike frequency, and duty cycle) to match extracellularly recorded endogenous bursting in pharmacologically isolated heart interneurons (Cymbalyuk et al. 2002; Tobin and Calabrese 2005). We chose these characteristics as appropriate to the neurons/network studied. First, these measures are easily repeatable and represent native activity that is as unbiased as possible by the recording method. Second, because we sought to use these models in network simulations of the heartbeat pattern generator, period and spike frequency are important output parameters in the living system (Norris et al. 2006).

The specific target values chosen for the genetic algorithm were arrived at empirically. Trial runs indicated that constraining maximum spike frequency, rather than mean spike frequency, would more likely yield realistic intraburst spike frequency profiles. Such trial runs also indicated that spike frequency was more easily matched than period, leading us to weight period matching more heavily. Moreover, our experience with hand-tuning the single-compartment model on which the present morphological models are based (Hill et al. 2001) suggested that obtaining long periods was more difficult than short ones, so we picked a period target near the end of the observed range in living pharmacologically isolated extracellularly recorded heart interneurons (4.8–8.6 s; Cymbalyuk at al. 2002; Tobin and Calabrese 2005). Duty cycle was only loosely constrained, with equal fitness assigned for a broad range of duty cycle 40–80% because we reasoned that in half-center bursting strong mutual inhibition would constrain duty cycle to the nearly 50% observed in the living system. Because we used three criteria (period, maximum spike frequency, and duty cycle) for assigning fitness values, we did not expect to match all values exactly and we assigned weights for period and spike frequency in a two-tiered fashion, such that those values outside the biological range were penalized more heavily in their deviation from the desired value than those within the range of biological variability. Not unexpectedly we obtained model activity that did not closely match any of the target values, but exhibited biologically realistic activity. During the selection process, bursting activity was assessed over 80 s of simulation time for computational efficiency. This decision necessitated a further assay of stability; each of the eight sets independently chosen by the genetic algorithm was run for 300 s of simulation time. Only four sets produced sustained bursting and were thus chosen for further analysis. For these four sets the period was uniformly below and the maximum spike frequency uniformly above the target values and duty cycle was uniformly at the low end of the target range (Table 1). Judging from this limited number of searches, it appears that the searches converge on a competitive compromise of the target values.

Our experience suggests that in parameterizing a neuronal model, such optimization criteria and weighting should be chosen with regard to the system of interest and the desired goal. Characteristics other than those chosen here, such as amplitude of slow waveform or number of spikes per burst, in rhythmically active neurons, or spike timing, spike shape, or response to stimulation in input driven neurons may be considered important characteristics defining the activity and function of a neuron. Additionally, where possible, recordings from multiple sites in a neuron may be particularly useful in constraining multicompartmental models with nonhomogeneous conductance distributions (Keren et al. 2005).

Use of a Reduced Model for parameter optimization of a Full Model

The Reduced Model of an oscillator heart interneuron was created as a computationally efficient representation of the general morphology and the passive and active properties of the detailed morphology Full Model. The geometry of the Reduced Model was constrained to match both passive and active properties of the Full Model. Because these properties are determined by the conductance densities and distributions of the Reduced and Full Model, the ge