Journal of Neurophysiology

Alternative to Hand-Tuning Conductance-Based Models: Construction and Analysis of Databases of Model Neurons

Astrid A. Prinz, Cyrus P. Billimoria, Eve Marder

Abstract

Conventionally, the parameters of neuronal models are hand-tuned using trial-and-error searches to produce a desired behavior. Here, we present an alternative approach. We have generated a database of about 1.7 million single-compartment model neurons by independently varying 8 maximal membrane conductances based on measurements from lobster stomatogastric neurons. We classified the spontaneous electrical activity of each model neuron and its responsiveness to inputs during runtime with an adaptive algorithm and saved a reduced version of each neuron's activity pattern. Our analysis of the distribution of different activity types (silent, spiking, bursting, irregular) in the 8-dimensional conductance space indicates that the coarse grid of conductance values we chose is sufficient to capture the salient features of the distribution. The database can be searched for different combinations of neuron properties such as activity type, spike or burst frequency, resting potential, frequency–current relation, and phase-response curve. We demonstrate how the database can be screened for models that reproduce the behavior of a specific biological neuron and show that the contents of the database can give insight into the way a neuron's membrane conductances determine its activity pattern and response properties. Similar databases can be constructed to explore parameter spaces in multicompartmental models or small networks, or to examine the effects of changes in the voltage dependence of currents. In all cases, database searches can provide insight into how neuronal and network properties depend on the values of the parameters in the models.

INTRODUCTION

The spontaneous firing pattern of a neuron and how it responds to inputs from other neurons is crucially determined by the densities and dynamics of the ion channels in the neuron's membrane. These membrane conductances have a nonlinear dependence on the membrane potential, which itself is changed by the currents flowing through the conductances. A neuron with even a small number of membrane conductances is a complex dynamical system, and predicting the behavior of a cell with a physiologically realistic set of currents becomes very difficult.

Faced with ongoing channel turnover, neurons must constantly adjust their membrane currents to maintain their electrical identity (Marder and Prinz 2002). Experiments and simulations have shown that even small changes in one or a few currents can dramatically alter the activity of a neuron (De Schutter and Bower 1994; Goldman et al. 2001). On the other hand, similar activity can be achieved with widely varying sets of conductances in biological and model neurons (Bhalla and Bower 1993; De Schutter and Bower 1994; Foster et al. 1993; Goldman et al. 2001; Golowasch et al. 1999a), and conductance averages computed from neurons with similar behavior may not reproduce that behavior (Golowasch et al. 2002). Through a combination of mechanisms (Desai et al. 1999; Golowasch et al. 1999a,b; LeMasson et al. 1993; Liu et al. 1998; MacLean et al. 2003; Thoby-Brisson and Simmers 2002; Turrigiano et al. 1994), neurons maintain their conductances in a range that allows them to generate the desired activity. To understand how tightly a neuron's conductances need to be regulated and how this is achieved, it is essential that we know how the activity and response properties of nerve cells depend on their membrane conductances.

The question of how a cell's conductances influence its electrical behavior is usually addressed by varying one conductance in a given biological or model neuron (Golowasch and Marder 1992; Golowasch et al. 1992; Guckenheimer et al. 1993; McCormick and Huguenard 1992). In biological neurons this has been done by genetically modifying the density of one type of ion channel (Kupper et al. 2002; MacLean et al. 2003), by pharmacologically blocking one of the conductances (Debanne et al. 1997; Hoffman et al. 1997), or by adding a conductance with the dynamic clamp (Turrigiano et al. 1996). The behavior of the cell before and after the manipulation is then compared, and the observed change is interpreted as an effect of the altered conductance. However, the behavior of a nerve cell is the outcome of interactions between all of its membrane conductances (De Schutter and Bower 1994; Foster et al. 1993). It is therefore not obvious that changing the same conductance will have similar effects on the behavior of different neurons, irrespective of their other conductances. The result of varying one conductance at a time may thus depend on the particular choice of neuron that is studied (Foster et al. 1993).

We suggest that studying the relationship between a neuron's membrane conductances, its electrical activity, and its response properties is best done by generating a large number of computational model neurons that cover an extended range of values of all of the neuron's conductances. Here, we describe the construction of a database of about 1.7 million model neurons that were generated by independently varying the 8 maximal conductances of a realistic conductance-based model. The database contains information about each model neuron, including a reduced version of the neuron's voltage trace, its resting potential, spike frequency or burst frequency, frequency–current relation (fI curve), and, for bursters, the number of spikes per burst, burst duration, and phase-response curve. Because of this wealth of information, the database can be used to address a wide variety of questions, some of which are discussed here.

To our knowledge, the database method presented here is the first attempt to systematically explore a multidimensional neuronal parameter space and make extensive information about the steady-state voltage trace, the spontaneous behavior, and the response properties of every neuron available for future studies. Goldman et al. (2001) determined the type of activity for model neurons on a grid covering a 5-dimensional parameter space of a model related to ours, with the number of grid points corresponding to about 1% of the number of neurons contained in our database. Their study examined which combinations of membrane conductances determine whether a neuron is silent, fires tonically, or bursts. However, the conductance combinations explored in that study were constrained to a subspace of the whole model and the stored simulation results were too reduced to allow studies that go beyond the authors' original question. Similarly, Bhalla and Bower (1993) identified model neurons that mimic the behavior of olfactory bulb neurons by a “brute force” search of their model's parameter space. Again, however, their simulation results were not suitable for unrelated studies on the same model.

An obvious application of a model neuron database is the identification of neurons that reproduce a behavior observed in a biological cell. Traditionally, researchers have identified suitable model neurons for their biological neuron of interest by manually adjusting the parameters of a model (De Schutter and Bower 1994; Nadim et al. 1995; Traub et al. 1991). This process requires skill and experience and is often tedious and not always successful. If it fails, it is not clear whether the manual exploration was done in the wrong part of parameter space or whether the desired behavior cannot be achieved by the chosen model. We demonstrate how a database of model neurons can be searched for neurons that approximate a desired set of properties and can thus facilitate, complement, or replace manual searches of high-dimensional parameter spaces. Finally, we describe other potential database applications.

METHODS

Model

The model we used, based on experimental data obtained from lobster stomatogastric neurons (Turrigiano et al. 1995), was previously described (Prinz et al. 2003). Related models have been used before (Goldman et al. 2001; Liu et al. 1998). Each of the model's membrane currents is described by Math, where Ei is the reversal potential, A = 0.628 × 10–3 cm2 is the membrane area, and ḡi is the maximal conductance. To construct the database, the maximal conductances of all 8 currents in the model were varied independently to generate different model neurons (i.e., different combinations of maximal conductances). Reversal potentials were +50 mV for Na+, –80 mV for K+, –20 mV for IH, –50 mV for Ileak, and the Ca2+ reversal potential was determined from the momentary intracellular Ca2+ concentration using the Nernst equation and an extracellular Ca2+ concentration of 3 mM. The exponents pi are given in Table 1. The activation and inactivation variables mi and hi change according to Math with time constants τm and τh and steady-state values m and h as given in Table 1 (where voltages are in mV and times are in ms). Plots that show the voltage dependence of m, h, τm, and τh for all the voltage-dependent currents in the model are available on-line as supplemental material. IH was based on Huguenard and McCormick (1992) and all other currents are from Liu et al. (1998). All time constants were doubled to compensate for a temperature difference between the voltage clamp experiments by Turrigiano et al. (1995) and Huguenard and McCormick (1992) and the recordings reported here.

View this table:
table 1.

Voltage dependence of model currents

Together with the input current, the membrane currents govern the potential V according to Math where C = 0.628 nF is the membrane capacitance. For current step simulations, the input current was stepped from zero to a depolarizing DC current of 3 or 6 nA. For phase-response curves, the input consisted of pulses of inhibitory synaptic current Isyn = ḡsyn(V – Esyn), where ḡsyn is an instantaneously activating synaptic conductance. The synaptic current reversed at Esyn = –80 mV.

The intracellular Ca2+ concentration changes according to Math where τCa = 200 ms is the Ca2+ buffering time constant, f = 14.96 μM/nA translates the Ca2+ current into an intracellular concentration change (Liu et al. 1998), and [Ca2+]0 = 0.05 μM is the steady-state intracellular Ca2+ concentration if no Ca2+ flows across the membrane.

The differential equations for the membrane potential and the intracellular calcium concentration were integrated with the exponential method described in Dayan and Abbott (2001); those for the activation and inactivation variables were integrated with the Euler method. The numerical integration time step was 50 μs.

A previous study distinguished between tonically active neurons with narrow spikes and tonically active neurons with a broad shoulder after each spike (Goldman et al. 2001). The authors called the latter category “one-spike bursters,” arguing that the amount of transmitter such neurons release from a graded synapse is closer to that released by a bursting neuron than to that released by a spiker. To make a similar distinction, we continuously computed the integral T under the voltage trace between –40 and –15 mV Math The increase in T from one local voltage maximum to the next corresponds to the area under the peak between –40 and –15 mV. Goldman et al. (2001) interpreted the increase in T per discharge as a measure for the amount of transmitter the discharge would release from an instantaneous graded synapse with a transmission threshold of Vth = –40 mV and a saturation voltage of Vsat = –15 mV.

Model implementation

The model was implemented in C++ programs and executed on a shared Dell Beowulf cluster with 36 1.2-GHz processors. The computation time required for the simulation and classification of the spontaneous activity of all neurons was about 1 wk and resulted in 3 GB of data. The current injections took >1 mo and produced 23 GB, and the phase-response simulations were finished within a few days and generated 0.2 MB. With some code optimization and the growth in computer speed, we are confident that the computation time can be reduced in future database projects.

Phase-response curves

To obtain phase-response curves (PRCs) of the bursting model neurons, we simulated square pulses of synaptic conductance ḡsyn at different times during the rhythm with ḡsyn = 0 at all times before and after the pulse (Demir et al. 1997; Prinz et al. 2003). We assign phase 0 to the 1st maximum in each burst (burst onset) (Pinsker 1977). For a pulse beginning at time ΔT after the preceding burst onset, the stimulus phase is defined as ΔT/P, where P is the free-run burst period (Winfree 1980).

The period change ΔP caused by the pulse was the time difference between the 1st burst onset after the start of the pulse and the time at which this burst would have started if the pulse had not occurred. If the next following burst occurs earlier than in the free-run rhythm, ΔP < 0; if the burst starts later than in the unperturbed rhythm, ΔP > 0. When the normalized period change ΔP/P attributed to a pulse is plotted against the stimulus phase ΔT/P, the resulting curve is the classical PRC (Pinsker 1977).

Access to the database

The complete database described here is available for potential users. Because of the size of the database, downloading it from a website is impractical. We will therefore mail it on a set of DVDs on e-mail request to prinz{at}brandeis.edu. The DVDs contain conductance densities, spontaneous discharge patterns and discharge patterns under current injection, classification results, and spike or burst frequencies for all neurons, as well as burst durations, duty cycles, PRCs, and numbers of maxima per burst for the bursting neurons and resting potentials for the silent neurons. In addition, the DVDs provide documentation about the file formats used in the database and an executable that allows the simulation of a voltage trace for any neuron in the database.

RESULTS

Simulation and initial classification of spontaneous activity

We generated a database that contains information about the electrical activity of about 1.7 million model neurons with different maximal conductances of 8 Hodgkin–Huxley-type membrane currents. We call a particular set of conductances a “model neuron” or “neuron” and refer to the entire class of model neurons as the “model” (Goldman et al. 2001). The currents in this single-compartment model are based on those of lobster stomatogastric ganglion (STG) neurons (Turrigiano et al. 1995) and include a Na+ current, INa; 2 Ca2+ currents, ICaT and ICaS; a transient K+ current, IA; a Ca2+ -dependent K+ current, IKCa; a delayed rectifier K+ current, IKd; a hyperpolarization-activated inward current, IH; and a leak current, Ileak. The voltage dependence and dynamics of the currents and an intracellular calcium buffer included in the model are described in methods. To generate the database, we independently varied each conductance over 6 equidistant values, ranging from 0 mS/cm2 to a conductance-specific maximum, and collected results for all 68 = 1,679,616 possible combinations of these values. Maximum values were (in mS/cm2): 500 for INa, 12.5 for ICaT, 10 for ICaS, 50 for IA, 25 for IKCa, 125 for IKd, and 0.05 for IH and Ileak.

The spontaneous activity of each model neuron on this grid in conductance space was simulated and classified. Based on local voltage extrema (referred to as “extrema,” “maxima,” and “minima” from here on), the neurons were initially classified in 4 categories: silent, tonically active, bursting, and nonperiodic.

When simulating and classifying electrical activity for many neurons, it is crucial to minimize the amount of simulation per neuron. This is difficult because different neurons require different times to reach their steady state after initialization. In addition, the duration of spontaneous activity required for a correct classification varies between cells because of different discharge frequencies and activity types. We solved this problem by simulating spontaneous activity in epochs of 1 s and attempting to classify the activity after each epoch, stopping the simulation as soon as an unambiguous classification was achieved.

The adaptive algorithm for the simulation and 1st classification of the spontaneous activity for each model neuron is outlined in Fig. 1. For each neuron, the 13 dynamic model variables were initialized to V = –50 mV, [Ca2+] = 0.05 μM, m = 0 for activation, and h = 1 for inactivation variables. From this initial state, the approach of each neuron to its steady-state activity was simulated and local voltage maxima were counted. The approach to steady state was terminated when 500 maxima had been counted or after 10 s of simulated time, whichever occurred first (Fig. 1). The simulation time t and maxima count n were then reset to 0 and the spontaneous activity of the neuron was simulated and classified (Fig. 1). Spontaneous activity was simulated in epochs of 1 s while local voltage extrema were stored (gray box in Fig. 1). After each epoch, the stored extrema were used to attempt classification of the neuron as silent, tonically active, or bursting (Fig. 2).

fig. 1.

Flowchart for simulation of spontaneous activity. Simulation and classification algorithm was designed to minimize simulation time per neuron. n is number of voltage maxima, t is simulation time. After ≤10 s of simulated time allotted for reaching steady state, spontaneous activity was simulated in epochs of 1 s and classification was attempted after each epoch. Simulation was terminated as soon as a neuron was identified as tonically active or bursting. For silent neurons, simulation was continued until 20 s of spontaneous activity had passed without voltage extrema. To avoid misclassification of neurons that had not reached their steady state after 10 s or 500 maxima, t and n were reset after 20 epochs (one pass through gray box) and the algorithm in box was repeated. Neurons that had not been classified as silent, tonically active, or bursting after 4 passes through box were classified as nonperiodic.

fig. 2.

Classification of spontaneous activity. A–D: examples of steady-state activity categories used to classify all neurons in database. Horizontal gray bar indicates –50 mV, and time and voltage scales are identical for all traces. A: silent neurons have a steady-state voltage trace without extrema. Maximal conductances of this neuron are (in mS/cm2): Na 500, CaT 0, CaS 0, A 40, KCa 0, Kd 75, H 0.01, leak 0. B: tonically active neurons have extrema and interval between any 2 consecutive maxima differs by <1% from average interval. Maximal conductances (in mS/cm2): Na 100, CaT 0, CaS 10, A 40, KCa 50, Kd 75, H 0.02, leak 0.03. C: bursting neurons have periodic sequences of maxima and homologous intervals in consecutive periods differ by <1%. Maximal conductances (in mS/cm2): Na 100, CaT 0, CaS 4, A 0, KCa 15, Kd 50, H 0.02, leak 0.03. D: neurons that were neither silent, tonically active, nor bursting were classified as nonperiodic. Maximal conductances of this neuron (in mS/cm2): Na 300, CaT 0, CaS 10, A 20, KCa 20, Kd 125, H 0.05, leak 0.01. E: sample traces of tonically active neurons with (top) and without (bottom) a shoulder after peak. Insets: area (gray) under trace used to distinguish between spiking neurons and one-spike bursters. Numbers in insets give area and peak potential. Maximal conductances (in mS/cm2) are Na 100, CaT 0, CaS 4, A 10, KCa 10, Kd 75, H 0.01, leak 0.03 for spiker and Na 0, CaT 12.5, CaS 10, A 20, KCa 5, Kd 75, H 0.04, leak 0.03 for one-spike burster. F: area under peak, plotted against peak voltage for all tonically active neurons. There are 2 distinct groups of neurons. Those with a small area and a high peak potential were classified as spiking neurons; all others were classified as one-spike bursters.

A neuron was tentatively classified as silent if no extrema had been detected. If >10 maxima had been stored and all intermaximum intervals (IMIs) differed by <1% from their average, the neuron was classified as tonically active (Komendantov and Canavier 2002). The neuron was identified as a burster if n > 10 maxima had been stored and the voltage trace showed a periodic sequence of <n/2 IMIs with homologous intervals in consecutive periods differing by <1%.

Once a neuron was classified as tonically active or bursting, the simulation was terminated. For seemingly silent neurons, the simulation was continued by further epochs to ensure that the lack of detected extrema was not the result of a low firing rate. The simulation was also continued if ≤10 maxima had been stored or if the neuron could not be classified as tonically active or bursting in spite of >10 stored maxima.

After 20 s of simulated time or when 1,000 maxima had been stored, the simulation of spontaneous activity was halted. Neurons that still had no stored extrema were classified as silent.

Some neurons had not been classified as silent, tonically active, or bursting at this point. They belonged to 3 groups. 1) One group had such low discharge rates that ≤10 maxima occurred in the 20 s of simulated time and a reliable classification was not possible. 2) A second group was silent, tonically active, or bursting, but had not reached their steady state after the 10 s allotted for it. These neurons did not fulfill the classification criteria because of maxima stored before the steady state was reached. 3) A 3rd group did not fulfill the criteria for tonic activity or bursting because of truly nonperiodic firing patterns.

To separate group 2 from the others, the simulation time and maxima count were reset to 0 for all unclassified neurons. The gray box in Fig. 1 was then repeated up to 3 more times for these neurons, corresponding to a total of 90 s of simulated time or ≤4,500 maxima. After this, all neurons in group 2 had reached their steady state and could be classified.

Of the remaining unclassified neurons, those with >10 maxima were classified as nonperiodic. The spontaneous activity of those with ≤10 maxima was further simulated until 100 maxima had been stored and the activity could be classified.

Spontaneous activity data

Once a neuron was classified, a list with the time t, voltage V at the peak or trough, and synaptic release measure T for each stored extremum was saved to ASCII files along with a conductance code that uniquely identifies each neuron. Although these extrema contain salient information about the neurons' activity patterns, their storage does not require the prohibitive amounts of memory that would be needed if a complete voltage trace were saved for each neuron. In a separate file, the activity type—silent, tonically active, bursting, or nonperiodic—was saved for each neuron.

The approach to steady state was sometimes lengthy. At the end of the simulation of spontaneous activity, we therefore saved a snapshot of the 13 dynamic variables (the membrane potential, the intracellular calcium concentration, and the 7 activation variables and 4 inactivation variables) for each neuron. This allowed us to later reset all dynamic variables to their saved value and resume simulation at the same point, which greatly reduces the computation time necessary for additional simulations.

Reclassifications

The spontaneous activity classification results were checked for a random subset of neurons by inspecting their extrema. During this, we encountered several classes of misclassified neurons. One misclassified group consisted of silent neurons that took minutes to reach steady state. During this time their membrane potential showed damped oscillations with regularly spaced extrema (Fig. 3A), such that these neurons were mistaken as tonically active based on their IMIs. To fix these misclassifications, we searched all tonically active neurons for monotonically decreasing oscillation amplitude. Neurons for which the voltage difference between each saved maximum and its preceding minimum was monotonically decreasing were identified as potential slowly damped silent neurons. Their dynamic variables were reset to the stored snapshot values and the simulation was continued until the silent steady state was reached. The neurons were then reclassified as silent.

fig. 3.

Neurons that were difficult to classify. A: after initialization, this neuron showed damped oscillations that required >30 min to decay. Maximal conductances (in mS/cm2) are Na 0, CaT 0, CaS 4, A 40, KCa 10, Kd 100, H 0.02, leak 0.01. B: traces from a spiking neuron that initially showed no periodicity. After about 1 min, activity switched to regular spiking and remained regular as long as simulation was continued. Times in A and B refer to simulated time, not simulation time. Maximal conductances of this neuron are (in mS/cm2): Na 100, CaT 0, CaS 10, A 50, KCa 5, Kd 75, H 0.05, leak 0.03. C–D: traces from 2 neurons at far ends of a continuum of nonperiodic neurons. C: nonperiodic neuron reclassified as an irregular burster. Number of spikes in each burst is indicated. Discharge pattern looks regular at first sight, but close examination reveals no strict periodicity. Maximal conductances (in mS/cm2): Na 400, CaT 0, CaS 8, A 50, KCa 20, Kd 50, H 0.04, leak 0. D: irregular neuron. This model neuron showed no periodicity, even by less-strict criterion used to detect irregular bursters. Maximal conductances (in mS/cm2): Na 100, CaT 0, CaS 10, A 50, KCa 20, Kd 100, H 0.04, leak 0.02. Horizontal bars indicate –49 mV in A and –50 mV in BD.

Similar misclassifications occurred for a group of tonically active neurons that also reached their stable activity pattern only after minutes of simulated time and showed nonperiodic activity as they approached it (Fig. 3B). To identify such neurons, all nonperiodic neurons were reclassified based on the last 100 of their saved maxima. Neurons that showed tonic or bursting activity in this last part of the voltage series were reclassified accordingly. However, this reclassification is likely incomplete because it is possible that some neurons reached a periodic activity pattern so late that they appeared to be nonperiodic for all of the simulated 90 s or 4,500 maxima. In principle, this possibility remains for all irregular neurons, which make up about 0.5% of the neurons in the database.

A 3rd group of reclassifications was motivated by the fact that many nonperiodic neurons appeared like bursters with only slightly irregular firing patterns, whereas others showed no tendency to be periodic (Fig. 3, C and D). To distinguish them, we used a less-strict periodicity criterion: nonperiodic neurons were classified as irregular bursters if their IMI sequence showed bursts of maxima with the time between consecutive burst onsets differing by <10% from the average period (Fig. 3C). All other nonperiodic neurons were classified as truly irregular. This distinction between irregular bursting and truly irregular firing is somewhat arbitrary.

Many tonically active neurons had narrow spikes, whereas others showed pronounced shoulders after their discharge (Fig. 2E). To distinguish between the 2 groups, we computed the area under the peak between –40 and –15 mV for each tonically active neuron and plotted it against the peak voltage (Fig. 2, E and F and methods). The neurons fall into 2 distinct groups: one type of neuron has small areas under the discharge in spite of relatively depolarized peak potentials, whereas the other type has larger areas and shows a wide range of peak potentials. We defined tonically active neurons with an area <0.4 mVs and a peak >0 mV as spiking neurons and reclassified all other tonically active neurons as one-spike bursters (Goldman et al. 2001). Given that 2 simple features of the voltage trace, the area under the trace and the peak potential, are sufficient to distinguish between 2 types of neurons, we believe that this classification reflects 2 different discharge mechanisms, rather than artificially splitting a continuum of neurons into subcategories (Fig. 2F).

Additional spontaneous activity data processing

To facilitate searches of the database, we extracted several features from the stored extrema and saved them: for silent cells we saved the resting potential. For tonically active neurons we saved the discharge frequency (inverse of the IMI), and the transmitter release per period (increase in the integral T per period) at a graded synapse. For regular bursters we saved the burst period, the number of maxima per period, and the transmitter release per period. We also saved the number of spikes (maxima >0 mV) per burst period, the burst duration [burst period minus largest interspike interval (ISI)], and the duty cycle (burst duration divided by period). For nonperiodic neurons we saved the average discharge frequency (inverse of the average IMI).

We reduced the size of the database by removing redundant information from the saved extrema lists. The lists were reduced to the last 3 periods of extrema for tonically active or bursting neurons and the last 2,000 extrema for nonperiodic neurons.

Overview of model neurons in the database

Figure 4A illustrates the content of the database. Of the approximately 1.7 million neurons, 17% were silent, 16% were spiking, 67% were bursting (including 19% one-spike bursters and 3% irregular bursters), and 0.5% were irregular. Figure 4 shows that each of these categories contains a wide variety of neurons.

fig. 4.

Database contains a wide variety of neurons. A: pie chart showing percentage of database neurons in each category. B: histogram of resting potentials for silent neurons with logarithmic vertical scale. Silent category includes neurons with depolarized resting potentials. C: frequency histogram for spiking neurons. Insets: percentages of slow spikers (left) and fast spikers (right) with a given conductance density for ICaT, ICaS, and IKCa. Slow spikers have no ICaT and are likely to have small amounts of ICaS, whereas fast spikers tend to have a lot of ICaS and little or no IKCa. D: burst duration as determined from voltage maxima >0 mV, plotted against burst period for all regular bursters with >1 maxima per burst. E: histogram of number of voltage maxima per burst for all regular bursters. Note logarithmic vertical scale. Gray distribution was determined by counting all maxima per period; black distribution is based on maxima >0 mV. Insets: voltage traces for neurons with few (left) and many (right) maxima per period. Maximal conductances (in mS/cm2) are Na 400, CaT 2.5, CaS 4, A 50, KCa 25, Kd 75, H 0, leak 0.04 for left inset and Na 300, CaT 7.5, CaS 8, A 0, KCa 10, Kd 125, H 0.01, leak 0.03 for right inset.

Because the set of possible values each maximal membrane conductance could assume included 0 mS/cm2, the database contains pathological neurons such as ones with only inward currents. Many such pathological combinations were silent. This is why Fig. 4B, a histogram of resting potentials for all silent neurons, includes cells with resting potentials in the depolarized range. In the physiologically meaningful hyperpolarized range, the silent neurons cluster around the reversal potentials of the underlying conductances. Presumably, this is because one conductance is so dominant that it keeps the membrane potential close to its reversal potential and prevents the neuron from spiking or bursting.

Figure 4C shows a spike frequency histogram for the spiking neurons. The spike frequencies in the database range from 0.03 to 124 Hz, covering more than 4 orders of magnitude. The histogram shows 2 distinct distributions around 3.5 and 60 Hz with almost no overlap. We inspected the conductances of these 2 groups and found that the low-frequency spiking neurons had no ICaT and low densities of ICaS, whereas the high-frequency spiking neurons had higher densities of ICaS, but no or little IKCa. These differences suggest that the 2 neuron populations reflect 2 different mechanisms for supporting spiking. Separate regions of conductance space with low-frequency and high-frequency spiking neurons that differ in their calcium conductance were previously reported (LeMasson et al. 1993).

The burst durations of the regular bursters with >1 voltage maxima per burst are plotted against the burst period in Fig. 4D. The durations and periods cover wide value ranges. Just below the diagonal, there is a band of about 300 ms width with low burster density, suggesting that there is a lower limit for the interburst interval.

The burst durations in Fig. 4D were determined based on the timing of the voltage maxima >0 mV. This threshold was introduced to separate spikes from small-amplitude voltage maxima that we would be reluctant to call spikes (e.g., see the last maximum in each burst in Fig. 2C). Because the overall distribution of voltage values at local maxima in the entire database suggests no one threshold for spike detection (data not shown), this choice of threshold is arbitrary. However, the plot in Fig. 4D looked similar when we used all voltage maxima for the determination of the burst duration. The distribution of burst durations therefore seems to be fairly insensitive to the choice of spike threshold.

The gray histogram in Fig. 4E shows the number of maxima per burst for all regular bursters. Although 99% of the bursters have <30 maxima per burst, the number of maxima per burst ranges up to almost 500. When we again counted only the voltage maxima >0 mV, we obtained the black histogram in Fig. 4E. This distribution is similar to the histogram based on all maxima, but contains fewer neurons with >250 maxima per burst. This is attributed to elliptic bursters in the database (Fig. 5E). The voltage traces of these elliptic bursters show low-amplitude oscillations between bursts (Bertram et al. 1995; Izhikevich 2000; Wang and Rinzel 2003) that do not contribute to the number of maxima per burst if a threshold of 0 mV is applied.

fig. 5.

Diverse electrical behavior of database neurons. Examples that illustrate diversity of firing patterns found in database. Horizontal lines at left edge of each trace indicate –50 mV. A: neuron with unusually large discharge amplitude and width. Maximal conductances (in mS/cm2): Na 200, CaT 0, CaS 2, A 0, KCa 15, Kd 0, H 0.03, leak 0.04. B: neuron firing spike triplets. Maximal conductances (in mS/cm2): Na 100, CaT 0, CaS 10, A 50, KCa 10, Kd 50, H 0.03, leak 0.05. C: neuron with pronounced plateau at end of burst. Maximal conductances (in mS/cm2): Na 400, CaT 2.5, CaS 10, A 20, KCa 5, Kd 25, H 0.04, leak 0.03. Scale bar and –50 mV marker are same for A, B, and C. D: parabolic burster. Parabolic interspike interval (ISI) profile is superimposed on 2nd burst. Maximal conductances (in mS/cm2) are Na 100, CaT 0, CaS 6, A 10, KCa 10, Kd 50, H 0.03, leak 0.05. E: elliptic burster. Inset: voltage oscillations in interburst interval, expanded 4-fold in time and 100-fold in voltage with respect to main trace. Maximal conductances (in mS/cm2): Na 100, CaT 12.5, CaS 0, A 30, KCa 0, Kd 50, H 0.04, leak 0.02. F: low-amplitude oscillations. Arrows indicate beginning and end of one oscillation period. Maximal conductances (in mS/cm2): Na 0, CaT 0, CaS 6, A 20, KCa 25, Kd 0, H 0.02, leak 0.05. G: bursting neuron that alternates between 2 slightly different burst shapes. Number of spikes in each burst is indicated. Maximal conductances (in mS/cm2) are Na 500, CaT 2.5, CaS 8, A 0, KCa 15, Kd 75, H 0.05, leak 0. Scale is identical in F and G.

The diversity of firing patterns found in the database is further illustrated by Fig. 5. By varying the maximal conductances of otherwise unchanged membrane currents we obtained a large variety of different discharge shapes, including very high and wide action potentials (Fig. 5A), spike multiplets (Fig. 5B), and unusual burst shapes (Fig. 5C). The database also contains parabolic and elliptic bursters (Bertram et al. 1995; Izhikevich 2000; Wang and Rinzel 2003). A parabolic burster is shown in Fig. 5D; the ISIs within its burst first decrease and then increase again, leading to an approximately parabolic ISI profile. Figure 5E shows an example for an elliptic burster with small voltage oscillations in the interburst interval. The amplitude of the oscillations increases dramatically during the burst, leading to a voltage envelope reminiscent of an ellipse.

A subset of database neurons showed higher-order discharge periodicity. Examples are shown in Fig. 5, F and G. The neuron in Fig. 5F shows subthreshold oscillations whose amplitude repeats only every 6th maximum. Similarly, the cell in Fig. 5G alternates between 2 different interburst intervals, resulting in a period that contains 2 bursts of spikes.

Accuracy of the numerical integration

To allow the simulation and classification of almost 1.7 million model neurons within an acceptable time frame, we had to make a compromise between numerical accuracy and computation speed by selecting a time step of 50 μs for the numerical integration. To assess the impact of our choice of integration time step on the database as a whole, we randomly selected 10,000 neurons from the database and simulated their spontaneous activity with 5-μs time resolution and the 2nd-order method described in Dayan and Abbott (2001). The activity of these model neurons with a time step of 5 μs was then classified as described above, and the results were compared with the original classification from the simulation with 50 μs. The comparison showed that 99% of the 10,000 neurons had the same type of activity with both integration methods. Of the 107 (1%) neurons that had a different activity type under the 2 conditions, 55 were irregular in the database but were classified as spiking or bursting when simulated at high time resolution. Correspondingly, only 0.1% of the 10,000 test neurons were irregular at high time resolution, as opposed to 0.5% in the database itself. This suggests that many of the irregular neurons in the database were classified as irregular because of numerical artifacts and that this is responsible for more than half of the type differences at low compared with high time resolution.

We further quantified the influence of the integration method on our results by analyzing how the activity of each type of neuron was affected. For each of the 1,710 (17%) test neurons that were silent in the database and at high time resolution, we computed the difference between the resting potentials in the 2 conditions. The mean and SD of this difference were –0.08 and 47.67 μV, indicating that the silent neurons were practically identical in the 2 conditions. This is not surprising because the numerical integration method should affect mainly the dynamics, but not the silent steady state of the simulated neurons.

For the 1,639 (16%) test neurons that were spiking at high and low time resolution we compared the spike periods. For 86% of these neurons, the periods differed by <1% of the spike period at high resolution. The difference was smaller than 2% for 95%, smaller than 3% for 98%, and smaller than 4% for 99% of the spiking neurons. Using the low instead of the high time resolution tended to slow down the spiking neurons: 1,334 (81%) were slower in the database than at high resolution.

Similarly, we compared the burst periods of the 6,533 (65%) test neurons that were bursting in the database and at high time resolution. For 95% of these neurons, the burst period saved in the database differed by <3% from the high resolution burst period. The remaining 5% showed larger period deviations. In most of these cases, the burst period saved in the database was an integer multiple of the period found at high time resolution, or vice versa. This can occur if the burst pattern in one of the 2 cases shows higher-order periodicity, with a repeating sequence of 2 or more only slightly different burst periods. Such a sequence can be incorrectly detected as one long period by the fairly strict periodicity detection algorithm we used. As with the spiking neurons, most (i.e., 87%) of the bursters were slightly slower in the database than at high time resolution.

We also compared the number of maxima per burst at low and at high resolution and found that 93% of the regular bursters had the same number of maxima per burst under both conditions and 98% differed by 2 or less maxima. The numbers of maxima per burst with peak voltages above or below 0 mV were similarly sensitive to the time resolution of numerical integration: 95% of the regular bursters had the same number of spikes (maxima with peaks >0 mV) per burst at high and low time resolution and 98% differed by one spike or less. Compared with that, 94% of the regular bursters had the same number of maxima <0 mV per burst under both conditions and 99% differed by one or no maximum <0 mV.

It is conceivable that small accumulated errors resulting from the 50-μs time resolution could contribute to the long settling times and response times observed in some neurons (Figs. 2A and 8F). Although we did not systematically explore this possibility, we have evidence that the choice of time resolution may have a small influence on some of these slow processes. For example, the extremely slowly damped neuron in Fig. 2A reaches its silent steady state slightly faster at high than at low numerical time resolution. On the other hand, the long delay after inhibitory input shown in Fig. 8F is independent of the time resolution, suggesting that not all slow processes are subject to this type of error.

fig. 8.

Phase-response properties of regular bursters. A: definition of phase-response curve (PRC) parameters. Voltage trace of regular burster before, during, and after pulse of inhibitory synaptic conductance. The 1,000-nS stimulus arrived with a delay ΔT after previous burst onset and lasted for 25% of unperturbed burst period. For PRC, relative period change (P′ – P)/P = ΔP/P is plotted against stimulus phase ΔT/P. Positive values are delays, negative values advances. Maximal conductances of this burster are (in mS/cm2): Na 100, CaT 0, CaS 8, A 0, KCa 25, Kd 100, H 0.05, leak 0.01. Horizontal bar indicates –50 mV. B: example PRCs of regular bursters. Gray arrows indicate neutral phases (see text) for PRC 3 and PRC 4. Maximal conductances (in mS/cm2): Na 100, CaT 0, CaS 2, A 10, KCa 5, Kd 25, H 0, leak 0 for PRC 1; Na 400, CaT 0, CaS 6, A 30, KCa 0, Kd 100, H 0, leak 0.01 for PRC 2; Na 100, CaT 5, CaS 0, A 0, KCa 25, Kd 75, H 0, leak 0.02 for PRC 3; Na 400, CaT 0, CaS 6, A 30, KCa 20, Kd 25, H 0.01, leak 0.02 for PRC 4. C: histogram of percentage of period where inputs caused delays. D: distribution of neutral phases obtained by linear intrapolation from 2 neighboring PRC values. All bursters respond with delays to inputs in last 20% of their period. E: histogram of minimum (black) and maximum (gray) phase response ΔP/P. Because of scale, some very long delays are not visible. F: voltage trace from burster that responds to brief inhibitory input (indicated by bar below trace) with a long delay of next burst. PRC 1 in B is from this bursting neuron.

Taken together, these results show that the influence of the numerical integration time step on the overall database composition is fairly subtle. However, we believe it is important to quantify and be aware of these differences. Even with substantially higher computer speeds, the simulation of large numbers of numerical models will always have to deal with a trade-off between simulation speed and accuracy. A compromise between these 2 goals will therefore continue to be a part of database construction.

Sampling of conductance space

Our database covers an 8-dimensional conductance space with a grid of 6 values for each conductance. We addressed how well such a sparse sample captures the distribution of different types of activity in this space in 2 ways.

First, we examined families of neurons that differ in only one conductance. In conductance space, each such family corresponds to a line parallel to one of the conductance axes (Fig. 6). To gain insight into the distribution of the types of activity in our conductance space, we consider 2 extreme theoretical cases: in a simple conductance space, neurons of a given type (silent, spiking, bursting, or irregular) would occupy a single, continuous region, and all such regions would be of simple shape. In such a space, increasing one conductance while keeping all others constant means moving along a line that never returns to a region of one activity type after crossing through regions of other types. Figure 6A illustrates this for a 2D slice through the conductance space of our database. In contrast, a conductance space with a complex distribution of activity types could have small, discontinuous, or warped regions. In such a space, a line parallel to a conductance axis would likely cross in and out of a region of a given type multiple times. Figure 6B shows a slice through the database that contains such lines.

fig. 6.

Exploring conductance space with a sparse grid. A: 2D slice through a “well-behaved” region of conductance space. Conductances of INa and IH were varied; other conductances were held constant at (in mS/cm2): CaT 0, CaS 8, A 50, KCa 0, Kd 75, leak 0. Black lines show sampling grid. Symbols on grid points indicate spontaneous activity type of corresponding neuron. Bold black arrow shows a well-behaved family line crossing through different regions without returning to an activity type previously left behind. Two symbols off grid are random neurons. Top: located in a square with bursters on all corners and is itself a burster. Bottom: differs from its nearest neighbor on grid because it is near a boundary between regions of different types. B: 2D slice through a less well behaved region of conductance space. ICaS and IA were varied; other maximal conductances were (in mS/cm2): Na 400, CaT 5, KCa 10, Kd 125, H 0.02, leak 0.03. Bold arrow is a family line that encounters a burster, then an irregular neuron, and then returns to a region of bursters.

We previously defined a family of neurons as those that differ in only one conductance whereas the other 7 conductances are the same. For example, there are 67 families whose members differ only in their amount of INa. Analogously, there are 67 families for each of the 8 membrane conductances, resulting in a total of 8 × 67 = 2,239,488 families in the database. There are more families than neurons in the database because each family contains 6 neurons, whereas every neuron belongs to 8 families. Each neuron has one of the 4 activity types described above. Given these 4 types, 98% of the family lines in our database never returned to an activity region they had previously left behind. This suggests that the distribution of activity types in our conductance space is well-behaved and does not contain many disconnected or warped regions. Our sparse sample of this space is therefore likely to capture the salient features of the dependence of electrical activity on the underlying membrane conductances.

This result depends on our choice of activity types. The most arbitrary choice we made is the distinction between truly irregular neurons and irregular bursters, and our decision to include the latter in the burster category. We therefore repeated the family line analysis, but tentatively grouped the irregular bursters in the irregular category, resulting in 95% “well-behaved” families. Similarly, 98% of the lines were well-behaved when we grouped all irregular neurons in the burster category. Irrespective of this choice of categories, the space covered by the database seems to be well-behaved, suggesting that its structure can be studied on the basis of a sparse sampling grid.

Many of the family lines that were not “well-behaved” contained irregular neurons; examples for this can be seen in Fig. 6B. We have shown above that many of the irregular neurons in the database turn into spikers or bursters when the numerical integration is performed at a 10 times higher time resolution, suggesting that their irregularity is an artifact of the low time resolution we had to use for reasons of feasibility. For example, one of the 12 irregular neurons in the database slice in Fig. 6B turned out to be a regular burster at high time resolution. It is therefore conceivable that the true percentage of well-behaved family lines would be even higher than reported above if the entire database were simulated at high time resolution.

For our 2nd approach to the sampling question, we covered the range of conductance space defined by our grid with a new data set consisting of 10,000 neurons with random membrane conductances. We refer to this set of neurons as “random data set” and to its members as “random neurons.” If the sparse database grid captures the distribution of activity types in the space, one would expect most random neurons to be of the same type as neurons on nearby grid points.

We classified the activity of the random neurons as described for the neurons on the grid. Each random neuron is located in an 8-dimensional hypercube of conductance space whose 28 = 256 corners are defined by all possible combinations of the 2 grid values that bracket the random value of each conductance. In each dimension, a random point is either closer to the grid value below it or to the grid value above it. The combination of the 8 grid values closest to a random point therefore defines a unique nearest neighbor on the grid. Using the 4 activity types described above, we found that 90% of the random neurons had a nearest grid neighbor that showed the same activity type.

Furthermore, the 10% of random neurons that were not of the same type as their nearest grid neighbor were all located in hypercubes with corners of more than one type. This indicates that the activity type of these random neurons may differ from that of their nearest neighbors because they are located close to boundaries between conductance space regions of different types. A random neuron in Fig. 6A illustrates this in 2 dimensions.

The random neurons yielded an additional indication that the distribution of activity types is not very complex or discontinuous: we inspected random neurons that were located in a hypercube with the same activity type on all corners. All such random neurons showed the same type of activity as the corners of the hypercube in which they were located. This again indicates that there are few or no small-scale pockets of one activity type in regions of another type. Figure 6A shows an example in 2 dimensions.

Current injections

We further characterized all neurons by their response to depolarizing current steps of 3 and 6 nA. To simulate a current step, all dynamic variables were set to the values stored in the snapshot and the simulation was continued while storing extrema. For tonically active or bursting neurons, the current was stepped up when the simulation was halfway through the largest IMI of the period. For irregular neurons, the current step started at a local voltage minimum. The injection was continued and extrema were saved until a new steady state was reached, but for ≥1 s after the step (Fig. 7, A and B). The steady-state activity in response to 3- and 6-nA current injection was then characterized and saved as it was for the spontaneous activity. We tried to classify the neurons that were tonically active under current injection as either one-spike bursters or spiking neurons as in Fig. 2, E and F. However, there was no obvious separation of the tonically active neurons under current injection. The equivalents of Fig. 2F for 3- and 6-nA current injection showed the same 2 populations as Fig. 2F (neurons with narrow spikes and high peak potentials and neurons with broad shoulders after the discharge), but at peak potential between 30 and 50 mV the 2 populations were connected by neurons with a continuum of areas under the discharge (data not shown). We therefore treat one-spike bursters and spiking neurons both as tonically active in the following analysis.

fig. 7.

Current injections. A and B: activity of 2 neurons before and after injection current is stepped to 3 nA (top) and 6 nA (bottom). Maximal conductances (in mS/cm2) are Na 0, CaT 5, CaS 4, A 10, KCa 20, Kd 100, H 0.02, leak 0.03 in A and Na 400, CaT 2.5, CaS 10, A 50, KCa 20, Kd 0, H 0.04, leak 0 in B. C: number of neurons of activity types under 0 and 3 nA injection. D: number of neurons of 4 activity types under 3 injection conditions. E: number of voltage maxima during 1st s after step to 3 nA, plotted against steady-state discharge frequency during 3-nA injection. Neurons above diagonal show frequency adaptation. Insets: examples of neurons above and below diagonal. Maximal conductances (in mS/cm2) are Na 400, CaT 12.5, CaS 6, A 50, KCa 15, Kd 0, H 0.01, leak 0.02 for the top inset and Na 100, CaT 0, CaS 8, A 0, KCa 25, Kd 100, H 0.05, leak 0.01 for the bottom inset. F: steady-state fI curve examples. Curves can be linear (1) or have increasing (2, 4) or decreasing slope (3, 5). Maximal conductances of corresponding neurons are (in mS/cm2): Na 400, CaT 5, CaS 4, A 50, KCa 0, Kd 25, H 0.05, leak 0 for curve 1; Na 400, CaT 5, CaS 2, A 50, KCa 0, Kd 25, H 0.03, leak 0 for curve 2; Na 0, CaT 2.5, CaS 4, A 50, KCa 25, Kd 50, H 0.02, leak 0.04 for curve 4; Na 0, CaT 2.5, CaS 6, A 20, KCa 5, Kd 50, H 0.01, leak 0 for curve 5. Curve 3 belongs to neuron in A.

A total of 970,553 (58%) of the neurons showed the same type of steady-state activity whether no current, 3 nA, or 6 nA was injected. This was the case for 34% of the silent, 77% of the tonically active, 53% of the bursting, but for none of the irregular neurons. That no neuron showed irregular activity under all 3 injection conditions means that the firing of all irregular neurons in this database can be regularized by DC current injection. The remaining 709,063 (42%) neurons switched to a different activity type under one or both of the current injections, and 43,922 (3%) neurons showed 3 different types of activity under the 3 injection conditions. Figure 7C summarizes this for the 3-nA step. The majority of tonically active and bursting neurons retained the same type of activity under 3 nA, whereas most of the silent and irregular neurons switched to tonic activity (gray boxes in Fig. 7C). Similar results apply for the 6-nA step (data not shown). Figure 7D shows the numbers of silent, tonically active, bursting, and irregular neurons in the 3 injection conditions. With increasing current, more neurons became tonically active and fewer neurons were silent or bursting.

During the injection simulation, we saved the number of maxima (regardless of their peak voltage) during the first 1 s of injection and the steady-state discharge frequency to allow the measurement of adaptation. Figure 7E shows the number of maxima in the first 1 s after the step to 3 nA, plotted against the steady-state frequency during the 3-nA injection. Points above the diagonal are neurons that discharge more often in the first 1 s than during 1 s of steady-state activity, thus showing frequency adaptation. The neurons below the diagonal fire less during the first 1 s than later on. Based on this measure, 81% of all neurons showed adaptation when the injection current was stepped to 3 nA, and 77% adapted after a step to 6 nA. However, whether a neuron shows adaptation depends on the length of the spike-counting time window after the step. This is illustrated by the bottom inset in Fig. 7E, where a window of 0.1 s instead of 1 s would have resulted in an initial discharge frequency higher than the steady-state frequency.

The current injection data allow the construction of initial and steady-state 3-point fI curves for each neuron. Figure 7F shows examples for steady-state fI curves. The fI curves in the database show a variety of slopes and shapes and thus indicate that the neurons respond differently to the same input. fI curves with only 3 points can, however, give only a rough idea of the response properties of the neurons. For example, it is not possible to decide whether a neuron shows class 1 or class 2 excitability (Hodgkin 1948) based on these fI curves, that is, whether it can fire with arbitrarily low frequency.

An obvious extension of the current step simulations described here are current pulse simulations in which the injection current returns to 0 nA after a limited injection duration. Comparing the activity of database neurons before and after such a current pulse could be used to systematically search the database for bistable neurons. We have not simulated such current pulses in database neurons on a large scale, but we have evidence that some bistable neurons exist in the database from current pulse simulations in a selected subset of spikers (data not shown).

PRCs

The response of an oscillator to inputs depends on their timing relative to the oscillation (Brown and Eccles 1934). The response can be characterized by a phase-response curve (PRC) that describes the change in single-cycle period caused by inputs arriving at different times during the rhythm (Pinsker 1977; Pinsker and Kandel 1977). For the spontaneously, regularly bursting neurons in the database, we simulated PRCs in response to instantaneous, inhibitory synaptic inputs of 1,000-nS amplitude that lasted for 25% of the intrinsic burst period of the neuron (Fig. 8A). Such inputs were simulated at 10 equidistant phases of the ongoing rhythm, resulting in PRCs with a phase resolution of 0.1 (Fig. 8B).

Figure 8, B–E illustrate the distribution of phase-response properties of the bursters in the database. Of the 1,065,225 PRCs we computed, 65% showed monotonically increasing delays with increasing stimulus phase (like PRCs 1, 2, and 4 in Fig. 8B). Over 99% of the PRCs showed advances at some phases and delays at others, whereas the remaining 7,197 PRCs had delays at all tested phases (Fig. 8C). No PRCs showed advances for the entire phase range, so inhibitory synaptic inputs, if properly timed, can delay the next discharge of any regular burster in this model. At the phase where the PRC crosses the horizontal axis, the stimulus has no effect on the timing of the next burst—we call this the “neutral phase point.” The distribution of neutral phases shown in Fig. 8D drops off sharply between phase 0.6 and 0.7. This indicates that all bursters can be slowed by inhibitory synaptic inputs that occur during the last 20% of the burst period.

The range of advances and delays that can be achieved by inhibitory inputs to regular bursters is illustrated by Fig. 8E, which is a histogram of the minimum and maximum phase response ΔP/P in the PRC of each burster. Most of the minimum responses (largest advances) are between –0.7 and 0, and most of the maximum reponses (largest delays) are between 0.2 and 0.6. However, we found bursters that responded with extreme delays: for 1,671 bursters (0.16%), inhibitory input delayed the next burst by >1 period at all phases; in 609 bursters (0.06%), the delay was >10 periods; and in 74 bursters (0.01%), it was even >100 periods at all examined phases. Figure 8B shows an example of a PRC with delays of about 100 periods at all phases, and Fig. 8F shows a voltage trace for the same neuron. The burster remains hyperpolarized after the brief inhibitory input, slowly depolarizes over the next 250 s, and finally resumes bursting with a delay of >100 periods. Neurons with such extreme response properties are worth further investigation because these neurons are in a sense retaining a “memory” of previous inhibition for durations that far exceed the time constants of the neuron's conductances.

Example application

One application of a database of model neurons is the identification of model parameters that reproduce the behavior of a specific biological neuron. To illustrate this approach, we use the database introduced above to identify model neurons that reproduce the behavior of the pyloric pacemaker of the lobster. This pacemaker consists of an anterior burster (AB) neuron that is electrically coupled to 2 pyloric dilator (PD) neurons. Figure 9C shows a recording from a PD neuron that was obtained as described in Prinz et al. (2003) and illustrates the bursting activity exhibited by this pacemaker in its modulatory environment in the STG.

fig. 9.

Using database to identify pyloric pacemaker models. A: circles representing number of candidate pacemaker models at different stages of search. B: pyloric dilator (PD) neuron PRC (black dots) and model neuron PRCs (lines) in response to inhibitory synaptic conductance pulses lasting 25% of burst period in model neurons and 20 –30% of burst period in experiments (n = 3). Best-fitting model PRC is black, 34 other PRCs with acceptable fits are dark gray, and 45 PRCs that fit data less well are light gray. Maximal conductances for best fit (in mS/cm2) are Na 500, CaT 10, CaS 0, A 40, KCa 0, Kd 100, H 0.01, leak 0.04. C: PD neuron voltage trace with typical burst period, burst duration, duty cycle, and slow wave. Gray bars in C, D, and E indicate –60 mV. D: voltage trace of model neuron that gave best PRC fit in B. E: voltage trace of one of remaining 34 bursters with acceptable PRCs. Maximal conductances (in mS/cm2) are Na 200, CaT 5, CaS 4, A 40, KCa 5, Kd 125, H 0.01, leak 0. F: circles representing number of candidate neurons with INa densities given next to each circle after period was required to be between 1 and 2 s (black) and after burst duration was required to be between 0.5 and 0.75 s (white). Nonzero INa is necessary for burst durations in desired range. G: circles representing number of candidate neurons with IH densities given next to each circle after period and burst duration were required to be in right range (black) and after duty cycle and PRC criteria were applied (white). For duty cycles and PRCs similar to those of pyloric pacemaker, density of IH has to be around 0.01 mS/cm2.

To identify database neurons that reproduce the pacemaker's behavior, we started with the entire database and progressively refined our selection criteria to narrow the search to progressively more suitable subsets of neurons (Fig. 9). We first required that the selected neurons be bursters, reducing the number of candidates to 1,120,235. Next, we specified that the burst period be between 1 and 2 s, the range typical for pyloric pacemakers. This reduced the number of suitable neurons to 200,986. When we further required that the burster be regular and have a burst duration between 0.5 and 0.75 s, 8,047 neurons were left in the candidate pool. Limiting the range of duty cycles (burst duration divided by period) to 0.3–0.4 left 80 neurons that showed the desired spontaneous activity.

Next, we used the phase-response properties of the pyloric pacemaker to further limit the set of candidate neurons. Figure 9B shows PRC data recorded from 3 synaptically isolated pyloric pacemakers in response to inhibitory synaptic inputs with durations between 20 and 30% of the pyloric period [PRCs were recorded as described in Prinz et al. (2003)]. Superimposed are the PRCs of the 80 candidate neurons. Although none of the simulated PRCs fits the experimental data perfectly, there is a subset of PRCs that shows higher curvature than the rest and agrees better with the general shape of the experimental PRC. Using a cutoff criterion for the chi-square deviation of the model PRCs from the experimental PRC, we identified the 35 model bursters whose PRCs best fit the pyloric pacemaker PRC. They are shown in dark gray in Fig. 9B, among them the best fit in black.

The voltage traces of the neuron that gave the best PRC fit and of one other neuron from the group of 35 are shown in Fig. 9, D and E. Although both neurons fulfill all criteria we used thus far, the traces look different. Whereas the spikes in Fig. 9E ride on a slow wave, the trace is almost flat between bursts in Fig. 9D. As Fig. 9C shows, the pyloric pacemaker output consists of spikes occurring during the depolarized part of a slow wave voltage oscillation. In PD neurons, the lowest point of the slow wave typically lies between –70 and –50 mV, the highest point is between –55 and –25 mV, and the slow wave amplitude (difference between the highest and the lowest potential) ranges from 10 to 30 mV. We determined the lowest point, highest point (approximated by the last voltage maximum in the burst), and amplitude of the slow wave for the previously selected 35 model neurons and determined whether they fell in the ranges given above. Of the 35 candidates, 9 neurons had highest and lowest point, and amplitude in the right range, including the neuron in Fig. 9E, but not the one in Fig. 9D. These 9 neurons represent our final selection for model neurons that approximate the behavior of the pyloric pacemaker.

Figure 9A shows how the first 4 criteria reduce the number of suitable neurons step by step. In principle, the final result of this search—9 pyloric pacemakers models— could have been achieved in a single step by applying the criteria described above all at once. Why perform the search step by step? If the criteria are applied one at a time, inspection of the conductance distribution after each step can reveal information about the role of the underlying conductances in shaping the neuron's electrical activity. Figure 9, F and G show examples for this. The black circles in Fig. 9F illustrate that the bursters with periods between 1 and 2 s are almost evenly distributed over the 6 sodium conductance values, including the lowest value, 0 mS/cm2. Adding just one additional criterion— burst durations between 0.5 and 0.75 s— eliminates all neurons with 0 Na+ conductance (see white circles), implying that a nonzero Na+ current is not necessary for bursting with a period similar to the pyloric period, but is necessary for bursting with a burst duration in the range seen in PD neurons. This is reassuring because burst durations between 0.5 and 0.75 s can be achieved only with multiple spikes per burst, a behavior one would not expect in neurons without reasonable amounts of INa. Similarly, the circles in Fig. 9G show that bursting in the required range of periods and burst durations can be achieved with a wide distribution of IH densities, but that additional constraints on the duty cycle and the PRC narrow this distribution to just one of the IH conductance values covered in the database. Both the duty cycle criterion and the PRC criterion contribute to this narrowing: the duty cycle requirement alone leads to a distribution in which 54% of the candidates have an IH conductance of 0.01 mS/cm2, with the remaining 46% distributed nearly evenly over the other 5 values (data not shown). Additional application of the PRC criterion then reduces the number of candidates with 0.01 mS/cm2 from 43 to 35 and eliminates all other candidates (Fig. 9G). These findings suggest that IH plays an important role in controlling both the duty cycle and the phase-response properties of pacemaker neurons. This is consistent with previous results indicating that IH shapes the response of neural oscillators to hyperpolarizing perturbations (Prinz et al. 2003).

The maximal conductances of the 9 pyloric pacemaker models and their similarity to one another are shown in Fig. 10. Two of the 8 membrane conductances, those of IH and Ileak, are identical in all 9 pacemaker models, but the remaining conductances cover ranges of 2 to 4 of their 6 possible values. Seven of the pacemaker models are a direct grid neighbor of at least one other pacemaker model, where direct neighbors are neurons that have identical values for 7 of their conductances and adjacent values for the 8th. These 7 neurons form 2 pairs and one triplet of direct neighbors. The 2 pairs are in turn connected by pairs that have identical values for 6 of their 8 conductances, and the group of 4 they form together is linked to the triplet by several pairs that have 5 identical conductances. The remaining 2 pacemaker models (the first and the last in the list in Fig. 10) differ in 2 or 3 conductances from the nearest other pacemaker model. Taken together, the results shown in Fig. 10 indicate that the 9 pacemaker models are located in the same region of conductance space but do not form a directly connected, tight group. However, it is conceivable that all 9 pacemakers would be part of a continuous (although probably not convex) subspace if the “empty” space between them were sampled.

fig. 10.

Conductances and similarity of pyloric pacemaker models. Maximal conductances (left) and similarity matrix (right) of pyloric pacemaker models found in database. For every pair of neurons, similarity matrix shows how many conductances are identical in both neurons. Neurons are in same order in conductance table and in matrix. Maximal conductances of IH and Ileak (0.01 and 0 mS/cm2, respectively) were identical in all 9 models and are therefore not listed. Not all models are a direct neighbor of (i.e., have 7 identical conductances with) another pacemaker model, but every model has 5 or more conductances in common with at least one other pacemaker model.

DISCUSSION

Adaptive simulation and classification algorithm

Our adaptive simulation and classification algorithm simulates the activity of a model neuron in epochs of 1 s and attempts to classify the activity at the end of each epoch. For tonically active or bursting neurons, this limits the simulation time to the smallest number of epochs that permit classification.

The question of when to stop simulating is more difficult for nonperiodic neurons because the failure to identify a periodic sequence of IMIs in the simulated data can be the result either of truly nonperiodic activity or of periodic activity with a period longer than the duration of the simulation. Without an estimate of the maximum possible burst period in our model, this problem is unsolvable. We therefore limited the simulation for nonperiodic neurons to 90 s of simulated time. The longest burst period in the database is 43.3 s, which is much longer than the longest time constant in the model, and 99.9% of all identified bursters have periods below 7.6 s. We therefore assume that the neurons we classified as nonperiodic actually are nonperiodic, rather than bursters with a very long period.

Similarly, we limited the simulation time for silent neurons to 30 s. Given that 99.9% of the bursters in the database had a longest IMI below 3.7 s and that 99.9% of the tonically active neurons had periods <1 s, it is again unlikely that neurons classified as silent were misclassified because of very low discharge rates.

An additional and fundamental classification problem is caused by the potential for bistable neurons. We know that at least some bistable neurons exist in the database. We classified each neuron based on the stable activity pattern it approached after initialization of its dynamical variables. For bistable neurons, which can discharge in either of 2 different, stable patterns, this is necessarily incomplete, given that only one of the 2 stable states of the neuron is classified. This problem is not easily solved because there is no single stimulus that can switch all bistable neurons from one activity pattern to another. Searching the database for bistable neurons by simulating transient current pulses (as described in results) will therefore not necessarily identify all bistable neurons present in the database.

Sampling distributions in conductance space with a sparse grid

We sampled the distribution of silent, spiking, bursting, and irregular neurons in an 8-dimensional conductance space with a relatively sparse sampling grid of only 6 different conductance values in each dimension. The distribution of activity types in a space with more than 3 dimensions is hard to grasp, and further work is necessary to elucidate how the underlying conductances shape neuronal behavior. However, our analysis of the distribution of activity types in conductance space indicates that the salient features of such a distribution can be captured with a coarse sampling grid. This is fortunate from a computational viewpoint because it means that this high-dimensional parameter space can be characterized with a tolerable amount of computation.

The well-behaved nature of this conductance space is also fortunate from the viewpoint of the neuron. Neurons face ongoing turnover of their molecular components and must constantly regulate the channel densities in their membrane to maintain a stable electrical identity (Desai et al. 1999; Golowasch et al. 1999a,b; Marder and Prinz 2002; Thoby-Brisson and Simmers 2002). If the distribution of a given type of activity in conductance space were complex and discontinuous, small fluctuations in channel densities could lead to dramatic changes in the activity of a neuron, thus requiring cells to monitor and regulate their conductances very tightly. In contrast, experimental evidence suggests that biological neurons can achieve similar firing patterns with a continuum of different membrane conductances and can maintain their activity by regulating their current densities (Golowasch et al. 1999b; MacLean et al. 2003). This is understandable if the dependence of the electrical activity on the current densities is well-behaved in the above sense.

For feasibility reasons, databases of model neurons in high-dimensional parameter spaces necessarily must use sparse sampling, although this means that groups of neurons that fulfill a given set of requirements—like the 9 pyloric pacemaker models identified above—are likely surrounded by many more suitable neurons not covered by the sparse grid. There are several conceivable ways to identify these off-grid neurons and to thereby increase the number of neurons that meet a set of criteria. One possibility is to cover the hypercube surrounding each suitable grid neuron with a finer grid and apply the database method to it. Another possibility is to use the solution neurons identified on the sparse grid as starting points for manual parameter tuning to explore their vicinity in parameter space. In a 3rd option, one could scan lines connecting the solution neurons on the grid to neighboring grid points outside the solution region at higher resolution, thus searching for the high-dimensional boundaries of the solution region. Any off-grid neuron within those boundaries would then also be part of the solution region.

On using single-compartment models

The PRCs in Fig. 9B reveal that even the best 35 model PRCs found in our database do not perfectly agree with the experimental data. Instead, they show phase responses larger than most of the PD responses, with larger advances at early phases and larger delays at late phases. This may be an example of the limitations when using a database of single-compartment model neurons to mimic an electrotonically extended biological neuron or, as here, a group of electrically coupled neurons. In the pyloric pacemaker group, a hyperpolarizing input applied to a PD cell body will likely be attenuated before reaching the burst generation site in the AB neuron. It may thus cause a smaller phase change than the same input simulated in a single-compartment model neuron.

Despite this limitation, we believe that the identification of single-compartment models that approximate the behavior of an electrically extended biological neuron is useful. Even in cases where more elaborate models are needed to fully reproduce the behavior of a biological neuron in all detail, such single-compartment model neurons and their distribution in conductance space can be valuable. These models can be used as a starting point for the manual or automated construction of model neurons with several compartments. In this manner, we believe that the database method can facilitate or complement the traditional method of identifying suitable model neurons by manually tuning cellular parameters.

Furthermore, the general method of simulating and classifying the behavior of a large set of model neurons is not limited to conductance space. Although the first database we introduce here varies only the maximal conductances of a single-compartment model, the technique can easily be expanded to models with several compartments and/or the variation of other membrane parameters such as midpoints of activation or time constants.

Other potential applications

The database described here contains information about the spontaneous activity of a large number of model neurons and the way they respond to input. This information can be used to identify model neurons that mimic a desired set of behaviors, although this is only the simplest of many applications of such a database. The database can be used to search for any combination of conductance values, activity type, resting potential, spike or burst frequency, spikes per burst, fI curve, PRC, and any property extractable from the saved voltage extrema. Most such applications fall into 3 categories. The 1st category includes applications that identify model neurons with a given set of properties, as in the example above. The 2nd category of applications is concerned with analyzing the role of a model parameter in a given behavior. An example of this would be an analysis of the role of the 8 membrane conductances in gain modulation. This could be studied by searching the database for pairs of neurons that differ in only one of their conductances, but have different fI curve slopes (or “gains”). The existence of such pairs would indicate that changing just one conductance can affect the gain of the fI curve—such a conductance could be a target for gain modulation (Chance et al. 2002; Holt and Koch 1997; Smith et al. 2002). A 3rd category of database applications is concerned with the distribution of different behaviors in conductance space.

Alternative approaches

Using a database of computational models is not the only approach that can shed light on the dependence of a neuron's behavior on its membrane conductances. A number of studies have used nonlinear systems theory and bifurcation analysis to classify neural behaviors (mainly for different types of bursting) and determine their dependence on parameters of the underlying model (Guckenheimer et al. 1993, 1997; Izhikevich 2000; Rush and Rinzel 1994). Such studies are rewarding in models with a small number of dynamic variables that operate on clearly separable time scales. Models with a large number of variables with overlapping time scales, like the model used here, have to be simplified to make them amenable to bifurcation analysis (Canavier et al. 1991). In addition, the insights obtained from bifurcation analysis are often qualitative. For the quantitative analysis of physiologically realistic, multivariable conductance-based models, large-scale numerical simulations may be more fruitful than studies based on dynamical systems theory.

Foster et al. (1993) used a stochastic search method to locate regions of parameter space that fit certain predefined criteria, such as fI curve shape, to within a given tolerance. This approach generates valuable insight into the relative importance of different parameters for the desired behavior and into relationships between these parameters. However, the data generated with this approach are specific to the initial question and, in contrast to the data generated with our database method, cannot be used to answer other questions about the same system.

Conclusion and future applications

We have described the construction and preliminary analysis of a database of about 1.7 million model neurons and demonstrated that such a database can be used for the identification of suitable model neurons that mimic a desired set of behaviors. The single-compartment model we used is based on data from lobster STG neurons (Turrigiano et al. 1995), and we varied the maximal conductances of the membrane currents and characterized the behavior of each resulting model neuron without external stimuli, under DC current injection, and—for bursting neurons—in response to phasic inputs. The approach described here can be generalized in several ways. 1) Biophysical measurements from any neuronal type in any organism can be used to construct a database of model neurons appropriate for modeling that preparation. 2) The same strategy can be used to study the effects of conductances in different neuronal compartments (Mainen and Sejnowski 1996). 3) The database method can be used to study the influence of intrinsic and synaptic properties on the dynamics of small networks. 4) Similar analyses can be employed to study the effects of other model parameters such as activation and inactivation thresholds and slopes.

The method described here has only recently become practical as computer speed has increased and costs decreased. However, in years to come, as computers become increasingly faster, this kind of brute force approach should become part of the compendium of techniques available to neuroscientists wishing to construct models to inform our understanding of cellular level computation in the nervous system.

DISCLOSURES

This work was supported by Deutsche Forschungsgemeinschaft (PR 649/1-1) and National Institute of Mental Health Grant MH-46742.

Acknowledgments

We thank L. F. Abbott for helpful discussions and careful reading of the manuscript.

Footnotes

  • The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

References

View Abstract