## Abstract

Neurons, and realistic models of neurons, typically express several different types of voltage-gated conductances. These conductances are subject to continual regulation. Therefore it is essential to understand how changes in the conductances of a neuron affect its intrinsic properties, such as burst period or delay to firing after inhibition of a particular duration and magnitude. Even in model neurons, it can be difficult to visualize how the intrinsic properties vary as a function of their underlying maximal conductances. We used a technique, called clutter-based dimension reordering (CBDR), which enabled us to visualize intrinsic properties in high-dimensional conductance spaces. We applied CBDR to a family of models with eight different types of voltage- and calcium-dependent channels. CBDR yields images that reveal structure in the underlying conductance space. CBDR can also be used to visualize the results of other types of analysis. As examples, we use CBDR to visualize the results of a connected-components analysis, and to visually evaluate the results of a separating-hyperplane (i.e., linear classifier) analysis. We believe that CBDR will be a useful tool for visualizing the conductance spaces of neuronal models in many systems.

## INTRODUCTION

In recent years, it has been recognized that the voltage-gated conductances of a neuron are subject to continual regulation, which helps maintain a target activity level (Desai et al. 1999; Golowasch et al. 1999; LeMasson et al. 1993; Liu et al. 1998; MacLean et al. 2003; Turrigiano et al. 1994). This realization has spawned a number of questions. First, how is this done? That is, what are the subcellular mechanisms underlying this regulation? Second, *why* is it done in this way? What are the critical properties of the neuron that are being maintained in a particular regime? One key to addressing the second question is understanding how the voltage-gated conductances determine the intrinsic electrophysiological properties of the neuron, such as its firing rate in isolation, whether it is endogenously bursting, and the magnitude of postinhibitory rebound (PIR) the cell exhibits.

Intrinsic properties are determined by the maximal conductances and kinetics of the different voltage-gated conductances in the cell membrane. However, the mapping from maximal conductances and kinetics to intrinsic properties can be complex and nonlinear. It is sometimes assumed that two neurons with similar intrinsic properties must have similar maximal conductances, but recent work has emphasized that this is not the case (Goldman et al. 2001; Golowasch et al. 2002; Prinz et al. 2003). Furthermore, a particular intrinsic property is often attributed to a particular voltage-gated conductance, but it is clear that several voltage-gated conductances typically act in concert to produce a given intrinsic property (Goldman et al. 2001; Golowasch et al. 2002; Harris-Warrick et al. 1995; MacLean et al. 2003, 2005; Prinz et al. 2003). Therefore the mapping from maximal conductances to intrinsic properties is complex, and understanding how changes in the conductances influence intrinsic properties is difficult.

As an alternative to hand-tuning the parameters of a neuronal model, Prinz et al. (2003) systematically conducted a coarse-grained scan of the maximal conductances of a model with eight different voltage- and calcium-gated conductances. The model was based on detailed voltage-clamp measurements in cultured crustacean stomatogastric ganglion (STG) neurons (Turrigiano et al. 1995). Prinz et al. (2003) thus constructed a “database” of ∼1.7 million model neurons, which sampled from the eight-dimensional space of maximal conductances on a uniformly spaced grid. For each model in the database, a variety of information was stored, including its spontaneous activity type (e.g., silent, tonically spiking, bursting) and its response to various perturbations. Thus the database contains a good deal of information about the gross form of the mapping from maximal conductances to intrinsic properties.

Unfortunately, extracting this information in a form that gives insight into the relationship between maximal conductances and intrinsic properties is difficult. There are many analytical tools that can be brought to bear on this problem, such as dynamical systems theory (Guckenheimer and Holmes 1983; Strogatz 1994) and statistical machine learning (Bishop 1995; Hastie et al. 2003; Vapnik 1998), but we wished to find a very simple method to lay out all the models in the database “at a glance,” prior to performing these or other more in-depth analyses. The space of maximal conductances is eight dimensional, and so it is difficult to visualize how a given electrophysiological property changes as the maximal conductances change. In this paper, we apply a visualization technique, called clutter-based dimension reordering (CBDR) (Peng 2005; Peng et al. 2004), to the problem of visualizing the properties of a neuronal model in this eight-dimensional conductance space. This is the first time, to our knowledge, that CBDR has been applied in the context of neuronal modeling. CBDR presents all of the sampled models in a single image, without averaging, and reveals patterns in the database that were not previously apparent. Furthermore, it serves as a useful complement to other analysis techniques, such as connected-components analysis and linear classification (i.e., hyperplane fitting, a particular form of statistical machine learning). We find it useful for conducting a preliminary survey of a large database of model neurons.

## METHODS

The model used in the Prinz et al. (2003) database was derived from a model of cultured crustacean STG neurons (Turrigiano et al. 1995). Previous models had also used variants of this model (Goldman et al. 2001; Liu et al. 1998). The model contains eight different voltage- and calcium-gated conductances: a fast sodium conductance (*g*_{Na}); a fast, transient calcium conductance (*g*_{CaT}); a slow, transient calcium conductance (*g*_{CaS}); a fast, transient potassium conductance (*g*_{A}); a calcium-dependent potassium conductance (*g*_{KCa}); a delayed-rectifier potassium conductance (*g*_{Kd}); a hyperpolarization-activated mixed-cation conductance (*g*_{h}); and a leak conductance (*g*_{leak}). To construct the database, the maximal conductance of each was independently varied between zero and some conductance-specific maximum. The maximum value (in mS/cm^{2}) was 500 for *ḡ*_{Na}, 12.5 for *ḡ*_{CaT}, 10 for *ḡ*_{CaS}, 50 for *ḡ*_{A}, 25 for *ḡ*_{KCa}, 125 for *ḡ*_{Kd}, 0.05 for *ḡ*_{h}, and 0.05 for *ḡ*_{leak}. (We denote the maximal conductance of the sodium channel, for example, by *ḡ*_{Na}.) Each maximal conductance was varied in six steps from zero to the maximum, resulting in 6^{8} (≈1.7 million) different models. Each model was simulated, and several types of information about the behavior of each model were collected.

The spontaneous activity of each model was automatically classified into one of four categories: silent, tonically spiking, bursting, or irregular (Fig. 1) (see Prinz et al. 2003 for details on the classification algorithm). In addition, the bursting models were subcategorized according to how many local maxima of voltage occur during a single burst. (Each local maximum could be either an action potential or a subthreshold “hillock.”) Thus we sometimes speak of the “*n*-maxima-per-burst bursting models” or the “*n*-MPB bursting models,” meaning that subset of the bursting models that display *n* local maxima per burst in the voltage trace. Tonically spiking models were distinguished from one-maxima-per-burst models based on the area under the voltage trace between successive voltage minima (see Prinz et al. 2003). More details about the model and the database can be found in Prinz et al. (2003).

All visualization and analysis done here was performed using Matlab (The Mathworks, Natick MA).

### Connected-components analysis

We were interested in the number of “islands” of each activity type in conductance space. Assuming the sampling grid is dense enough (and we believe it is—see Prinz et al. (2003); see also discussion), one can use connectedness on the grid to drive conjectures about true connectedness in the underlying continuous conductance space. We investigated the number of islands of each behavior type on the grid using a standard connected-components algorithm from the image processing literature (Jain et al. 1995). This algorithm determines whether any two models of a given class on the sampled grid are *connected*. Two models are connected if there is a path (on the sampling grid) from one model to the other, where every model on the path is also part of the class in question. Intuitively, each connected component is a separate island (Fig. 2). Essentially, the connected components algorithm colors each model according to which connected component it is in. Any two models with the same color are connected; any two models with different colors are not.

In the definition of connected components given in the preceding text, two points are in the same connected component if and only if there is a path from one to the other such that all points in the path are in the same class. A “path” is defined as a list of grid points, such that subsequent points in the list are neighbors of each other. Two definitions of “neighbor” are in common usage. We used a restrictive definition, considering models that differed by one grid step in one dimension to be neighbors, but not models that differed by one grid step in multiple dimensions. In two dimensions, the restrictive definition implies that a single pixel has four neighbors (north, south, east, and west, called the four-neighborhood), whereas the more permissive criterion implies that a pixel has eight neighbors (the four-neighbors, plus northeast, southeast, southwest, and northwest, and called the eight-neighborhood). We chose the more restrictive definition because we found that the different behavior types were highly connected even using the more restrictive definition (see results).

We should stress that the connected components algorithm runs on the full eight-dimensional dataset and that the dimensional stack image is only a way to visualize the output (Fig. 2). In particular, we note that dimensional stack images often make things look disconnected that are not.

### Hyperplane fitting

We attempted to fit the boundaries of the various spontaneous activity types using hyperplanes, following Goldman et al. (2001). Given a vector of maximal conductances, *ḡ* = (*ḡ*_{Na},…, *ḡ*_{leak}), the equation of a hyperplane is (1) where *i* ranges over the different conductances. **w** is a vector orthogonal to the hyperplane boundary, and *b* is equal to the distance from the origin to the hyperplane along the direction of **w**, times the length of **w**. This is a single linear equation and so in eight dimensions defines a seven-dimensional hyperplane (Bishop 1995). For a given boundary (e.g., silent/nonsilent, tonic/nontonic, bursting/nonbursting, etc.), we attempted to find a **w** that provided an accurate prediction of whether a given *ḡ* was in the class or not in the class.

Following Goldman et al. (2001), we quantified the accuracy of the hyperplane fit using the *success rate* under the assumption that points were equally likely to be drawn from either class. This is simply the average of two fractions: the fraction of points in class A that are on the correct side of the hyperplane and fraction of points in class B that are on the other side of the hyperplane.

Although the success rate is a good measure of the hyperplane's goodness of fit, it has the disadvantage that it is not a continuous function of **w**, which creates difficulties for many optimization techniques. To find a **w** that provided good classification performance, we therefore chose to optimize the *negative cross-entropy* (NCE) of the classifier rather than the success rate (Bishop 1995). The NCE is the negative log-likelihood of the data given the parameters, so choosing parameters that minimize NCE selects a hyperplane that is the “most consistent” with the data. The NCE is a continuous function of **w**, and generally classifiers with high NCE tend to provide a high success rate. (Empirically, we find this to be the case for our classification problem.) The maximization was done using a Broyden-Fletcher-Goldfarb-Shanno quasi-newton routine (the fminunc() function in Matlab) (Bishop 1995).

### Multi-class logistic regression

In the description in the preceding text, we assume that we want to find a separating hyperplane between two different classes of models (e.g., tonic neurons and nontonic neurons). More generally, however, we really want to find a set of boundaries that divide the space into multiple regions (e.g., silent, tonic, bursting, and irregular). One of the simplest approaches to this problem is a linear discriminant. This is a generalization of the approach described in the preceding text. For each class, one defines a linear *discriminant function*, and the class assigned to any given point in the space is the class for which the discriminant function is largest. The linear discriminant function for class *k* has the form (2) The NCE can again be defined for this case and can be used to find a set of discriminant functions that effectively partition the conductance space. The boundaries between classes in this case are found at places where two classes are “tied” for having the largest discriminant function. i.e. the two classes have equal discriminant functions, and all other classes have lower discriminant functions. The condition that two classes, *i* and *j*, have equal discriminant functions can be written (3) To be a boundary point between two classes, the point must satisfy *Eq. 3* but also be such that all other discriminants are lower than *y _{i}*(

*ḡ*) and

*y*(

_{j}*ḡ*). This implies that the boundaries between classes will be hyperplanes but bounded by the regions where they intersect the hyperplanes that divide other pairs of classes. If there are only two classes, the boundary will be a simple hyperplane, defined by

**w**=

**w**

_{i}−

**w**

_{j}and

*b*=

*b*−

_{i}*b*. Thus in the case of two classes, the equation for the boundary between classes reduces to

_{j}*Eq. 1.*

## RESULTS

### Construction of the dimensional stack image

We used CBDR (Peng 2005; Peng et al. 2004) to visualize the Prinz et al. (2003) database. CBDR, in turn, uses a technique called *dimensional stacking*. Dimensional stacking is a technique for viewing a high dimensional data set in only two dimensions (LeBlanc et al. 1990). It is similar to the visualization of a three-dimensional data set as a montage of two-dimensional “slices” (Fig. 2). However, in dimensional stacking, the process of making montages is iterated, yielding a final image that is a montage of montages. For concreteness, we will describe the construction of a dimensional stack image of the database described in Prinz et al. (2003), an eight-dimensional dataset with six sampled points in each dimension. The starting point of a dimensional stack image is a grid scanning two maximal conductances with all other maximal conductances set to their minimum values (Fig. 3*A*). In Fig. 3*A*, the KCa and Na maximal conductances are scanned with the rest of the maximal conductances set to zero. Each square in this grid would then be colored according to some property of that point in the database, e.g., the spontaneous activity type. A 6 × 6 montage is then made of images of this sort (Fig. 3*B*) with each individual image in the montage scanning *ḡ*_{Kd} and *ḡ*_{Na} and with two additional conductances (here, *ḡ*_{Kd} and *ḡ*_{CaT}) varying over the montage. The original image (Fig. 3*A*) corresponds to the bottom-left image in the montage because the original image has both *ḡ*_{Kd} and *ḡ*_{CaT} set to zero. A 6 × 6 montage of 6 × 6 images yields a 36 × 36 image, which scans four conductances. A “second-order montage” is then made of these montages (Fig. 3*C*) and scans an additional two conductances (here *ḡ*_{h} and *ḡ*_{CaS}). This is a 6 × 6 grid of 36 × 36 montages, thus 216 × 216 pixels total. A third-order montage is then made from the second-order montages (Fig. 3*D*) and scans the final two conductances (here, *ḡ*_{leak} and *ḡ*_{A}). This is a 6 × 6 grid of 216 × 216 images, thus 1,296 × 1,296 pixels total, or ∼1.7 million pixels, equal to the number of models in the database. This third-order montage is the final dimensional stack image. In this way, a high-dimensional data set can be viewed, in its entirety, in two dimensions.

Another way to describe dimensional stacking is as a linear projection. Each point in the database can be described by its eight maximal conductances, scaled so that each value is an integer between 0 and 5 (as in Fig. 3). The pixel coordinates of that model, in the dimension stack image shown in Fig. 3*D*, would then be given by (4) (5) where *k̄*_{Na}, for instance, is the scaled maximal conductance of the fast sodium channel. Although this is a linear projection, it should be noted that the coefficients of the projection depend on the number of samples per dimension in the database. If, for instance, a new database was run with ten sample points per dimension, the coefficients of *Eqs. 4* and *5* would have to be powers of 10 to prevent overlapping of points in the dimensional stack.

### Optimization of the stack order

When a dimensional stack image is constructed, choices must be made about which conductances are to be the “high-order” conductances and which are to be “low-order” conductances. In Fig. 3, *ḡ*_{leak} and *ḡ _{A}* are the highest-order conductances, and

*ḡ*

_{KCa}and

*ḡ*

_{Na}are the lowest-order conductances. In addition, at each level one must decide which conductance will be on the

*x*axis and which on the

*y*axis. This set of choices is called the “stack order.” The stack order has a large impact on the usefulness of the resulting visualization (see following text). Because there are eight dimensions, there are 8! = 40,320 possible stack orders. (These can, however, be arranged into pairs which differ only by a switch of the

*x*and

*y*axes, so arguably the number of “meaningfully different” stack orders is only 8!/2 = 20,160.)

The central idea of CBDR is to define an easy-to-compute measure of the “usefulness” of a dimensional stack image and to then use some sort of search procedure to find a stack order that maximizes this measure. The measure of usefulness used in the original work on CBDR (Peng 2005; Peng et al. 2004) was only defined for binary (black-and-white) images, so we developed a different measure, which we call “edginess.” This is based on the idea that stack orders that yield large expanses of solid color are useful, and thus stack orders that consist of many small patches of color are less useful. For a candidate stack order, we generate the dimensional stack image and then count the number of pixel edges that have different colors on either side of the edge. This number is the edginess of the image. Edginess is actually a measure of the “uselessness” of the image, and so we seek to minimize it not maximize it.

We use a local search optimization routine to find an image with low edginess. The optimization proceeds by examining the edginess of all stack orders that are “neighbors” of the candidate stack order. A stack order is considered a neighbor of the candidate stack order if it can be generated by swapping two of the dimensions in the candidate order. Thus each stack order has 8 · 7/2 = 28 different neighbors. After the edginess of all neighbors has been computed, the neighbor with the lowest edginess is kept and becomes the new candidate stack order. If no neighbor has a lower edginess, the candidate stack order is considered the optimal stack order. (Strictly speaking, it is only *locally* optimal, but in practice this scheme seems to work well. Sometimes we find it necessary to start the routine from several (≤5) different random stack orders and pick the best of the final orders found.) Thus the algorithm moves from a random initial stack order to one that has low edginess.

A dimensional stack image of every model in the Prinz et al. (2003) database is shown in Fig. 4 *A1*. Each pixel represents a model, and the pixel's color specifies the spontaneous activity type of that model with bursting models subcategorized according to how many maxima per burst (MPB) they display (see methods). The mapping between color and activity type is shown in Fig. 4*C*, which also shows the fraction of models in the database that are of each type. Although some “structure” can be seen in Fig. 4*A1*, such as the bands of gray representing silent models, it leaves something to be desired. This gives rise to the need for stack-order optimization.

A sequence of stack orders arrived at by the optimization is shown in Fig. 4*A*. The initial stack order (Fig. 4*A1*) is essentially arbitrary. One can see that each successive image looks more “structured” than the last. The first image is mostly “snow,” with few patterns visible, but the final image (Fig. 4*A6*, shown again at a larger scale in *B*) has clear patterns in it.

The intuition behind the method is that conductances that have a big effect on the behaviors visualized will end up being high order, whereas conductances that have a small effect will end up being low order. If a particular conductance does not affect the behavior type, it will tend to be moved to the bottom of the stack order because in this way, long runs of solid color may be generated. In a sense, the high-order conductances are the ones to which the behaviors examined are sensitive and the low-order conductances are the ones to which the behaviors are less sensitive. Note that which conductances are high order and which are low order will completely depend on what the categories being visualized are. If, for instance, one was to visualize the time-average membrane potential of the models and then optimize the stack order, one would likely end up with a very different optimal stack order.

A number of statements about the “layout” of models in conductance space can be made, given this visualization (Fig. 4*B*). First, most of the nonbursting models are those with zero calcium-activated potassium conductance (*ḡ*_{KCa}). These correspond to the large vertical band of blue and gray models on the left-hand side of the image. Second, most one-maxima bursters have zero fast sodium conductance (*ḡ*_{Na}). These correspond to the large horizontal band of bright-green models at the bottom of the image. Many one-maxima busters that have nonzero *ḡ*_{Na} have zero delayed-rectifier potassium conductance (*ḡ*_{Kd}). These correspond to the thin bright-green horizontal bands in the image. Third, there seems to be a regular gradation of bursters from few maxima per burst to many maxima per burst as one increases *ḡ*_{Kd} and decreases *ḡ*_{CaT} (the fast, transient calcium conductance). This is represented by the repeated pattern of a gradient from green to yellow to orange to red found in the image. Fourth, the irregular bursters are almost all models with low but nonzero *ḡ*_{KCa} and correspond to models that “would be” low maxima-per-burst bursters if their *ḡ*_{KCa} was higher. Thus the dimensional stacking technique allows one to see structure in the data that is difficult to perceive otherwise.

### Visualizing the distribution of a subset of model neurons

Often one would like to understand how some set of constraints on the intrinsic properties of a neuron are reflected as constraints on the maximal conductances of that neuron. For instance, does a constraint on the intrinsic properties of a neuron strongly constrain the range of possible maximal conductances? Dimensional stacking provides an easy way to get insight into such questions.

We examined the distribution in conductance space of an interesting subset of models, the one-spike bursters with burst periods between 0.9 and 1.1 s (Fig. 5). This group is a subset of the 1-MPB bursters because it consists of 1-MPB bursters the voltage maximum of which are more depolarized than 0 mV and that have a burst period with the given range. (Using a threshold of 0 mV for counting a maximum as a “spike” is somewhat arbitrary, but reasonable.) Although these constraints yield models that behave similarly, there are many combinations of maximal conductances that can give rise to this behavior (Fig. 5*B*). Examining the dimensional stack image of this set (Fig. 5*A*), we see that models in this set span the full range of *ḡ*_{Na}, *ḡ*_{Kd}, *ḡ*_{A}, *ḡ*_{K}_{C}_{a}, *ḡ*_{CaT}, and *ḡ*_{CaS}. (It also turns out that they span the full range of *ḡ*_{h} and *ḡ*_{leak}, but this hard to see in the image.) Furthermore, it is possible to see some of the constraints on the maximal conductances that must hold for the model to be in the set. For instance, the “banding” pattern seen in the data is present because of the constraint that the period must be between 0.9 and 1.1 s, which imposes a constraint on the maximal conductances, particularly *ḡ*_{CaS}, ḡ_{CaT}, *ḡ*_{A}, and *ḡ*_{leak}.

In addition to being useful for seeing structure in a database, CBDR is also a useful complement to other techniques for gaining such insight. In the following text, we illustrate the use of CBDR in combination with two other techniques: connected-components analysis and hyperplane fitting (i.e., linear classification).

### Connected-components analysis

It is often of interest to know whether a given category of model neuron occupies a single connected region in the parameter space or whether neurons of the category occur in multiple, distinct “islands.” In the Prinz et al. (2003) database, we wondered whether each of the different behavior types (silent, tonic, bursting, irregular) found in this conductance space are connected. For example, is there one contiguous region of bursters in the space or do bursters occur in multiple, separated “islands”?

Assuming the uniformly sampled grid is dense enough, one can make reasonable conjectures about this question using connected-components analysis on the sampling grid. Furthermore, dimensional stacking and CBDR are useful tools for visualizing the results of such an analysis. For the silent, tonic, and bursting models, we found that all or nearly all the models were in a single connected component (on the sampling grid). We visualized the results of this analysis using dimensional stacking (Fig. 6). All the silent models were in a single connected component (Fig. 6*A*). For the tonic models, 87 connected components were found, but 99.96% of the tonic models were in the largest component (Fig. 6*B*, see *inset* for examples of small connected components). For the bursting models, 30 connected components were found, but 99.99% of the bursting neurons were in the largest component. We were concerned that the preponderance of bursting models in the database (63% of models in the database are bursting; see Fig. 4*C*) might cause them to appear connected even if they were not. To guard against this, we examined the connectivity of a subset of the bursting models, the 6–10 maxima-per-burst bursters. These comprise 12% of the models in the database. For this subset, there are 451 connected components, but 99.33% of models in this class were in the largest component. This suggests that the connectivity of the bursting neurons is not caused merely by their prevalence in the database. Thus all of the silent neurons, and nearly all of the tonic and bursting neurons, were each contained in a single connected component on the grid.

The irregular models, in contrast, had multiple components with a nonnegligible fraction of models in them (Fig. 7). For this group, there were 370 different connected components with 84.5% of the irregular neurons in the largest one. Although this is still a large fraction of the total, it falls short of the fractions found for the silent, tonic, and bursting neurons.

Based on these results, we conjecture that for each of the three major activity types (silent, tonically spiking, and bursting), there exists one connected region that includes ≥99% of the models of that type (measured by volume in conductance space). Although this is only a conjecture, and is subject to a variety of caveats (see discussion), it nevertheless seems to us a reasonable conjecture, based on the results of the connected-components analysis.

### Hyperplane fitting

Another way of understanding the structure of regions in a high-dimensional space is to use the methods of statistical learning theory (Bishop 1995; Hastie et al. 2003; Vapnik 1998). In the context of discrete categories (e.g., silent, tonic, bursting), all these methods essentially try to fit the boundaries between categories with mathematical functions of some sort (perhaps implicitly). Dimensional stacking can be a useful tool for visualizing such fits and gaining insight into why a given function does or does not fit the data.

Some of the simplest of statistical learning methods are collectively referred to as *linear classification* methods. These attempt to fit the boundaries between regions with linear functions, often called hyperplanes because they are the high-dimensional analog of a plane in three-dimensional space. We used hyperplanes to fit the boundary between different activity types and visualized the resulting fits using dimensional stacking (Fig. 8). We did this not on the full database but on the subset in which the models have no maximal conductances set to zero (Fig. 8*A*). This is a population of 5^{8} ≈ 400,000 models. We chose to focus on this subset because the borders of the activity regions are simpler for it than for the full database.

We found that the boundaries between activity types were well fit by hyperplanes over the ranges of maximal conductances found in the database. We quantified the quality of fit by calculating the success rate, a measure of how likely it is that a model is on the correct side of the hyperplane (see methods). We found that a hyperplane was able to separate the silent models from the nonsilent models (i.e., the tonic, burster, and irregular models, considered as a single group) with a success rate of 98.4% (Fig. 8*B*, *top*). This implies that the border between silent and nonsilent models is well approximated by a hyperplane (again, over the range of maximal conductances in the database). Fitting a hyperplane to the tonic/nontonic boundary produced a success rate of 89.3% (Fig. 8*B*, *middle*). This implies that a hyperplane is still a good approximation to the border between tonic and nontonic models, although not as good as for the silent/nonsilent models. Projecting the data onto the plane spanned by the normals to these two hyperplanes confirmed that the silent/nonsilent boundary was sharp and nearly linear, whereas the tonic/nontonic boundary was more gradual and displayed some amount of curvature (data not shown). This is consistent with the higher success rate for the silent/nonsilent hyperplane than for the tonic/nontonic hyperplane.

The fact that both the silent and tonic models can be reasonably well fit by hyperplanes suggests that these two hyperplanes might be combined to yield a good fit to the bursting models. We constructed a classifier using the silent and tonic hyperplanes and classifying a model as bursting if it was on the negative side of the tonic hyperplane *and* on the negative side of the silent hyperplane (Fig. 8*B*, *bottom*). This yielded a classifier with a success rate of 91.0%. Thus the bursting neurons are well fit by two hyperplane boundaries over the range of maximal conductances found in the database.

Combining the classifiers for silent/nonsilent, tonic/nontonic, and bursting/nonbursting together, we constructed a three-way classifier to predict the activity type. This classifier predicted a model as silent if it was on the positive side of the silent/nonsilent hyperplane. If a model was on the negative side of the silent/nonsilent hyperplane, it was predicted as tonic if it was on the positive side of the tonic/nontonic hyperplane and as bursting if on the negative side. This classifier achieved a success rate of 91.7% (treating irregulars as “don't cares”). Thus the silent/nonsilent and tonic/nontonic hyperplanes provide a complete partitioning of the conductance space over the ranges examined, and this partitioning is accurate for predicting the activity type of models.

This description of the layout of activity type is consistent with the view provided by the dimensional stack image (Fig. 8*A*). The bursting models seem to occupy a middle ground, with the silent models on one side and the tonic models on the other (arrows). This was confirmed by fitting a hyperplane to the silent/tonic boundary, treating bursters as “don't cares” (using only the no-zero-conductances subset). This hyperplane was able to divide the silent models and the tonic models perfectly, achieving a success rate of 100%. Inspection of the hyperplane showed that the boundary passed roughly through the middle of the bursting models (data not shown). Furthermore, we found that there are no silent models that neighbor a tonic model (again in the no-zero-conductances subset). Taken together, these facts confirm that the bursters lie between the silent and tonic models in this subset of the database.

As a further check, we also fit a three-class linear discriminant to the data, treating the irregulars as “don't cares” (see methods). This achieved a success rate of 91.2%, similar to the 91.7% rate achieved by the composite classifier in the preceding text. Examining the boundaries determined by this classifier, we found that its silent/burster boundary was very similar to the silent/nonsilent boundary determined above (RMS relative difference of the orthogonal vector was 0.6%). The tonic/burster boundary was also similar to the tonic/nontonic boundary determined above (RMS relative difference of the orthogonal vector was 13.0%). Thus the three-class classifier and the composite classifier provide very similar partitionings of the conductance space.

As mentioned in the preceding text, the high-order conductances are the ones to which the behavior type is, in a sense, most sensitive. Another way of gauging the sensitivity of behavioral transitions to different conductances is by examining the different components of **w** for the hyperplane fit to that transition. We wondered how well these two notions of “sensitivity” would accord with one another. For the tonic/nontonic boundary, the **w** vector is [0.19 −0.18 0.29 −0.09 −0.89 0.23 0.003 −0.02], where the order is [Na CaT CaS A KCa Kd h leak] and where the maximal conductances are all normalized to the range 0–1. Ordering the conductances by the absolute values of their respective **w** components, the order is [KCa CaS Kd Na CaT A leak h], which is roughly in accord with the ordering provided by the stack order optimization (Fig. 8*A*). We found, however, that the **w** vector for the silent/nonsilent boundary did not agree very well with the stack order shown in Fig. 8*A*. This is presumably because the silent models make up only a small fraction of the models shown in Fig. 8*A* and thus are not as important for the determination of the stack order.

## DISCUSSION

It is often assumed that a particular intrinsic property of a neuron is due to the presence of a particular voltage-gated conductance. For instance, postinhibitory rebound in a neuron is sometimes said to be caused by the hyperpolarization-activated (*g*_{h}) conductance. However, it is now generally recognized that different voltage-gated conductances can have overlapping effects on the intrinsic properties of a neuron (Gillessen and Alzheimer 1997; Goldman et al. 2001; Golowasch et al. 2002; Harris-Warrick et al. 1995; MacLean et al. 2003, 2005; Magee 1998; Ngo-Anh et al. 2005; Prinz et al. 2003; Ramakers and Storm 2002; Vervaeke et al. 2006). For instance, both *g*_{h} and the fast transient potassium conductance (*g*_{A}) can contribute to postinhibitory rebound (Harris-Warrick et al. 1995; MacLean et al. 2003, 2005). The fact that conductances have overlapping roles in producing intrinsic properties means that understanding how the maximal conductances of a neuron determine its intrinsic properties is much more difficult than it would be if each conductance determined a single intrinsic property. This question is made more interesting by the recognition that even identified neurons do not have fixed maximal conductances, either from animal to animal or within an individual animal's lifetime (Desai et al. 1999; Golowasch et al. 1999; LeMasson et al. 1993; Liu et al. 1998; MacLean et al. 2003; Schulz et al. 2006; Turrigiano et al. 1994). If neurons regulate their maximal conductances to maintain particular intrinsic properties (Marder and Prinz 2002; Turrigiano and Nelson 2004), then understanding how the maximal conductances determine the intrinsic properties will provide insight into which conductances must be co-regulated.

Understanding how the intrinsic properties of a neuron are related to its maximal conductances is difficult, even in model neurons, where in principle “everything is known” about the system. Ideally, mathematical analysis of a model would provide closed-form mathematical expressions for model properties as a function of the maximal conductances. This would then allow for insight into how these properties vary as a function of the maximal conductances. However, most realistic neuronal models are nonlinear, and deriving closed-form expressions for interesting model properties is difficult. This means that visualization of such properties can be essential for developing intuition about how the property varies as the conductances vary. This paper is a preliminary effort at applying such methods to a model database in a high-dimensional conductance space.

We applied the dimensional stacking technique to the problem of visualizing such a model database and used CBDR for optimizing the stack order to bring out salient features of the model space. We used these techniques to visualize the “layout” of spontaneous activity type in this space. Dimensional stacking complemented other techniques for understanding the geometry of spontaneous activity type: connected-components analysis and linear classification. Together, these methods allowed us to conjecture that nearly all silent, tonic, and bursting models are in a single connected component (each) and that the silent-bursting and bursting-tonic boundary are both approximately linear when attention is restricted to models with no zero conductances.

CBDR can reveal structure in the data that would not be apparent otherwise (Fig. 4*B*). By examining which conductances are chosen as the high-order axes, it is possible to determine which conductances have the greatest influence (in some sense) on the quantity plotted. The method makes no strong assumptions about what type of patterns are likely to be present in the data, unlike such methods as principal-components analysis and linear classifier analysis. Nonetheless, it can be used in conjunction with these methods and can be used to visualize the outputs of these methods (Fig. 8, see following text).

### Comparison with other visualization techniques

The visualization method presented here has many similarities to established techniques. Our method is a variant of the clutter-based dimension reordering technique (CBDR) (Peng 2005; Peng et al. 2004). This technique (as applied to dimensional stacking) involves optimizing a measure of the usefulness of a dimensional stack image by changing the stack order. The method seeks to minimize the *clutter* in an image, which is a measure of how not-useful a dimensional stack image is. The measure Peng et al. (2004) propose only works on binary (black-and-white) images, however, and differs from our edginess measure even when the latter is applied to binary images. The idea of optimizing the stack order to produce a useful visualization is similar to the idea of *projection pursuit* (Friedman 1987; Friedman and Tukey 1974). Projection pursuit (the visualization technique, not the regression method) involves projecting a swarm of points in *n* dimensions down to fewer dimensions and then using some measure of the usefulness of this projection to optimize the projection. The general idea of visualizing a high-dimensional space via embedding of lower-dimensional spaces is also used in the techniques of Mihalisin and colleagues (Mihalisin et al. 1990), and in the *Worlds within Worlds* approach (Feiner and Beshers 1990). All of the preceding methods are categorized as *geometric* methods by Keim (2000), who also discusses a variety of other methods for multidimensional visualization.

One way to view CBDR as applied to dimensional stack images is as dimensional stacking plus projection pursuit. Dimensional stacking is a linear projection of the *n*-dimensional space to a two-dimensional one. Different stack orders correspond to different projections. But the projections used in dimensional stacking have a specific feature: they take a uniformly spaced grid of points in *n* dimensions to a uniformly spaced grid of points in two dimensions. Preservation of the grid structure makes the edginess computation simple and fast, and the fact that the set of stack orders is discrete and finite makes the optimization problem simpler.

One advantage of dimensional stacking is that it allows all the models in the database to be visualized simultaneously. Alternative techniques, such as fixing all maximum conductances but three and then plotting the activity type as a function of the three remaining conductances, force one to make choices about what values of the other conductances to use, which could give a misleading picture of the database as a whole. Another technique would be to *average* over all the models corresponding to a given choice of the plotted conductance. That is, instead of picking a single value for the unplotted conductances, average over all possible values, perhaps using a weighted average of the colors chosen to represent each activity type. Again, however, the averaging technique does not allow the inspection of all models simultaneously. Furthermore, averaging necessarily obscures any patterns present in the averaged-over dimensions, which dimensional stacking preserves. This could be especially misleading when examining categorical variables with highly nonlinear boundaries (Golowasch et al. 2002).

A variant of the averaging technique was used by Goldman et al. (2001) in an analysis of a five-dimensional conductance space (see *Relation to previous work*). Those authors used multicolored symbols at each point to indicate the fraction of models at that point that were of each activity type. Although this technique is very useful, one must still specify in advance which variables will be placed on the axes, and which will not, and poor choices will result in less useful visualizations. (Those authors chose which axes to plot based on a separating-hyperplane analysis.) Using our technique, these sorts of choices are made automatically, and thus a useful visualization can be generated before one attempts to mathematically fit whatever structure may be present in the data.

### Relation to bifurcation theory

One approach to understanding how the dynamics of a neuronal model change as the maximal conductances (or other parameters) change is *bifurcation theory* (Guckenheimer and Holmes 1983; Hirsch et al. 2004; Strogatz 1994). Bifurcation theory allows one to find manifolds (curved surfaces) in the conductance space at which the system behavior changes qualitatively. For bifurcations of fixed points (essentially, silent models), in models of neurons, these manifolds can usually be described analytically in terms of the roots of a certain characteristic polynomial (Guckenheimer et al. 1993, 1997). Bifurcations of periodic orbits (e.g., from tonic to bursting behavior) must typically be computed numerically, and methods for doing so are an active area of research (Govaerts et al. 2005; Guckenheimer and Meloon 2000). However, an analysis of fixed-point bifurcations for a system under consideration can sometimes give indirect information about non-fixed-point bifurcations (Guckenheimer et al. 2005).

The approach taken here is not intended as a replacement for bifurcation theory. Rather it is intended to be used after one has generated a database that samples the parameter space on a coarsely sampled grid. One motivation for generating a database of this sort is that one wants to find models with particular properties but wants to put off decisions about what exactly those properties will be. With such a database in hand, CBDR is a fast method for seeing how model properties change as parameters are varied. It does not allow one to determine the locations of qualitative behavior changes with the precision offered by bifurcation theory. Another disadvantage of the database method is that the amount of computation required scales exponentially in the number of parameters. The database used here did not include a characterization of multistability that may be present in the models, something that bifurcation analysis handles quite readily. One advantage of the database method, however, is that it allows one to examine how any computable function of the model behavior varies with the parameters, not only the qualitative behavior type. (Although bifurcation analysis can provide quantitative information about, for instance, the frequency of tonic firing near a bifurcation.) Another advantage of the database method is its simplicity, since it requires comparatively little mathematical sophistication. Thus it is arguably better suited to the working biologist who wants to examine the range of behaviors possible in a complex model in a situation where computational resources are inexpensive and plentiful.

### Generalizations of the techniques

The techniques used here can be used on a wide variety of models. In this work, the free parameters that we examined were all maximal conductances, but this need not be so. Any set of parameters could be used to construct and optimize a dimensional stack image, including half-maximal activation voltages, time constants, etc. Similarly, our work examined a single-compartment model, but a multicompartmental model could also be subjected to dimensional stacking. There will in general be more free parameters in a multi-compartment model than in a single-compartment model, and so choices will have to be made about which parameters to do the dimensional stack on. In multicompartment models, however, the number of parameters usually scales more slowly than the number of compartments. This is because large groups of compartments typically have identical conductances. For instance, all axonal compartments or all the dendrites of a given class might have identical conductance densities, thus making the number of parameters considerably less than the number of compartments (Contreras et al. 1997; Traub et al. 1994, 2005).

In this work, we have focused on understanding the mapping from conductances to discrete variables. In Fig. 4, the quantity mapped is a combination category variable (activity type) and integer (number of maxima per burst). However, CBDR could easily be extended to deal with continuous-valued dependent variables. Such quantities could be represented as grayscale values, for instance, and the stack-order optimization could be based on minimizing squared differences between neighboring pixels. Such an analysis would provide a useful alternative to the parameter explorations often performed on conductance-based models (Hill et al. 2001, 2002; Olsen et al. 1995). Such analyses usually involve varying one parameter at a time around a central canonical model and examining how model properties change as a function of this. CBDR would enable the examination of simultaneous changes to multiple parameters and would allow for a determination of what parameters affected the model properties most strongly.

### Connected components

Knowing the number of connected components for a given behavior type is interesting for at least two reasons. *1*) If there are multiple components for a given behavior, they may correspond to distinct mechanisms for generating the behavior. *2*) If the system is subject to homeostatic regulation (Liu et al. 1998), it would presumably be easier to perform regulation within a connected component, whereas “jumping” from one connected component to another might be more difficult. Our discovery that nearly all silent, tonic, and bursting models in the database occupy a single connected component (each) is thus reassuring. In addition to finding that silent, tonic, and bursting models each occupy a single connected component, we found that, in a subset of the database, their borders are fairly simple, being well approximated by hyperplanes. This reinforces the results of the connected-components analysis and underscores the relative simplicity of the boundaries between activity types, at least over the ranges we examined.

It is important, however, to recognize what the connectedness results are *not*. We have determined the connectedness of models on a uniform grid in the conductance space. We have reason to believe that this grid samples the underlying continuous space reasonably well—Prinz et al. (2003) found that 90% of models sampled at random within the continuous space had the same activity type as the nearest grid point. But we have not *proved* anything in a mathematical sense about the connectedness of these regions in the underlying continuous space. It is possible, indeed likely, that some neighboring points on the grid that share the same activity type actually have regions of another activity type along the line segment connecting them. Conversely, some grid points that are not connected to each other on the grid could be connected in the underlying continuous space and only appear disconnected because we have not sampled finely enough. But given that the grid seems to sample the space reasonably well (see preceding text; see also Prinz et al. 2003), it seems likely that most of the models in the three main activity types are connected to one another in the continuous conductance space. It is still possible that there are regions of the space that contain many (perhaps infinitely many) small disconnected islands of a given activity type, but it nonetheless seems that *most* of the conductance space is not like this.

Another important caveat is that the database only contains information about a limited range of maximal conductances. Although the ranges were chosen to be roughly consistent with the biology, the results of the connected-components analysis and the hyperplane-fitting analysis must be interpreted with this caveat in mind. It is possible that there are separate “islands” of particular activity types that are outside the ranges scanned by the database. Similarly, we found that the boundaries between activity types are only linear over a particular range of conductances. We left out models with zero conductances when doing hyperplane fits precisely because hyperplanes provided poorer fits when the zero-conductance models were included. This indicates that the boundaries between activity types display nonnegligible curvature as the maximal conductances approach zero. However, even excluding the zero-conductance models, the maximal conductances scanned cover a fivefold range and the fact that the boundaries are well-approximated by hyperplanes over this range is noteworthy.

### Relation to previous work

The results reported here extend those of Goldman et al. (2001). In that work, the authors showed that there were approximate hyperplane boundaries between silent, tonic, and bursting model neurons, in a model family similar to the one we used. They did this in a five-dimensional conductance space in which *ḡ*_{Na}, *ḡ*_{CaT}, *ḡ*_{A}, *ḡ*_{KCa}, and *ḡ*_{Kd} were free to vary. (*ḡ*_{CaS} was set to 80% of *ḡ*_{CaT}, *ḡ*_{leak} was fixed, and *ḡ*_{h} was set to 0). Those authors also found that hyperplanes were able to effectively delineate different spontaneous activity types. This work extends that analysis to more dimensions and applies the CBDR for the visualization of this higher-dimensional model space.

### General conclusions

In this work, we have used CBDR for gaining insight into how the properties of model neurons vary as their parameters are varied. We found this tool to be useful in the context of a database of model neurons that sampled an eight-dimensional space of maximal conductances. Furthermore, we found that it provided a way of visualizing the “answers” given by such techniques as connected-components analysis and linear classification. Because the methods described here are easy to implement, they may be ideally suited for neurophysiologists wishing to gain intuitions about models constructed to capture the dynamics of specific neuronal types and thus may be useful adjuncts to other mathematical approaches to understanding the relationship between model parameters and the resulting neuronal dynamics.

## GRANTS

This work was supported by National Institute of Neurological Disorders and Stroke Grants NS-07292 and NS-50928 to A. L. Taylor, Burroughs-Wellcome Fund Career Award at the Scientific Interface to A. A. Prinz, and National Institute of Mental Health Grant MH-46742 and the McDonnell Foundation to E. Marder.

## Acknowledgments

The authors thank J. Langton, D. K. Wittenberg, E. Gifford, and J. Guckenheimer for helpful discussions.

## 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.

- Copyright © 2006 by the American Physiological Society