| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
Volen Center and Biology Department, Brandeis University, Waltham, Massachusetts 02454
Submitted 5 July 2003; accepted in final form 26 August 2003
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMETHODSRESULTSDISCUSSIONDISCLOSURESACKNOWLEDGMENTSREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONDISCLOSURESACKNOWLEDGMENTSREFERENCES |
|---|
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 |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONDISCLOSURESACKNOWLEDGMENTSREFERENCES |
|---|
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 
, where Ei is the reversal potential, A = 0.628 x 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
 |
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.
|
Together with the input current, the membrane currents govern the potential V according to
 |
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
 |
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
 |
) 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 |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONDISCLOSURESACKNOWLEDGMENTSREFERENCES |
|---|
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).
|
|
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.
|
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.
|
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.
|
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.
|
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.
|
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 x 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.
|
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.
|
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
