Abstract
Most neurons have large numbers of voltage and timedependent currents that contribute to their electrical firing patterns. Because these currents are nonlinear, it can be difficult to determine the role each current plays in determining how a neuron fires. The lateral pyloric (LP) neuron of the stomatogastric ganglion of decapod crustaceans has been studied extensively biophysically. We constructed ∼600,000 versions of a fourcompartment model of the LP neuron and distributed 11 different currents into the compartments. From these, we selected ∼1300 models that match well the electrophysiological properties of the biological neuron. Interestingly, correlations that were seen in the expression of channel mRNA in biological studies were not found across the ∼1300 admissible LP neuron models, suggesting that the electrical phenotype does not require these correlations. We used cubic fits of the function from maximal conductances to a series of electrophysiological properties to ask which conductances predominantly influence input conductance, resting membrane potential, resting spike rate, phasing of activity in response to rhythmic inhibition, and several other properties. In all cases, multiple conductances contribute to the measured property, and the combinations of currents that strongly influence each property differ. These methods can be used to understand how multiple currents in any candidate neuron interact to determine the cell's electrophysiological behavior.
Introduction
The mixture of voltage and calciumdependent channels and their distribution across the surface of a neuron collectively determine the electrophysiological properties of the neuron, including its resting membrane potential, input impedance, firing rate in response to current, and the extent to which it exhibits neuronal behaviors such as bursting, plateauing, and postinhibitory rebound (Marder, 1998). For a neuronal network to function properly, the neurons in that network must exhibit networkappropriate behaviors.
Recent work has shown that although a neuron may behave similarly from animal to animal, the densities of voltagegated channels can vary severalfold from animal to animal (Golowasch and Marder, 1992a; Golowasch et al., 2002; MacLean et al., 2005; Schulz et al., 2006, 2007). How, then, is reliable neuronal behavior maintained in the face of this variability? One possible answer is that although conductances are variable, they are variable in a correlated manner, and these correlations preserve the function of the neuron. This idea was given additional plausibility by the discovery that many channels are expressed in a correlated manner, at least at the level of mRNA (Schulz et al., 2006, 2007).
We sought to examine the hypothesis that correlations between conductances are required to conserve the electrophysiological function of neurons. We generated a large population of models of a single identified neuron, the lateral pyloric (LP) neuron, found in the stomatogastric ganglion (STG) of the crab Cancer borealis. We generated this population by choosing model parameters at random and filtering out models that did not behave appropriately. We then measured a number of electrophysiological properties in the resulting models and looked for correlations between model parameters that were induced by our constraints on electrophysiological properties.
Historically, a number of different approaches have been used to assess the role of a given channel in a given electrophysiological property of a neuron. One such approach is pharmacological blockade (Hille, 2001). Although this can be done successfully, it presupposes both that a given drug's cellular targets are known and that it is selective. Often, neither of these conditions is met. More recently, many have studied the effects of genetically altering or deleting currents on neuronal dynamics (Raman et al., 1997; Khaliq et al., 2003; MacLean et al., 2003; Swensen and Bean, 2005). Again, although these manipulations are often informative, their interpretation is complicated by compensatory mechanisms that occur during development. A third method is to use the dynamic clamp to add or subtract a current from a neuron (Prinz et al., 2004a).
The most common computational technique used to assess the role of a given parameter on the behavior of a neuron is sensitivity analysis (Butera et al., 1999; Hill et al., 2001, 2002; Jezzini et al., 2004). Here, we used a different approach to understand the extent to which combinations of neuronal conductances influence a number of different neuronal behaviors and compare this approach to sensitivity analysis.
Materials and Methods
Electrophysiology
Experiments were performed essentially as described previously (Goaillard et al., 2004). All preparations had two intact superior esophageal nerves and included the esophageal ganglion, but they had from zero to two intact inferior esophageal nerves.
Modeling voltagegated currents
The LP neuron models contained 10 voltage and calciumgated currents. These were the delayedrectifier potassium current (I_{Kd}); the transient potassium current (I_{A}); the calcium current (I_{Ca}); the calciumactivated potassium current (I_{KCa}); the hyperpolarizationactivated inward current (I_{h}); the proctolinactivated current (I_{pr}), a mixed inward current that is activated by several neuromodulators, including proctolin (Swensen and Marder, 2000); the fast sodium current (I_{Na}); the axonal delayedrectifier potassium current (I_{Kda}); and the axonal transient potassium current (I_{Aa}). The I_{A} current was modeled as the sum of a fast and a slow component (I_{Af} and I_{As}), but the maximal conductances of these currents were held in a fixed ratio (0.885:1), because they model a single voltagegated current as characterized physiologically. The model contained two leak currents, I_{leak} (the somatoneuritic leak) and I_{leak,ax} (the axonal leak conductance), both described by the usual formalism but with different maximal conductances.
The mathematical description of most of these currents (I_{Kd}, I_{A}, I_{Ca}, I_{KCa}, I_{h}, and I_{pr}) was based on descriptions in previous work (Buchholtz et al., 1992; Golowasch and Marder, 1992b), but with several revisions. In all cases, we tuned the halfactivation voltages and kinetics of these currents to provide good fits to voltageclamp data (Golowasch and Marder, 1992a; Turrigiano et al., 1995; Schulz et al., 2006). The axonal currents (I_{Na}, I_{Kda}, and I_{Aa}) were based on those in the ConnorStevens model of crab walking leg axon (Connor et al., 1977), but the halfactivation and kinetic parameters of these currents were tuned to have the desired currentclamp behavior in the baseline model (see below, Construction of baseline model). Finally, there were three synaptic currents in the final LP neuron models, representing synaptic inputs from other members of the pyloric network: I_{synAB}, I_{synPD}, and I_{synPY}. They are described below.
All model currents except I_{Ca} were modeled as ohmic currents, described by the following equation: where ḡ is the maximal conductance, w is the product of the gating variables raised to the appropriate powers, v is the membrane potential, and E is the reversal potential. The reversal potential of each current in the model depended on the permeant ions. The reversal potential was +55 mV for sodium currents, −80 mV for potassium currents (including I_{synPD}), −25 mV for I_{h}, −10 mV for I_{pr}, and −70 mV for I_{synAB} and I_{synPY}.
The calcium current was modeled using the GoldmanHodgkinKatz equation, and was given by the following: where P̄_{Ca} is the maximal calcium permeability, m_{Ca} is the calcium channel activation, h_{Ca} is the calcium channel inactivation, N_{A} is Avogadro's number (6.022 × 10^{23} mol^{−1}), q_{e} is the elementary charge (1.602 × 10^{−7} pC), z_{Ca} is the valence of a calcium ion (+2), [Ca^{2+}]_{i} is the internal calcium concentration, [Ca^{2+}]_{o} is the external calcium concentration (13 mm), k is Boltzmann's constant (1.381 × 10^{−23} J/K), and T is the temperature in Kelvin (283.15°K = 10°C).
In the model, we assumed that I_{Ca} and I_{KCa} channels were clustered together, and each cluster was associated with a calcium microdomain (Marrion and Tavalin, 1998; Fakler and Adelman, 2008). We know of no evidence either for or against this hypothesis in the STG, but it seems to be common in other systems (Fakler and Adelman, 2008). The maximal permeability of I_{Ca} was related to the percluster maximum permeability by the following: and the maximal conductance of I_{KCa} was related to the percluster maximal conductance by the following: where ρ_{clust} is the density of clusters (in clusters/nF), P̄_{Ca}^{1} is the percluster maximal permeability of I_{Ca}, and ḡ_{KCa}^{1} is the percluster maximal conductance of I_{KCa}. We simulated the accumulation of Ca^{2+} in an (average) small volume associated with each microdomain. The equation governing intracellular calcium was as follows: where [Ca^{2+}]_{∞} is the steadystate calcium concentration in the absence of calcium current (20 μm), τ_{Ca} is the time constant of calcium buffering (70.4 ms), and V_{Ca} is the volume associated with each microdomain (6.49 μm^{3}).
In generating our population of models, we picked values of P̄_{Ca} and ḡ_{KCa} at random and solved Equations 5 and 6 for ρ_{clust} and ḡ_{Ca}^{1}, assuming a fixed P̄_{Ca}^{1} of 1.1675 × 10^{−3} μm^{3}/ms. One consequence of this scheme is that an increase in P̄_{Ca} alone will not change the dynamics of calcium accumulation (because the calcium accumulation is done on a permicrodomain basis, and increasing P̄_{Ca} increases the density of microdomains but not the dynamics within a microdomain) and therefore will not change the dynamics of I_{KCa}. Another approach would be to make ρ_{clust}, P̄_{Ca}^{1}, and ḡ_{Ca}^{1} all be independently tunable parameters, but this seemed overly complex given the experimental uncertainties about the mechanics of calcium accumulation in this system.
Equations for the steadystate (in)activation and time constants of each voltage and calciumgated current are given in Table 1.
Modeling synaptic currents
During the pyloric rhythm, the LP neuron receives inhibition from three main sources: the anterior burster (AB) neuron, the pyloric dilator (PD) neurons, and the pyloric (PY) neurons. In the model, we approximated these inputs by fixed rhythmic patterns of synaptic conductance.
The synapses from AB and PY to LP are glutamatergic, predominately chloridemediated, and rapid (Hartline and Gassie, 1979; Eisen and Marder, 1982; Marder, 1984). The synapse from PD to LP, in contrast, is cholinergic, potassiummediated, and slow (Eisen and Marder, 1982; Marder and Eisen, 1984). To reflect this, the model AB and PY synapses had instantaneous dynamics and a reversal potential of −70 mV, and the PD synapse had slower dynamics (with a fixed time constant of 50 ms) and a reversal potential of −80 mV (the K^{+} reversal potential in the model). Each synapse was driven by a fixed presynaptic voltage waveform with a typical pyloric period of 1 s. These were based on smoothed recordings of AB (for the AB and PD synapses) and PY (for the PY synapse), which were then subjected to Fourier analysis. The presynaptic voltage used for the AB and PD synapses was as follows: where v_{AB/PD} is in millivolts, t is in seconds, and f_{0} = 1 Hz, the fundamental frequency of the simulated rhythm. Similarly, the presynaptic voltage used for the PY synapse was as follows: For all synapses, the synaptic conductance was then given by the following equations: where g_{syn} is the synaptic conductance, ḡ_{syn} is the maximal synaptic conductance, s is the fractional activation of the synapse, τ_{syn} is the synaptic time constant, v_{pre} is the presynaptic voltage, v_{thresh} is the synaptic threshold, and v_{scale} sets the voltage sensitivity of the synapse. s_{∞}(x) is the NakaRushton function, given by 1/(1 + x^{−2}) for x > 0, and 0 otherwise. For the AB and PD synapses, v_{thresh} = −58 mV and v_{scale} = 15 mV. For the PY synapse, v_{thresh} = −53 mV and v_{scale} = 3 mV. For the AB and PY synapses, τ_{syn} = 0 ms, and for the PD synapse, τ_{syn} = 50 ms.
Construction of baseline model
Once we had established reasonable descriptions of the intrinsic and synaptic currents in the LP neuron, we next wanted to establish a morphology for the model that mimicked the morphology of the biological LP and captured some rudimentary aspects of its electrophysiology. Additionally, we feared that a “brute force” search of possible morphologies and possible conductance distributions would be prohibitively time consuming. We therefore built a “baseline” model that was similar to LP in a number of ways: it had a morphology grossly consistent with the biological LP, it had spikes that were generated distally and were greatly attenuated in the soma, its input impedance was similar to that of LP (data not shown), its spikefrequency versus voltage curve was similar to that of LP, and it had voltagegated currents similar to those of LP.
In developing the baseline model, we started with the simplest possible model and added complexity as needed to agree with the data. We eventually settled on a fourcompartment model (Fig. 1) that combined a good fit to the impedance over a range of frequencies (0.1–200 Hz), a reasonable spike height as recorded in the soma (∼9 mV), and a reasonable resting membrane potential (approximately −40 mV).
We then added the reformulated voltagegated conductances to the soma and neurite compartments of the fourcompartment model (Fig. 1). Somatoneuritic currents (I_{leak}, I_{Kd}, I_{A,} I_{Ca}, I_{KCa,} I_{h}, and I_{pr}) were present at uniform density in the soma and neurite compartments and were absent elsewhere. The synaptic currents were present at uniform density only in the neurite compartments, and the axonal currents (I_{leak,ax}, I_{Na}, I_{Kda}, and I_{Aa}) were present only in the axon compartment. All of these distributions are simplified versions of the true channel distributions in biological LP neurons (Hartline and Graubard, 1992; Baro and HarrisWarrick, 1998; Baro et al., 2000; French et al., 2004), but they are likely to be reasonable approximations at the level of detail we are interested in at present. The maximal conductances of most voltagegated conductances were set to values close to those seen in the voltageclamp recordings. The single exception to this was ḡ_{KCa}, which was set to a value ∼40% of the value seen in voltageclamp recordings. We found that increasing ḡ_{KCa} beyond this value led to endogenous bursting, which is not usually observed in the isolated LP in C. borealis. (It is not entirely clear why this is the case, but we suspect that the published voltageclamp recordings were traces chosen for their large currents and that single LP neurons may not have all currents simultaneously large.)
Because the resulting model had physiological behavior that approximately matched that of the biological LP, with grossly reasonable maximal conductances, we declared this model an acceptable baseline model for later sampling of the parameter space. We note that the baseline model is not itself an admissible LP model. In particular, it bursts too late in the cycle in response to pyloriclike input, and the spike frequency during the burst is too low. This necessitated the use of sampling to generate admissible LP models (see Results, Production of LP model population).
All simulations of the baseline and sampled models were performed using NEURON (Carnevale and Hines, 2005) as a simulation engine, which was called from Matlab (Mathworks). Unless noted otherwise, the CVODE integration method was used with an absolute tolerance of 25 μV for voltage, 10 nm for internal [Ca^{2+}], and 0.001 for gating variables.
“Measuring” electrophysiological properties of models
We measured the electrophysiological properties of models under three different conditions: a small voltageclamp step, to assess input conductance; no input, to assess spontaneous activity; and a pattern of rhythmic synaptic inhibition, to mimic the input the LP neuron receives during the ongoing pyloric rhythm. These inputs were delivered to each model cell, and the output of the cell was analyzed to measure the various electrophysiological properties. All of the codes for managing the simulations, checking whether their responses had come to steady state, and characterizing the resulting model behaviors were written in the Matlab scripting language.
First condition: small voltageclamp step.
The cell soma was clamped to −60 mV and allowed to come to steady state. The soma potential was then stepped to −65 mV for 150 ms, a typical step duration used in experimental protocols to measure input conductance. The difference in steadystate current, ΔI, between the holding potential and the test potential was then divided by the potential difference, Δv = −5 mV, to obtain the input conductance, G_{in}. For this stimulus, the sodium conductance, proctolinactivated conductance, and synaptic conductances were all zeroed, to simulate the conditions in which these measurements were made experimentally. Additionally, an electrodeinduced shunt (with a conductance of 15 nS and a reversal potential of −35 mV) was included in the models for this protocol only, because an electrodeinduced shunt in the experiments will clearly have a strong effect on measurements of G_{in}. Because simulation of voltage clamp necessitates the solution of a stiff set of differential equations (Carnevale and Hines, 2005), this protocol was run using the backward Euler integration method, and a time step of 25 μs.
Second condition: no current injection.
This condition simulates the spontaneous activity of a neuron that has been isolated from pyloric inputs but still receives descending modulatory inputs. Therefore, the synaptic conductances were zeroed, but the proctolinactivated conductance and the sodium conductance were not. We simulated this condition until the model came to apparent steady state and computed a number of features based on its steadystate waveform. These include the qualitative type of activity (silent, periodic spiker, aperiodic spiker, periodic nonspiker, or aperiodic nonspiker); the minimum, maximum, mean, and SD of the soma potential, axon potential, and somatic Ca^{2+} concentration; the number of spikes; the spike frequency; the mean spike duration; the mean and SD of the interspike intervals; the duration of the burst; the minimum, maximum, mean, and SD of the slowwave somatic potential (the somatic potential if spikes are deleted, interpolated over, and the resulting waveform is lowpass filtered); the maximum spike height; and the cycle period. Some of these are only applicable for some qualitative activity types.
Third condition: periodic synaptic inhibition.
Each model was subjected to a pattern of rhythmic synaptic inhibition that approximates what the LP neuron receives during the ongoing pyloric rhythm (see above, Modeling synaptic currents). In this condition, no conductances were zeroed. The synaptic currents were injected into the model until it reached apparent steady state, and the same scalar properties described above were computed based on the steadystate waveform in this condition.
In more detail, the synaptic inputs were provided to each model for several cycles, using the same initial conditions as used for the spontaneous activity. After each cycle, the last cycle was compared with the penultimate cycle to determine whether the model had come to steady state. To be considered at steady state, four quantities had to be approximately equal on the last and the penultimate cycle: mean somatic voltage, SD of somatic voltage, mean somatic [Ca^{2+}], and SD of somatic [Ca^{2+}]. Voltages (both mean and SD) were considered “approximately equal” if they differed by <0.1 mV, and Ca^{2+} concentrations were considered approximately equal if they differed by <0.01 μm.
Statistics
To test whether two parameters were significantly correlated in the population of 1304 LP models, we calculated the correlation coefficient (r) between the two parameters and used a permutation test to determine whether the calculated r differed significantly from zero. A significance level of α = 0.001 was used for each pair of parameters to facilitate comparison with Schulz et al. (2007).
Cubic fits
To develop a simple description of how a particular electrophysiological property depends on the model parameters, we fit cubic functions to a set of “data” consisting of a population of LP models. Each model had a set of parameters (maximal conductances, etc.) and a number of electrophysiological properties (resting membrane potential, slowwave amplitude in response to rhythmic inhibition, etc.). One property at a time was fit, and each was treated as a scalar function of the parameters. To put all parameters on a common scale, each parameter was “zscored” (the mean value of the parameter in the LP model population was subtracted off, and the SD was then divided out) before fits were performed.
A quadratic polynomial in 17 variables has 17 + 2 choose 17 = 171 terms, each with an associated coefficient, and a full cubic model has 17 + 3 choose 17 = 1140 terms/coefficients. A polynomial with a small number of coefficients will generally have good generalization performance but may not be able to fit the target function well. A polynomial with a large number of coefficients will be better able to fit the target but may suffer from overfitting and therefore not generalize well (Hastie et al., 2001). To have finergrained control over the number of coefficients (and thus the fitquality/generalization tradeoff), we used a subset of the 1140 possible terms, selected using a crossvalidation procedure.
The LP models were first divided into a training set (90%) and a test set (10%). We used fivefold crossvalidation on the training set to determine a set of terms that were likely to have an optimal fit on novel models. We varied the number of terms, starting with a single constant term, and added terms until we found a set of terms that gave the lowest crossvalidation error. The next term to add was determined by temporarily adding all terms that could be obtained by increasing the degree of an existing term by 1 (not going beyond cubic terms), fitting the resulting polynomial to all of the training data, and retaining the most significant novel term (as determined by a t test). For this type of procedure, the crossvalidation error at first decreases as the number of terms increases, but after a critical number of terms are added, the crossvalidation error starts to increase. We then selected the simplest polynomial that had error not significantly worse than the minimum error, as assessed by a t test, and used this as our final polynomial. This was then fit to all of the training set to establish a final set of coefficients, and this fit was evaluated on the test set as a final metric of the quality of the fit. Determining the cubic coefficients involves finding the leastsquares solution of a set of equations that is linear in the coefficients. In all cases, the final number of coefficients was much less than the number of data points in the training set (largest number of coefficients was 158, number of training data points was 1174) (see Table 4), so the bestfit cubic coefficients were uniquely determined. SEs for the coefficients were calculated using the usual normaltheory formulas for linear regression (Draper and Smith, 1998).
Results
The LP neuron has a complex morphology (Fig. 2A), and it plays a critical role in the pyloric rhythm of decapod crustaceans (HarrisWarrick et al., 1992). The pyloric rhythm consists of three distinct phases: (1) the AB and PD neurons fire a burst of spikes (because of their strong electrical coupling, these neurons are collectively referred to as the “pacemaker kernel”); (2) the LP neuron fires a burst of spikes; and (3) the PY neurons fire a burst of spikes (Fig. 2B). The phase at which the LP neuron fires is remarkably consistent from animal to animal, despite substantial variation in the pyloric frequency across animals (Bucher et al., 2005) (J. M. Goaillard, A. L. Taylor, D. J. Schulz, and E. Marder, unpublished observation). In the crab C. borealis, the waveform of the LP neuron reflects the three phases of the pyloric rhythm: during the first phase, when the pacemaker kernel is active (Fig. 2B, left gray bar), it is deeply hyperpolarized; in the second phase, it fires a burst of spikes; and in the third phase, it is inhibited by the PY neurons, but not as much as during the pacemaker activity (Fig. 2B, right gray bar). The bursting of LP is primarily network driven: when it is synaptically isolated from pyloric inputs, LP fires tonically in C. borealis (Fig. 2C) (Golowasch and Marder, 1992a).
Previous work has shown that mRNA levels for a number of voltagegated channels in LP are correlated (Schulz et al., 2007). Furthermore, for two of three channel types examined, these subunit mRNA levels correlate with the conductance of the channel (Schulz et al., 2006). To gain insight into whether these sorts of correlations might arise from electrophysiological constraints on LP, we constructed a population of multicompartment models of LP, in which maximal conductances were varied and their effects on electrophysiological properties were examined.
Production of LP model population
To generate the LP model population, we first developed a baseline model (see Materials and Methods) and then generated ∼6 × 10^{5} models with randomized parameters in a large neighborhood around this baseline model (Fig. 3). For maximal conductances and maximal permeabilities, the range of possible values for each was from zero to approximately twice the value in the baseline model, so that each would span a realistic range. [Maximal conductances are found to vary threefold to fourfold in this system (Golowasch and Marder, 1992a; Goldman et al., 2001; Golowasch et al., 2002; Schulz et al., 2006, 2007).] Leak reversal potentials were varied over a 10 mV range centered on their baseline values, somewhat arbitrarily. The halfactivation potential of the proctolinactivated conductance (v_{1/2,pr}) was varied over a 20 mV range centered at −45 mV, approximately consistent with the range observed experimentally (Golowasch and Marder, 1992b; Swensen and Marder, 2000) (Goaillard, Taylor, Schulz, and Marder, unpublished observation). The parameters and their ranges are shown in Table 2. Each parameter was chosen from a uniform distribution within the given range, and each parameter was chosen independently of all other parameters. The capacitances and axial resistances of the model were fixed. This approach is similar to that taken in previous work (Foster et al., 1993; Goldman et al., 2001; Prinz et al., 2003, 2004b; Achard and De Schutter, 2006; Tobin et al., 2006; Hobbs and Hooper, 2008; Weaver and Wearne, 2008).
We chose 594,910 parameter sets in this way and performed simulations on each of the resulting models to measure their electrophysiological properties (see Materials and Methods). Each model was assessed in three different conditions. First, the model was subjected to a small voltageclamp step. Second, the model was given no input, and its spontaneous behavior was observed. Third, a periodic synaptic input was delivered to mimic the pyloric inputs that LP receives during the ongoing pyloric rhythm. In the voltageclamp condition, the model's input conductance was measured. In the noinput condition, several properties were measured, including the model's spontaneous activity pattern, its resting membrane potential, and its resting spike rate. In the pyloricinput condition, measurements included the phase of burst onset, the phase of burst offset, the spike rate in the burst, the slowwave amplitude, the peak slowwave membrane potential, and the coefficient of variation of the interspike interval in the burst. All properties were automatically measured by the simulation code (see Materials and Methods).
We then examined this population of ∼6 × 10^{5} models and selected only those with properties that were in a range that agreed with those of the biological LP neuron, using the bounds shown in Table 3. To be an acceptable LP model, the model had to (1) spike periodically when deprived of input; (2) have an interspike interval coefficient of variation <0.01 when deprived of input, to ensure that it spiked approximately tonically; (3) fire in a reliable pattern from cycle to cycle in response to the pyloriclike input; and (4) have all properties listed in Table 3 within the bounds given. In most cases, the bounds were chosen to contain the central ∼85% of the experimental data points (see Fig. 4, dashed lines). Resting spike rate was higher than those seen in experiments in which the LP neuron is isolated from pyloric network inputs, but similar to the spike rates within an LP burst. Input conductance was restricted to a range generally lower than previously reported experimentally [because the impedance measurements (data not shown) indicated a DC impedance generally larger than those reported previously], but the ratio between highest and lowest allowed values was similar to the ratio between the upper and lower bound of the central 85% of the experimentally measured values (Golowasch and Marder, 1992a).
Filtering the population of ∼6 × 10^{5} models using these criteria yielded 1304 acceptable models of LP. Traces of model LP activity look reasonable, as shown in Figure 3. Both the spontaneous activity of these models (compare Figs. 2C, 3A) and their activity in response to pyloriclike input (compare Figs. 2B, 3B) look similar to recordings of biological LP neurons. In response to pyloriclike inputs, the model LPs generally displayed a burst phase, a phase in which they were mildly hyperpolarized (during which LP receives PY inhibition), and a phase of deep hyperpolarization (during which LP receives pacemaker inhibition). These same three phases are present in recordings from biological LP neurons (Fig. 2B). Although strongly constrained electrophysiologically, the parameters of these models displayed wide variability (Fig. 3C) (also see below).
Histograms of the constrained properties are generally similar between model and experiment (Fig. 4). Consistent with this, the resting membrane potential (data not shown) was similar to that of the biological LPs (experiment: −46.8 ± 6.5 mV, n = 5; model: −44.1 ± 2.1 mV, n = 1304; mean ± SD). Two additional properties that we did not explicitly constrain at this stage were also consistent with the biology. The height of spikes as measured in the soma was similar between model and experiment (experiment: 7.7 ± 1.8 mV, n = 5; model: 8.0 ± 0.6 mV, n = 1304; mean ± SD). Also, the model impedances all matched the qualitative shape of the experimental impedance from 0.1 to 200 Hz (data not shown).
Weak correlations between pairs of parameters in the population of LP models
We examined whether the constraints imposed on the electrophysiological behavior of model LP neurons was reflected in the distributions of model parameters. Histograms of the model parameters for the 1304 models showed wide variability in most parameters (Fig. 5A). For most parameters, models could be found that spanned the sampled range. That is, when the range was binned into 10 bins, at least one model was found in all bins. Other parameters had restricted ranges in the LP population. The most strongly restricted parameter was ḡ_{Aa}, which was restricted to about half the sampled range. ḡ_{Kda} and ḡ_{synPY} also had restricted ranges.
Does the population of 1304 model LP neurons exhibit any correlations between pairs of model parameters? This was of particular interest because previous experimental results demonstrated correlations between different voltagegated conductances, some biophysically and some at the level of channel mRNA values (MacLean et al., 2005; Schulz et al., 2006, 2007). Therefore, we calculated correlation coefficients for each pair of parameters and tested these for significance (Fig. 5B) (see Materials and Methods). We found that correlations between different parameters in the LP population were generally weak (Fig. 5B) (red and blue backgrounds show significant correlations at α = 0.001 for individual comparisons). Although a rather stringent significance level of α = 0.001 was used for individual comparisons [and to facilitate comparison to Schulz et al. (2007)], this corresponds to a familywise significance level of α = 0.13 (n = 136 comparisons, assuming comparisons are independent), which is rather generous. Of the 136 pairs, 126 (93%) had estimated correlation coefficients between −0.2 and +0.2. The strongest correlation was between ḡ_{Na} and ḡ_{Kda}, which had a correlation coefficient (r) of 0.49 ± 0.04 [95% confidence interval (CI)].
To determine whether the observed correlations would likely be detected using typical experimental techniques, we picked a random subset of 20 models and tested for significant correlations (this time using a t test, rather than the more computationally expensive permutation test). n = 20 was chosen because it is a reasonable value of n for these types of experiments. This was repeated 2000 times, each time using a different random subset. A significant correlation (at α = 0.001) between ḡ_{Na} and ḡ_{Kda} was found in 12.5 ± 1.5% (95% CI) of trials, and all other pairs were lower than this. Of the 136 pairs examined, 131 (96%) had a chance of discovery of <1%. Together, these results imply that the constraints imposed on the electrophysiological properties of the model LP neurons were not sufficient to induce the strong pairwise correlations observed experimentally. This implies that, at least in the model, strong correlations like those measured experimentally are not required for a population to have the electrophysiological properties of an LP neuron.
The LP models are part of a single “island” in parameter space
We were curious to know whether all of the LP models were part of a single connected island of LP cells in parameter space. That is, is it possible to continuously vary the parameters of one LP model and thereby transform it into any other LP model, and to do so in such a way that all intermediate models are also LP models? We cannot prove mathematically that this is the case, but we have strong circumstantial evidence that all of the LP models are part of a single connected volume in parameter space.
To establish this, we drew a line segment between pairs of models in parameter space and tested four equally spaced points on this line to see whether they also satisfied the criteria for an LP model. If they did, we considered the two models to be connected. In addition, if two models are each connected to a third model, then we considered the original two models to be connected. By performing this test on many pairs of models, we were able to establish that all 1304 LP models are connected to one another, at least provisionally. (Note, however, that it does not imply that the set is convex, and, in fact, we found many pairs of models that had nonLP models on the line between them but that were connected via a path involving one or more intercalated models.) This provides strong evidence that our population of LP cells are all part of a single connected volume of LP cells in parameter space, rather than being in a number of islands of LP models separated by one or more “oceans” of nonLP models.
Polynomial fits of the property functions
It is highly desirable to have a quantitative description of how each parameter in a neuronal model influences each electrophysiological property. The computational model provides such a description only implicitly. A more direct description would be provided by a single equation that computes the magnitude of the electrophysiological property given the parameters. To provide such a description, we performed a multiple regression on each electrophysiological property, using the model parameters as predictors. This approach is similar to approaches used previously (Rabitz and Alis, 1999; Feng et al., 2004, 2006; Feng and Rabitz, 2004; Feng, 2005), but those approaches are not well suited to our LP population because our LP model population was not generated via independent variation of parameters (because of the filtering step). We first tried linear fits and quadratic fits to each electrophysiological property, but neither provided acceptable fits. A linear fit has the property that a given change in parameters will have the same effect on output, independent of the starting parameters; for instance, increasing ḡ_{A} by 3 μS/nF and ḡ_{pr} by 1 μS/nF will always increase the slowwave amplitude by 2 mV, regardless of the starting parameters. This will obviously not be true for many parametertoproperty maps arising in complicated biological systems, and allowing quadratic and higherorder terms in the polynomial fit allows a fit to avoid this constraint.
We then tried cubic fits to each electrophysiological property and found that these performed well: we obtained R^{2} values >0.85 for six of the eight properties we examined (Table 4). This indicates that these six properties, as functions of the parameters, were well approximated by cubic polynomials.
As an example, the coefficients of the fit to one of these, the slowwave amplitude, are shown in Figure 6. The quality of the fit is shown graphically in Figure 6A, which plots the true slowwave amplitude in each model versus the slowwave amplitude “predicted” by the fit. The points are near the identity line, indicating a good fit. Insight into which parameters influence the slowwave amplitude, and in what way, can be gained by examining the coefficients of the fit. The coefficients of the linear terms in the cubic fit show large contributions from ḡ_{A}, ḡ_{synPD}, ḡ_{pr}, P̄_{Ca}, ḡ_{Aa}, ḡ_{synAB}, and ḡ_{leak,ax} (Fig. 6B). The sign of these coefficients reveals whether the contribution of that term tends to increase the slowwave amplitude as that parameter is increased (positive sign), or decrease it (negative sign). Some of these effects are expected: because the PD synapse onto LP acts to hyperpolarize LP during the trough of its oscillation, one would expect that increasing ḡ_{synPD} would increase the slowwave amplitude. In addition, because the somatoneuritic I_{A} current would act to limit the rebound from inhibition, one might expect increasing somatoneuritic ḡ_{A} to decrease the slowwave amplitude. Other coefficients were more surprising: the linear coefficient of ḡ_{Aa} (the axonal form of I_{A}, which has generally faster kinetics than the somatoneuritic I_{A}) was large and positive, indicating that increasing ḡ_{Aa} tends to increase the slowwave amplitude. This was unexpected and is especially surprising in light of the fact that increasing the somatoneuritic I_{A} has the opposite effect. We found that this effect was caused primarily by a hyperpolarization of the trough potential in the axon (the lowest point of the slow wave), itself caused by the release from inactivation of I_{Aa} after the LP burst, which was then passively propagated to the soma. This constitutes a prediction of our model.
The quadratic coefficients of the cubic fit reveal nonlinear effects of changing model parameters (Fig. 6C). Two different kinds of quadratic terms in the cubic fit are possible: one kind involves the square of a single parameter (along the diagonal in Fig. 6C), and the other involves the product of two different parameters (off the diagonal in Fig. 6C). In the case of the slowwave amplitude fit, the largest quadratic coefficient was for the square of v_{1/2,pr}, the halfactivation of the proctolinactivated current. The coefficient for this term was negative, meaning that either increasing or decreasing v_{1/2,pr}, relative to its mean value will tend to decrease the slowwave amplitude, counteracting somewhat the linear effect of decreasing v_{1/2,pr}.
Cubic terms in the polynomial also describe nonlinear effects of changing model parameters, but involving products of three parameters, rather than two. These three parameters can all be the same, resulting in terms containing single parameters to the third power, or combinations of different parameters. In the slowwave amplitude fit, the largestmagnitude cubic coefficient is for the cube of ḡ_{synPD}, reflecting a nonlinear influence of ḡ_{synPD} on slowwave amplitude that somewhat counteracts the linear ḡ_{synPD} term (because the cubic coefficient is negative whereas the linear coefficient is positive and has a larger magnitude). These sorts of insights into what parameters affect other properties can also be obtained by examining the coefficients of those fits, at least in cases where the fit is good.
Using polynomial fits to quantify the effect of each parameter on each electrophysiological property
We wanted to quantify the overall contribution of each model parameter to determining each of the neuron's electrophysiological properties. We did this by starting with the coefficients of the final fit to a property, and setting all of the model parameters to their mean values. We then set them back to their actual values one parameter at a time and noted how much of the variance in the property was explained with the addition of each parameter, as a fraction of the total variance explained by the complete fit. Typically, the amount of variance accounted for by each parameter varies depending on the order in which parameters are added. Therefore, we repeated this process 3000 times, with a random order each time, and averaged the variance added over all repeats. This yielded a measure of the influence exerted by each parameter on the final fit such that the influence of all parameters summed to 1. These “influence factors” are shown in Figure 7 for each property we examined. Interestingly, we found that essentially all properties are affected by many parameters and that the relative importance of parameters varied from property to property.
Comparing polynomial fitting to traditional sensitivity analysis
Fitting a polynomial to a model property can be viewed as an alternative to more traditional sensitivity analysis. In sensitivity analysis, one typically starts with a single best model, and then each model parameter is varied, one at a time, and the change in the model property is plotted versus each parameter (Butera et al., 1999; Hill et al., 2001, 2002; Jezzini et al., 2004). In theory, fitting a polynomial gives the same information but also tells one how the model responds to combinations of parameter changes. Nevertheless, we were curious to see how the polynomial fits compared with sensitivity analysis. For individual models in the test set, we varied each parameter ±25% of its original value and performed simulations to see how model properties changed. Voltage parameters were varied ±5 mV. We compared the values in these simulations to the values predicted by the polynomial fit. An example is shown in Figure 8A, for the polynomial fit to the spike rate in the LP burst. This example displays good agreement between the cubic fit and the sensitivity analysis results.
In many cases, the polynomial fit lines matched the sensitivity analysis curves fairly well (Fig. 8B). To quantify the quality of these fits, we calculated R^{2} values for the cubic fit to each property in the neighborhood of each test model, using data from the sensitivity analysis around each test model. For this analysis, we were mainly interested in how well the cubic fits predicted changes from the property value at the test model, so we subtracted the fit property value at the central model from the fit values and similarly subtracted the simulation property value at the central model from the simulation values. (Failure to perform this step yielded R^{2} values that were generally much lower.) This analysis yielded a range of R^{2} values, but the median across test models was >0.5 for all but one property (phase of burst end, the property for which the cubic fit performed worst overall). This indicates that the cubic fits were able to predict the results of sensitivity analysis fairly well in most cases. However, in some cases, the cubic fits did not predict the sensitivity curves well, as indicated by low or even negative R^{2} values (a negative R^{2} value indicates a fit worse than would be achieved by simply fitting the data by a constant). Although the cubic fits are not a substitute for traditional sensitivity analysis, on the whole the ability of the cubic fits to predict the results of sensitivity analysis were striking, especially because the test models were not used in constructing the fits.
Discussion
How does the combination of voltagegated conductances in a neuron determine its electrophysiological properties? In this study, we analyzed a population of 1304 model neurons, all of which meet stringent criteria that match those of a specific biological neuron. We found (1) that a population of neurons can have a well constrained electrophysiological phenotype but still display disparate underlying maximal conductances with only weak correlations between conductances and (2) that each property of an identified neuron with well defined behavior is affected by a different set of maximal conductances. These results have implications for the interpretation of experimental results about the effect of neuronal parameters on electrophysiological properties.
Are correlations between maximal conductances required for a given electrophysiological phenotype?
The population of LP models did not show strong correlations between pairs of model parameters, including maximal conductance parameters. This is interesting, because experimental work found strong correlations between mRNA copy numbers for voltagegated channels in LP and other STG cells (Schulz et al., 2006, 2007). Although in some cases the mRNA levels correlate strongly with conductance magnitudes measured in voltage clamp (Schulz et al., 2006), it is not known whether mRNA copy numbers are well correlated with maximal conductances for all kinds of channels in all neuron classes.
There are two possible explanations for the lack of strong correlations in the model population: either the model fails to capture something important about the neuron's behavior or the neuron's electrical firing properties do not require the correlations seen in the molecular studies. Consistent with the latter explanation, it is possible that the correlations observed experimentally do not arise from electrophysiological constraints, since it is possible to construct LP neurons with conductance ratios other than those observed experimentally and still “act” like LP neurons electrophysiologically. If this were the case, the ratios observed experimentally might be a consequence of the transcriptional regulation of channel genes (Latchman, 1997; Fan et al., 2000; Brivanlou and Darnell, 2002; Rosati and McKinnon, 2004). For instance, having subunits of multiple channel types under the control of a common transcription factor might be simpler than independent transcriptional control of each channel type. This could result in different subunits being expressed in fixed ratios. Thus, the presence of correlations between conductances does not automatically imply that these correlations are required for proper electrophysiological behavior.
Using polynomial fits to capture the relative roles of multiple conductances on neuronal activity
The cubic fits allowed the quantification of the influence of each parameter on each electrophysiological property (Fig. 7). Each property was affected, at least to some extent, by many parameters in the model, but different properties were affected by various parameters to different extents. Recent work on a model of globus pallidus neurons yielded a similar result (Gunay et al., 2008), suggesting that this is a feature common to many neurons.
The large influence of the axonal conductances, and particularly the axonal leak and I_{A} conductance (ḡ_{Aa}), on the input conductance was somewhat unexpected. This suggests that, at least in the model, the axon is not as electrically irrelevant to the current–voltage relationship at the soma as might be supposed given the extreme attenuation of spikes in the soma. Although the contributions of multiple parameters on different intrinsic properties were generally different, those for the resting spike rate and the spike rate in the burst were similar. This reflects the fact that these two properties were strongly correlated in the LP model population (data not shown).
Small influence factors should be interpreted with caution; the synaptic conductances (ḡ_{synAB}, ḡ_{synPD}, and ḡ_{synPY}) had small, but nonzero, influence factors for the input conductance, resting membrane potential, and resting spike rate, despite the fact that all of these conductances were set to zero when these properties are measured in the model. This is possible because of highorder correlations between these parameters and parameters that do affect these properties. These correlations are exploited by the fitting routine to give significantly better fits than if the synaptic conductance parameters are not used (data not shown). These potentially misleading influence factors are uniformly small and do not detract from the utility of the influence factors overall.
The finding that each electrophysiological property is determined by a combination of several different conductances (Fig. 7) has important implications for how one interprets experiments involving pharmacological blockers, genetic deletions, and neuromodulation. When the levels of conductances vary from animal to animal, and when an electrophysiological property depends on multiple conductances, blocking one of these conductances may have highly variable effects on that property, because in some animals the conductance may already be low and be compensated for by the other conductances to keep the property at a normal level. A similar argument holds for genetic deletions, but with the added complication that if the deletions are nonacute, the animal may compensate for the knockout by regulating levels of other conductances (Guo et al., 2000; Zhang et al., 2002; Xu et al., 2003; Zhou et al., 2003). Similarly, neuromodulator effects may vary from preparation to preparation because neuromodulators act on a particular set of conductances, which may be high or low depending on the individual animal.
Approaches to understanding the effect of parameters on electrophysiological properties
In a model, the simplest approach to understanding the effect of parameters on electrophysiological properties is traditional sensitivity analysis: identifying one “best” set of parameters, varying one parameter at a time, and plotting the electrophysiological property of interest against the varied parameter. This approach has been fruitfully applied (Butera et al., 1999; Hill et al., 2001, 2002; Jezzini et al., 2004) but has two limitations: it is based on examining departures from a single model, and it is based on varying one parameter at a time, and therefore yields no information about simultaneously changing multiple parameters.
Another approach is to generate a population of models, either by picking random values for parameters or by systematically scanning the parameter space (Foster et al., 1993; Goldman et al., 2001; Prinz et al., 2003, 2004b; Achard and De Schutter, 2006; Tobin et al., 2006; Hobbs and Hooper, 2008; Weaver and Wearne, 2008). One can then use visualization techniques to determine how electrophysiological properties change as parameters are varied (Ward et al., 1994; Keim, 2000; Taylor et al., 2006), although these techniques are often qualitative and choosing a good method of visualization can be difficult.
In some cases, one can derive closedform expressions for model properties as functions of parameters (Abbott, 1990), and this is arguably the best solution, when it is possible.
In this work, we used cubic functions to fit electrophysiological properties as a function of model parameters. This allowed a relatively simple description of model behavior over a population of models (not just in the vicinity of a single model), and it described how these models respond to multiple simultaneous parameter changes. Similar methods have been used in the chemistry literature (Rabitz and Alis, 1999; Feng et al., 2004, 2006; Feng and Rabitz, 2004; Feng, 2005) but have not been commonly used to understand neuronal dynamics. We focused here on electrophysiological properties that were real functions of the parameters. For example, the slowwave amplitude is described by a single real number for each possible choice of parameters. Other cases of interest include properties that are integer (e.g., number of spikes per burst), category (e.g., activity type), or boolean (e.g., spiking or not) functions of the parameters. Fitting polynomials might be useful for integer properties if followed by a rounding operation but would not be ideal for category properties. For these sort of problems, an approach analogous to the one used here would be to use statistical classification (as opposed to regression) techniques (Hastie et al., 2001).
A related problem is finding the set of points in parameter space that all result in the same value for an electrophysiological property (the level set). Olypher and Calabrese (2007) used a numerical approximation to find the level sets of several electrophysiological properties. The method used here offers a more general answer to how parameters affect electrophysiological properties as long as the polynomial fits are sufficiently accurate. The cubic fits used here effectively describe all of the level sets of a particular electrophysiological property.
Relation to previous models of STG neurons
These LP models capture many important features of the biological LP neuron in C. borealis, including the response of LP to pyloric inputs, the attenuation of spikes in the soma, and LP's input impedance. Previous models of LP have captured some of these features but not all of them (Buchholtz et al., 1992; Nowotny et al., 2008). One difference between the behavior of the model and that of the biological LP neuron (Golowasch and Marder, 1992a) is that the model fires more rapidly than the biological neuron at rest, but without this we were not able to obtain the appropriate firing rates in response to inhibitory inputs.
Previous multicompartment models of STG neurons have had only two compartments, one to model the generation of spikes and another that contains the currents typically measured in voltageclamp experiments (SotoTrevino et al., 2005; Nowotny et al., 2008). Additional compartments were needed to accurately model both the input impedance of LP (data not shown) and the attenuation of spikes when measured in the soma.
Conclusions
In this work, we found that a population of neurons can have a well constrained electrophysiological phenotype but still display disparate underlying maximal conductances with only weak correlations between conductances, and that each property of an identified neuron with well defined behavior is affected by a different subset of several maximal conductances. These results have important implications for how one interprets experiments involving pharmacological blockers, genetic deletions, and neuromodulation.
Footnotes

This work was supported by National Institutes of Health Grants NS50928 (A.L.T.) and MH46742 and by the McDonnell Foundation (E.M.). The confocal imaging shown in Figure 2A was done in collaboration with Dirk Bucher.
 Correspondence should be addressed to Adam L. Taylor, Volen Center, Room 306, Brandeis University, MS 013, Waltham, MA 02454. altaylor{at}brandeis.edu