| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
INNOVATIVE METHODOLOGY
1Volen Center, 2Biology Department, and 3Computer Science Department, Brandeis University, Waltham, Massachusetts; and 4Biology Department, Emory University, Atlanta, Georgia
Submitted 7 April 2006; accepted in final form 26 April 2006
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
; 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 |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
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 (gNa); a fast, transient calcium conductance (gCaT); a slow, transient calcium conductance (gCaS); a fast, transient potassium conductance (gA); a calcium-dependent potassium conductance (gKCa); a delayed-rectifier potassium conductance (gKd); a hyperpolarization-activated mixed-cation conductance (gh); and a leak conductance (gleak). To construct the database, the maximal conductance of each was independently varied between zero and some conductance-specific maximum. The maximum value (in mS/cm2) 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 68 (
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).
|
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.
|
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) |
). 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) |
 | (3) |
) and yj(). 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 = wi – wj and b = bi – bj. Thus in the case of two classes, the equation for the boundary between classes reduces to Eq. 1. |
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
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. 3A). In Fig. 3A, 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 x 6 montage is then made of images of this sort (Fig. 3B) 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. 3A) corresponds to the bottom-left image in the montage because the original image has both Kd and CaT set to zero. A 6 x 6 montage of 6 x 6 images yields a 36 x 36 image, which scans four conductances. A "second-order montage" is then made of these montages (Fig. 3C) and scans an additional two conductances (here h and CaS). This is a 6 x 6 grid of 36 x 36 montages, thus 216 x 216 pixels total. A third-order montage is then made from the second-order montages (Fig. 3D) and scans the final two conductances (here, leak and A). This is a 6 x 6 grid of 216 x 216 images, thus 1,296 x 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.
|
 | (4) |
 | (5) |
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. 4C, which also shows the fraction of models in the database that are of each type. Although some "structure" can be seen in Fig. 4A1, 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.
|
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. 4B). 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. 5B). Examining the dimensional stack image of this set (Fig. 5A), we see that models in this set span the full range of Na, Kd, A, KCa, 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.
|
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. 6A). For the tonic models, 87 connected components were found, but 99.96% of the tonic models were in the largest component (Fig. 6B, 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. 4C) 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.
|
|
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. 8A). This is a population of 58 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.
|
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. 8B, 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. 8A). 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. 8A). We found, however, that the w vector for the silent/nonsilent boundary did not agree very well with the stack order shown in Fig. 8A. This is presumably because the silent models make up only a small fraction of the models shown in Fig. 8A and thus are not as important for the determination of the stack order.
|
DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
; 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 gh and the fast transient potassium conductance (gA) 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. 4B). 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 |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
ACKNOWLEDGMENTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
FOOTNOTES |
|---|
Address for reprint requests and other correspondence: A. L. Taylor, Volen Center, Rm 306, Brandeis University, MS 013, 415 South St., Waltham, MA 02454
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
Contreras D, Destexhe A, and Steriade M. Intracellular and computational characterization of the intracortical inhibitory control of synchronized thalamic inputs in vivo. J Neurophysiol 78: 335–350, 1997.
Desai NS, Rutherford LC, and Turrigiano GG. Plasticity in the intrinsic excitability of cortical pyramidal neurons. Nat Neurosci 2: 515–520, 1999.[CrossRef][Web of Science][Medline]
Feiner S and Beshers C. Worlds within worlds: metaphors for exploring n-dimensional virtual worlds. Proceedings of the Third Annual ACM Symposium on User Interface Software and Technology, edited by Olsen D, Snowbird UT, 1990, p. 76–83.
Friedman JH. Exploratory projection pursuit. J Am Stat Assoc 82: 249–266, 1987.[CrossRef][Web of Science]
Friedman JH and Tukey JW. A projection pursuit algorithm for exploratory data analysis. IEEE Trans Comput C 23: 881–890, 1974.
Gillessen T and Alzheimer C. Amplification of EPSPs by low Ni(2+)- and amiloride-sensitive Ca2+ channels in apical dendrites of rat CA1 pyramidal neurons. J Neurophysiol 77: 1639–1643, 1997.
Goldman MS, Golowasch J, Marder E, and Abbott LF. Global structure, robustness, and modulation of neuronal models. J Neurosci 21: 5229–5238, 2001.
Golowasch J, Abbott LF, and Marder E. Activity-dependent regulation of potassium currents in an identified neuron of the stomatogastric ganglion of the crab Cancer borealis. J Neurosci 19: 1–5, 1999.[Medline]
Golowasch J, Goldman MS, Abbott LF, and Marder E. Failure of averaging in the construction of a conductance-based neuron model. J Neurophysiol 87: 1129–1131, 2002.
Govaerts W, Kuznetsov YA, and Dhooge A. Numerical continuation of bifurcations of limit cycles in MATLAB. Siam J Sci Comput 27: 231–252, 2005.[CrossRef]
Guckenheimer J, Gueron S, and Harris-Warrick RM. Mapping the dynamics of a bursting neuron. Philos Trans R Soc Lond B Biol Sci 341: 345–359, 1993.[Web of Science][Medline]
Guckenheimer J and Holmes P. Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. New York: Springer-Verlag, 1983.
Guckenheimer J and Meloon B. Computing periodic orbits and their bifurcations with automatic differentiation. Siam J Sci Comput 22: 951–985, 2000.[CrossRef]
Guckenheimer J, Myers M, and Sturmfels B. Computing Hopf bifurcations I. Siam J Numer Anal 34: 1–21, 1997.
Guckenheimer J, Tien J, and Willms A. Bifurcations in the fast dynamics of neurons: Implications for bursting. In: Bursting: The Genesis of Rhythm in the Nervous System, edited by Coombes S and Bressloff PC: World Scientific, 2005, p. 89-122.
Harris-Warrick RM, Coniglio LM, Levini RM, Gueron S, and Guckenheimer J. Dopamine modulation of two subthreshold currents produces phase shifts in activity of an identified motoneuron. J Neurophysiol 74: 1404–1420, 1995.
Hastie T, Tibshirani R, and Friedman JH. The Elements of Statistical Learning. New York: Springer, 2003.
Hill AA, Lu J, Masino MA, Olsen OH, and Calabrese RL. A model of a segmental oscillator in the leech heartbeat neuronal network. J Comput Neurosci 10: 281–302, 2001.[CrossRef][Web of Science][Medline]
Hill AA, Masino MA, and Calabrese RL. Model of intersegmental coordination in the leech heartbeat neuronal network. J Neurophysiol 87: 1586–1602, 2002.
Hirsch MW, Smale S, and Devaney RL. Differential Equations, Dynamical Systems, and an Introduction to Chaos. San Diego, CA: Academic, 2004.
Jain R, Kasturi R, and Schunk B. Machine Vision. New York: McGraw Hill, 1995.
Keim D. An introduction to information visualization techniques for exploring large databases (tutorial given at Visualization ’00, Salt Lake City, UT), 2000.
LeBlanc J, Ward MO, and Wittels N. Exploring n-dimensional databases. Proceedings of the First IEEE Conference on Visualization ’90, edited by Kaufman A, 1990, p. 230–237.
LeMasson G, Marder E, and Abbott LF. Activity-dependent regulation of conductances in model neurons. Science 259: 1915–1917, 1993.
Liu Z, Golowasch J, Marder E, and Abbott LF. A model neuron with activity-dependent conductances regulated by multiple calcium sensors. J Neurosci 18: 2309–2320, 1998.
MacLean JN, Zhang Y, Goeritz ML, Casey R, Oliva R, Guckenheimer J, and Harris-Warrick RM. Activity-independent coregulation of IA and Ih in rhythmically active neurons. J Neurophysiol 94: 3601–3617, 2005.
MacLean JN, Zhang Y, Johnson BR, and Harris-Warrick RM. Activity-independent homeostasis in rhythmically active neurons. Neuron 37: 109–120, 2003.[CrossRef][Web of Science][Medline]
Magee JC. Dendritic hyperpolarization-activated currents modify the integrative properties of hippocampal CA1 pyramidal neurons. J Neurosci 18: 7613–7624, 1998.
Marder E and Prinz AA. Modeling stability in neuron and network function: the role of activity in homeostasis. Bioessays 24: 1145–1154, 2002.[CrossRef][Web of Science][Medline]
Mihalisin T, Gawlinksi E, Timlin J, and Schwegler J. Visualizing a scalar field on an n-dimensional lattice. Proceedings of the First IEEE Conference on Visualization ’90, edited by Kaufman A, 1990, p. 255–262.
Ngo-Anh TJ, Bloodgood BL, Lin M, Sabatini BL, Maylie J, and Adelman JP. SK channels and NMDA receptors form a Ca2+-mediated feedback loop in dendritic spines. Nat Neurosci 8: 642–649, 2005.[CrossRef][Web of Science][Medline]
Olsen OH, Nadim F, and Calabrese RL. Modeling the leech heartbeat elemental oscillator. II. Exploring the parameter space. J Comput Neurosci 2: 237–257, 1995.[CrossRef][Web of Science][Medline]
Peng W. Clutter-Based Dimension Reordering in Multi-Dimensional Data Visualization (Master's thesis). Worcester MA: Worcester Polytechnic Institute, 2005.
Peng W, Ward MO, and Rundensteiner EA. Clutter reduction in multi-dimensional data visualization using dimensional reordering. Proceedings of the IEEE Symposium on Information Visualization 2004, edited by Keahey A, Austin TX, 2004, p. 89–96.
Prinz AA, Billimoria CP, and Marder E. Alternative to hand-tuning conductance-based models: construction and analysis of databases of model neurons. J Neurophysiol 90: 3998–4015, 2003.
Ramakers GM and Storm JF. A postsynaptic transient K(+) current modulated by arachidonic acid regulates synaptic integration and threshold for LTP induction in hippocampal pyramidal cells. Proc Natl Acad Sci USA 99: 10144–10149, 2002.
Schulz DJ, Goaillard JM, and Marder E. Variable channel expression in identified single and electrically coupled neurons in different animals. Nat Neurosci 9: 356–362, 2006.[CrossRef][Web of Science][Medline]
Strogatz SH. Nonlinear Dynamics and Chaos. Reading, MA: Addison-Wesley, 1994.
Traub RD, Contreras D, Cunningham MO, Murray H, LeBeau FE, Roopun A, Bibbig A, Wilent WB, Higley MJ, and Whittington MA. Single-column thalamocortical network model exhibiting gamma oscillations, sleep spindles, and epileptogenic bursts. J Neurophysiol 93: 2194–2232, 2005.
Traub RD, Jefferys JG, Miles R, Whittington MA, and Toth K. A branching dendritic model of a rodent CA3 pyramidal neurone. J Physiol 481: 79–95, 1994.
Turrigiano G, Abbott LF, and Marder E. Activity-dependent changes in the intrinsic properties of cultured neurons. Science 264: 974–977, 1994.
Turrigiano GG, LeMasson G, and Marder E. Selective regulation of current densities underlies spontaneous changes in the activity of cultured neurons. J Neurosci 15: 3640–3652, 1995.[Abstract]
Turrigiano GG and Nelson SB. Homeostatic plasticity in the developing nervous system. Nat Rev Neurosci 5: 97–107, 2004.[CrossRef][Web of Science][Medline]
Vapnik V. Statistical Learning Theory. New York: John Wiley, 1998.
Vervaeke K, Hu H, Graham LJ, and Storm JF. Contrasting effects of the persistent Na+ current on neuronal excitability and spike timing. Neuron 49: 257–270, 2006.[CrossRef][Web of Science][Medline]
This article has been cited by other articles:
|
|
 |
|
  A. L. Taylor, J.-M. Goaillard, and E. Marder How Multiple Conductances Determine Electrophysiological Properties in a Multicompartment Model J. Neurosci., April 29, 2009; 29(17): 5573 - 5586. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 N. Keren, D. Bar-Yehuda, and A. Korngreen Experimentally guided modelling of dendritic excitability in rat neocortical pyramidal neurones J. Physiol., April 1, 2009; 587(7): 1413 - 1437. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 E. Hong, F. Gurel Kazanci, and A. A. Prinz Different Roles of Related Currents in Fast and Slow Spiking of Model Neurons From Two Phyla J Neurophysiol, October 1, 2008; 100(4): 2048 - 2061. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
A. V. Olypher and R. L. Calabrese Using Constraints on Neuronal Activity to Reveal Compensatory Changes in Neuronal Parameters J Neurophysiol, December 1, 2007; 98(6): 3749 - 3758. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
R. J. Calin-Jageman, M. J. Tunstall, B. D. Mensh, P. S. Katz, and W. N. Frost Parameter Space Analysis Suggests Multi-Site Plasticity Contributes to Motor Pattern Initiation in Tritonia J Neurophysiol, October 1, 2007; 98(4): 2382 - 2398. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
V. Matveev, R. Bertram, and A. Sherman Residual Bound Ca2+ Can Account for the Effects of Ca2+ Buffers on Synaptic Facilitation J Neurophysiol, December 1, 2006; 96(6): 3389 - 3397. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| Visit Other APS Journals Online |
