| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
Faculty of Life Sciences and the Leslie and Susan Gonda Interdisciplinary Brain Research Center, Bar-Ilan University, Ramat-Gan, Israel
Submitted 20 April 2005; accepted in final form 7 August 2005
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
; Johnston et al. 2003; Migliore and Shepherd 2002; Stuart et al. 1999), mainly as a result of patch-clamp recordings from visually identified dendrites in brain slices (Stuart et al. 1993) and novel imaging techniques (Antic 2003; Antic et al. 1999; Denk et al. 1994; Lasser-Ross et al. 1991; Tsien 1989). Action potentials initiated at or near the soma actively back-propagate into the dendritic tree (Bischofberger and Jonas 1997; Chen et al. 1997, 2002; Häusser et al. 1995; Spruston et al. 1995; Stuart and Sakmann 1994). Furthermore, dendrites generate complex regenerative Ca2+ and Na+ spikes (Amitai et al. 1993; Antic 2003; Ariav et al. 2003; Bischofberger and Jonas 1997; Johnston et al. 1996, 2003; Magee et al. 1996, 1998; Martina et al. 2000; Migliore and Shepherd 2002; Schiller et al. 1997; Zhu 2000), modulate synaptic potentials (Magee 1999; Magee and Johnston 1995), contain both electrically and chemically defined compartments (Bekkers 2000b; Hoffman et al. 1997; Korngreen and Sakmann 2000; Larkum et al. 1999, 2001; Magee 1999; Schiller et al. 1997, 2000), and influence the induction and expression of synaptic plasticity (Golding et al. 2002; Magee and Johnston 1997; Markram et al. 1997). This exciting new information revived the discussion about the computational properties of single neurons (Häusser and Mel 2003; Mel 2002; Poirazi and Mel 2001; Polsky et al. 2004).
Clearly, understanding the computations carried out by complex neurons requires functional models. These models must take into account several types of voltage-gated conductances, some of them nonhomogeneously distributed on the neuronal membrane. For example, at least ten voltage-gated conductances were isolated from layer 5 (L5) pyramidal neurons (Bekkers 2000a; Berger et al. 2001; Brown et al. 1993, 1994a,b; Colbert and Pan 2002; Fleidervish and Gutnick 1996; Fleidervish et al. 1996; Foehring 1996; Foehring and Armstrong 1996; Foehring and Scroggs 1994; Foehring and Waters 1991; Foehring et al. 2000; Korngreen and Sakmann 2000; Sayer et al. 1990, 1993). Some of these conductances were observed to have nonhomogeneous density distributions along the apical dendrite of the neurons (Bekkers 2000b; Berger et al. 2001; Korngreen and Sakmann 2000; Williams and Stuart 2000). Finally, as a result of profound experimental constraints, the density and properties of the voltage-gated conductances are known only for the soma and apical dendrite. Conductance densities along basal dendrites and the axon are not known. Therefore many parameters in models of L5 pyramidal neurons were estimated by hand (Mainen et al. 1995; Rhodes and Llinás 2001; Schaefer et al. 2003). Because of the large number of free parameters current models of L5 neurons are ill constrained (Mainen et al. 1995; Rhodes and Llinás 2001; Schaefer et al. 2003). Similar problems have been observed in other cell types.
Previous studies show that parameter search methods can be effective in finding matches between single-neuron models and a target data set (Eichler West 1996; Vanier and Bower 1999). Not all parameter search methods perform equally well, and the relative performance of the different methods depends on the model being optimized. Stochastic algorithms such as simulated annealing (Kirkpatrick et al. 1983) and genetic algorithms (GAs) (Mitchell 1996) have been found to be the most successful in constraining single-neuron compartmental models because such algorithms are less greedy and rarely fall into local minima, unlike gradient-descent methods that easily move into a local minimum in parameter space and stay there (Vanier and Bower 1999). For a search with a large number of parameters, the GA method, which can simultaneously search many points in parameter space, may be the most effective. The intrinsic parallelism of this method also makes it very well suited for implementation on parallel computers (Eichler West 1996).
Finding the best cost function to mach the data in the parameter space is one of the most important conditions for effective convergence of any optimization algorithm. It influences the speed of convergence, the performance and execution time of the process, and determines the overall behavior of assigning the right parameters into the model while skipping out of possible local minima. Thus two major questions were addressed in this study: What may be a good cost function for the estimation of the free parameters in a compartmental model? Is it possible to constrain a model containing nonhomogeneous distributions of ion channels once multiple dendritic and axonal recordings are added to the data set? To address these questions we used a single-neuron model with 19 free parameters and tested several cost functions using a well-defined genetic algorithm. The optimizations were performed on target data simulated with a known set of parameters, thus allowing us to directly assess the deviations of the best-fit parameter sets from the target parameter set under different conditions. We show that a cost function combining several modes of data analysis provides the best fit between the simulated data set and the model. Moreover, we show that increasing the data set by addition of dendritic and axonal recording sites improved the fit even further.
|
METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
The compartmental model and the various cost functions were programmed using NEURON 5.3 (Hines and Carnevale 1997), whereas the genetic algorithm was written in the Perl computer language. Simulations were run on a variety of computers running Linux on Pentium I, II, and III computers or SunOS 5.8 on SPARC workstations. We parallelized the process using a cluster of 15 to 32 machines sharing the same network file system (NFS), in which one of the machines functioned as a master—submitting and managing the jobs while the rest were used as slaves reading and writing their information from a shared directory in the NFS. Ion channel models were implemented using the NMODL extension of NEURON (Hines and Carnevale 2000). Results were analyzed using IgorPro 4.01 and a variety of Perl scripts.
Compartmental model
One model was used for all the simulations presented in this work. The model was constructed based on a previously published model of L5 neurons (Mainen et al. 1995) that was modified using recent experimental findings applying several assumptions:
THE MORPHOLOGY OF THE NEURON IS KNOWN.
This assumption relies on the excellent advances in morphological analysis in recent years allowing production of a good reconstruction of the neuron under study. Thus we used a simplified morphology (Fig. 1) containing the functional compartments of soma, axon, and dendrite instead of a full morphology. This model had 57 compartments consisting of a soma (length 20 µm, diameter 20 µm), an apical dendrite (length 800 µm, diameter 5 µm), and an axon that was modeled based on a previous axonal model suggested for L5 neurons (Mainen et al. 1995).
|
· cm2, Ra = 100 · cm, and Cm = 1.5 µF/cm2 with a passive reversal potential of Epassive = –65 mV. In the dendrite Cm was multiplied by a factor of two and Rm was divided by the same factor to account for dendritic spines (Stuart and Spruston 1998).
DENDRITIC CHANNEL GRADIENTS ARE A CONTINUOUS FUNCTION OF THE DISTANCE FROM THE SOMA.
In recent years many investigations have shown that the densities of dendritic conductances vary along the dendrite (Bekkers 2000b; Bekkers and Delaney 2001; Berger et al. 2001; Hoffman et al. 1997; Korngreen and Sakmann 2000; Magee 1999; Magee and Johnston 1995; Stuart and Häusser 1994; Williams and Stuart 2000). In all these cases the conductance gradient has been described by a continuous function of the distance along the somatodendritic axis. Thus we defined the dendritic gradient using three parameters: the somatic conductance density, the dendritic conductance density, and a distance parameter. The conductance changes linearly from the somatic value to the dendritic value over the distance along the dendrite set by the distance parameter. Distal to this distance from the soma the dendritic conductance is set to be constant and equal to the dendritic conductance density. This approach could be expanded to describe more complex channel gradients by defining several segments.
THE KINETICS OF THE CHANNELS IS KNOWN.
This assumption substantially reduces the number of free parameters, leaving mostly conductance density values to be constrained. We used the kinetics of conductances that have been experimentally characterized in sufficient detail to provide a kinetic model that could be implemented using the NMODL extension of the NEURON simulation language (Hines and Carnevale 2000). The model contained four active conductances: Kf, Ks, Ih, and Na. The two K+ conductances were based on recordings made from nucleated patches (Korngreen and Sakmann 2000). The voltage-gated Na+ conductance was also based on recordings from nucleated patches from L5 neocortical pyramidal neurons (A Korngreen and B Sakmann, unpublished results). The hyperpolarization activated nonselective cation conductance (Ih) was modeled based on cell-attached somatic and dendritic recordings (Berger et al. 2001; Williams and Stuart 2000). All the conductances were modeled using Hodgkin–Huxley-type models (Hodgkin and Huxley 1952).
Voltage-gated Na+ conductance
 |
 |
 |
 |
 |
). Therefore two additional parameters were used as target parameters for the genetic algorithm: the shift in the voltage of half activation (Nashift,act) and the shift in the voltage of half inactivation (Nashift,inact).
Fast inactivating K+ conductance (Kf)
 |
 |
 |
 |
 |
Slow inactivating K+ conductance (Ks)
 |
 |
 |
 |
 |
 |
 |
 |
Hyperpolarization activated cation conductance (Ih)
 |
 |
 |
The original set of parameters used in the model is listed in Table 1 along with the parameter ranges used by the genetic algorithm. The reversal potentials of the voltage-gated conductances were set to be 60, –80, and –30 mV for Na+, K+, and Ih, respectively. The same artificial data set was used for all the simulations consisting of six traces of membrane potential recordings induced by six different current injections (starting at –0.35 nA and increasing in steps of 0.15 nA).
|
A genetic algorithm is an optimization algorithm based on the mechanisms of Darwinian evolution, which uses random mutation, crossover, and selection operators to breed better models or solutions (individuals) from an originally random starting population (Mitchell 1996). In this study we started each minimization with a random population of 100 vectors, each describing a parameter set, and the model was evaluated for each one of them. The order of the parameters in the vectors was defined according to the role of the parameters in the model. Parameters that relate to the same neuronal function and can strongly influence each other, such as those that define a specific channel or passive ones, were placed close to each other. The population was sorted according to the fit score and a new generation was created. The creation of the new generation was done by point mutations applied to some of the values of the parameters sets (a point mutation will be issued with specific, defined probability and change a current value with a random one), and by a crossover operation in which two individuals were chosen randomly and exchange some of their values in the parameters vector. The best individual in the population was transferred to the next generation with no change, whereas the others were subjected to crossover (mixing parameters between sets) and mutation (randomly changing parameters inside a set) with a probability of 0.1 for both mutation and crossover. Those iterations continued until the termination criterion was met. Ideally, the termination criterion should be that the cost function reaches a value of zero. In practice this is not possible because the run time of the process is limited and reaching this score can take a substantial amount of time. In addition, the cost function will not reach a value of zero once noisy target data set will be used. In this work the process was terminated after 2,000 generations for all simulations. On average, a typical run of the genetic algorithm took about 1 wk to complete on our Linux cluster. Therefore only two runs were performed for each tested cost function or target data combination.
As mentioned earlier, simulations were carried out using two programming languages, Perl and NEURON. To narrow the parameter space, we bounded each parameter with a reasonable physiological value (Table 1). For purely technical reasons the Perl algorithm used variables bounded between 0 and 1 to represent each parameter as routinely performed when applying genetic algorithms to problems containing continuous parameters (Haupt and Haupt 1998). These values were linearly transformed to the real parameters of the compartmental model (bounded by the values given in Table 1) for each evaluation of the model. The analysis of the results was carried out using the real parameters and not those used internally by the Perl script.
Cost functions
Several functions were used to compare between the neuron's response to somatic current injection when the compartmental model was supplied with parameters varied by the genetic algorithm (test data set) and its response to current injection when the compartmental model was supplied with the fixed set of parameters shown in Table 1 (target data set).
INTERSPIKE INTERVAL (ISI).
This cost function computed the difference in the AP peak times for each trace of the simulation. A mean square calculation was made for each AP time
 | (1) |
WAVEFORM.
This cost function calculated the difference between the target and the test membrane potential traces and summed its mean squares yielding a cost function value that expresses the distance of the test from the target data set
 | (2) |
ISI + WAVEFORM.
This cost function was the simple sum of the above two functions (
2 + ISI2).
TRAJECTORY-DENSITY.
This cost function excluded the time domain by plotting the membrane potential V(t) versus its first time derivative 
(t), in which the periodicity of the signal was reflected by a closed loop that could be geometrically analyzed (LeMasson and Maex 2001). The phase plane was divided into small squares and the total number of points of the trajectory falling into each square was calculated. The cost is calculated by subtraction of two phase plane density matrices, one calculated from the test data (t) and one calculated from the target data set (T)
 | (3) |
COMBINED.
This cost function combined the Waveform analysis with the Trajectory-density and ISI. For the passive currents we used the Waveform calculations alone, whereas for the AP trains we also added those of the Trajectory-density and ISI; this sum was further penalized on the peak amplitude of the APs
 | (4) |
We used the normalized mean geometric distance of the test parameters from the target parameters to quantify the deviation of the parameter set obtained by the genetic algorithm
 | (5) |
|
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
The ability of genetic algorithms to converge to a global minimum is strongly influenced by the type of problem under investigation and the cost function used to quantify the distance between the target set of data and the data set generated by the genetic algorithm (Haupt and Haupt 1998; Mitchell 1996). Thus the first simulated experiment was designed to compare the performance of the genetic algorithm using five different cost functions (defined in METHODS). The target data set for these simulations was produced by somatic current injections and included six membrane potential traces recorded at the soma after different levels of current injection (Figs. 2 and 3, top). Previous optimizations of compartmental models using genetic algorithms have applied the ISI cost function (Eq. 1) to evaluate the quality of the fit (Eichler West 1996; Vanier and Bower 1999). Nevertheless, the test data set used in this study, simulated using the best parameter set obtained by the genetic algorithm using the ISI cost function, did not fit the target data set well (Fig. 2Ai) because, although the spiking pattern was faithfully reproduced, the test APs did not reach the full amplitude of the target APs and the passive traces, which are essentially meaningless in the ISI cost function's definition, thus missing the target data set curves (Fig. 2Ai).
|
|
; Vanier and Bower 1999), within a small number of generations, leveling off at an almost constant score for the rest of the algorithm's runtime (Fig. 2Aii). Furthermore, in both runs of the genetic algorithm the ISI cost function caused an overestimation of the values of GKs,node, GKs,soma, and GKf,node (Fig. 2Aiii). This overestimation was probably responsible for the small amplitude of the APs produced by the best parameter set (Fig. 2Ai). Additionally, the value of the axial resistance (Ra) had the largest error of the passive parameters, indicating a smaller passive space constant and thus considerable changes in the propagation of the AP into the axon and dendrite.
It has been predicted that the Waveform cost function (Eq. 2) may not be considered an appropriate cost function for algorithms attempting to constrain compartmental models (LeMasson and Maex 2001). Indeed, large deviations of both the passive responses and AP trains from the target traces (Fig. 2Bi) were observed when the Waveform cost function was used by the genetic algorithm to evaluate the quality of the fit. The most striking deviation was the convergence of the algorithm to a parameter set producing a test data set with a single AP, whereas the target data set expressed a train of five APs for the same current injection (Fig. 2Bi). The convergence of the best solution obtained by the genetic algorithm using the Waveform cost function also displayed a pattern of convergence different from that of the simulations using the ISI cost function (Fig. 2Bii); the two runs of the genetic algorithm produced a relatively small reduction in the fit value of the best individual as a function of the number of generations. Additionally, GNa,dend was overestimated. However, the minimizations using the Waveform cost function showed similar behavior to the runs using the ISI cost function in the overestimation of GKs,node, GKs,soma, GKf,node, and Ra (Fig. 2Biii).
To investigate whether it was possible to merge the different properties of the ISI and Waveform into one cost function, two runs of the genetic algorithm were performed using a simple sum of these two functions as one cost function (ISI + Waveform). As expected, using the ISI + Waveform cost function resulted in a different behavior. Similarly to the ISI cost function, the locations of the APs were faithfully constrained, whereas the passive traces were ill fitted (Fig. 2Ci). Additionally, the convergence of the best solution obtained by the genetic algorithm (Fig. 2Cii) displayed what appeared to be a combination of the convergence of the ISI (Fig. 2Aii) and Waveform functions (Fig. 2Bii); that is, the value of the best individual diminished considerably during the first 20 generations and then continued to decay at a slower rate until it leveled off after approximately 500 generations. Unlike the parameter sets obtained by separately using the ISI (Fig. 2Aiii) and Waveform (Fig. 2Biii) cost functions, the parameter set obtained using the ISI + Waveform did not display large deviations in parameters related to the AP axonal initiation zone (GKs,node and GKf,node, Fig. 2Ciii). However, the two runs of the genetic algorithm produced parameter sets displaying different degrees of variation of the same parameters. Thus although the final score of the genetic algorithm was similar in both runs of the genetic algorithm, the parameter vectors indicated that the two runs converged to two different subregions of parameter space.
To increase the weight of the APs in the target data set we used the Trajectory-density (Eq. 3) cost function (LeMasson and Maex 2001). Using this cost function the genetic algorithm constrained APs that had larger amplitudes than those observed with the previously tested cost functions (Fig. 3Ai). However, the APs were shifted relative to the APs in the target data set. Moreover, the convergence of the best solution obtained by the genetic algorithm using the Trajectory-density cost function was slower than that observed for the ISI and ISI + Waveform cost functions (Fig. 3Aii). Additionally, like the runs using the ISI and Waveform cost functions the value of GKs,node, was overestimated (Fig. 3Aiii); and, in one of the runs, the value of GNa,dend was overestimated as well. Another similarity to the minimizations using the ISI cost function is the overestimation of the value of Ra (Fig. 3Aiii). Finally, the ISI, Waveform, and Trajectory-density cost functions were assembled together to form the Combined cost function (Eq. 4). Using this new cost function, the genetic algorithm produced the best visual fit between the somatically recorded test and target data sets (Fig. 3Bi). This was manifested by the best solution obtained by the genetic algorithm, using the Combined cost function, converging slower than the previously tested cost functions (Fig. 3Bii). Furthermore, although the parameter set obtained by the two runs displayed large deviations in some of the parameters, these deviations were significantly smaller than those obtained with the other cost functions (Fig. 3Biii).
For further comparison of the performance of the different cost functions we calculated the mean geometric distance (Eq. 5) between the parameters obtained by the genetic algorithm and the target parameters (Fig. 4). The highest average deviation by best test parameters from the target parameters was observed when using the ISI cost function (129%) followed by the Waveform cost function (108%) and, finally, the ISI + Waveform cost function (90%). The best test parameter sets produced by the genetic algorithm using the Trajectory-density and Combined cost functions had similar mean deviations from the target parameters (75 and 73% respectively). To further investigate the different performance of the various cost functions, we tested the sensitivity of the three major cost functions (ISI, Waveform, and Trajectory-density) to variations in the parameter set. Using the target parameter set (Table 1) we generated 5,000 new parameter sets by randomly modifying all parameters within ±30% of their target value. The value of each cost function was calculated using these parameter vectors and the distribution of the results was displayed as logarithmic histograms (Fig. 5). We surmised that the shape of the results' histogram should manifest the sensitivity of a cost function to parameter variation because a sensitive cost function should generate a larger spectrum of results than a nonsensitive one. We found that the histogram of the results obtained using the ISI cost function displayed the smallest sensitivity to the variability of the parameter sets (Fig. 5A). We also found that the distribution of the results obtained using the same parameter vectors with the Waveform cost function (Fig. 5B) was broader than that obtained with the ISI function. However, the Waveform cost function displayed a sharp cutoff toward the low range of the score scale. Furthermore, the Trajectory-density cost function displayed a similar distribution of score results to that of the Waveform function, also producing results in the low range of the score scale (Fig. 5C). Thus similar to the ranking of the cost function using the mean geometric deviation (Fig. 4), the Trajectory-density function provided a greater sensitivity to changes in the parameter vector, in all probability allowing the genetic algorithm to better constrain the target parameter set.
|
|
Visual evaluations (Fig. 3) and numerical analyses of the obtained parameters (Figs. 4 and 5) suggested that the genetic algorithm using the Combined cost function was able to constrain the model better than the other cost functions tested in this study. We thus turned to investigate the second question raised in this study: Is it possible to constrain a model containing nonhomogeneous distributions of ion channels once multiple dendritic and axonal recordings are added to the data set? The results of the minimization performed using the Combined cost function and somatically recorded membrane potential traces (Fig. 3B) are shown again in Fig. 6A (to reduce the complexity of the figure only the membrane potential traces generated by 0.4- and –0.35-nA current injections are displayed). The membrane potential traces recorded in the initial segment of the axon (20 µm away from the soma) and at 400 µm along the dendrite, which were not used as part of the target data set, are shown in addition to the membrane potential recorded at the soma (Fig. 6A). Although the somatic test data set substantially fitted the target traces there were clear deviations of the traces at both the dendritic and axonal recording sites. However, an apparent visual improvement to the fit could be observed (Fig. 6B) when the target data set was increased by including the membrane potential traces recorded at 400 µm along the dendrite in addition to the somatic membrane potential traces. Although the fit of the dendritic and somatic data sets improved by this doubling of the target data set, the membrane potential traces at the axonal recording location still displayed visual deviations (Fig. 6B). Nevertheless, the addition of the dendritic recordings to the data set changed the deviations of the obtained parameter set from the target parameter set (Fig. 7A). Also, the error in some of the parameters decreased in comparison to the parameter set obtained when the target data set consisted of only somatic membrane potential traces. Surprisingly, the error in GNa,dend, a parameter predicted to be better constrained by the addition of dendritic data, increased. Notwithstanding, the mean geometric distance of the best test parameters obtained by the genetic algorithm from the target parameters for this simulation was slightly shorter than that obtained when the target data set contained only somatically recorded membrane potential traces (68 and 73% respectively, Fig. 7B).
|
|
|
DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
Cost function comparison
The genetic algorithm is an optimization method based on the principles of genetics and natural selection guiding the evolution of a population composed of many individuals to a state that maximizes the "fitness" of the population. Several criteria govern the selection of the cost function that estimates the fitness of the population found by the genetic algorithm. To begin with, the cost function should, as much as possible, be sensitive to changes in the parameters of the model; it is clear that a cost function that is insensitive to changes in one of the model's parameters will not enable a genetic algorithm, or any other optimization algorithm, to find an optimal value for this parameter. Subsequently, especially in minimization problems with interdependent parameters, the entire parameter set might miss the global minimum. Thus a rapid convergence or greediness of the best solution as the generations progress may be a hallmark of a poorly designed genetic algorithm or of a cost function that is not sensitive to the model's parameters (Haupt and Haupt 1998).
Of the five cost functions tested in this study the ISI cost function displayed the smallest degree of sensitivity to variations in the parameter set in the vicinity of the target parameter set (Fig. 5A). Indeed, the convergence of the genetic algorithm using the ISI cost function was rapid (Fig. 2Aii) and produced the largest mean deviation from the target parameter set (Fig. 4). Moreover, previous investigations using the ISI cost function in combination with a genetic minimization algorithm have reported similar convergence patterns (Eichler West 1996; Vanier and Bower 1999) and considerable deviations from the target parameters (Vanier and Bower 1999). Thus the ISI cost function can be considered an inappropriate cost function for the optimization of compartmental models. However, because the ISI analysis reduces the data from a large number of membrane potential points to a few timing points—forcing vital information to be lost—our findings are not surprising. It is also possible to predict that other ad hoc reductions of AP trains will not serve as good cost functions.
The Waveform cost function, the most widely used cost function for fitting a model to a set of data (Press et al. 1992), was found to be ill adapted for fitting AP trains (Fig. 2B). This is probably explained by the low applicability of the Waveform cost function to highly periodic signals (LeMasson and Maex 2001). Merging the ISI and Waveform functions did not combine the good qualities of the two functions either (Fig. 2C); it appears that the greediness of the ISI cost function guided the genetic algorithm into a subregion of parameter space by phase-locking the AP trains to the target data set, from which the tuning of the AP shape by the Waveform cost function was not powerful enough to enable the genetic algorithm to escape in the direction of the target parameter set.
By using the Trajectory-density cost function the genetic algorithm outperformed its combination with the ISI and Waveform cost functions (Fig. 3Ai). However, the APs simulated using the best test parameter set were shifted relative to the APs in the target data set. This temporal shift of the APs obtained by the genetic algorithm was probably a result of the missing time domain in the Trajectory-density cost function (Eq. 3). Moreover, the Trajectory-density cost function does not give significant weight to the passive traces arising from the relatively small changes in their first temporal derivative. It was thus clear that the Trajectory-density cost function overweighed the contribution of the APs in the target data set. Eventually, the combination of this property with the phase-locking ability of the ISI cost function and the Waveform cost function's ability to increase the weight of the passive responses provided the best visual fit (Fig. 3B) between the test and target data sets.
Several conclusions could be drawn from this part of the study. First of all, more research into the design of cost function for the minimization of compartmental models is called for; we were surprised to see how much influence the cost functions used for constraining compartmental models have over the behavior of the genetic minimization algorithm. Furthermore, all the functions tested in this study were off-the-shelf functions that had been used in similar studies (Eichler West 1996; LeMasson and Maex 2001; Vanier and Bower 1999) or were generally used to compare models with a data set (i.e., the Waveform cost function); that is, none of these functions was tailor made to fit the problem at hand. Compartmental models consist of a large set of coupled differential equations. This set of equations is being numerically integrated to obtain the spatial and temporal changes in the neuron's membrane potential. The integration process masks the dynamics of many parameters because it reduces the space defined by a large set of differential equations to a smaller set of values representing the membrane potential. This masking is probably the most intrinsic problem facing researchers attempting to constrain free parameters of a compartmental model. Consequently, further understanding of the loss of information between parameter spaces to the membrane potential is required to develop more efficient cost functions for compartmental models.
Second, the estimate of the robustness of a compartmental model relies heavily on the analysis method. In theory the best way to investigate the robustness of a compartmental model is to test all possible parameter combinations (Goldman et al. 2001). This technique rapidly becomes inapplicable—as a result of both computational reasons and the limitation of the human mind to grasp such a complex hypercube—as the number of free parameters increases. Alternatively, some information about the robustness can be gained by varying single or pairs of parameters (Bhalla and Bower 1993) or by randomly drawing parameter sets from one parent set (Golowasch et al. 2002) as we have performed here. For example, it can be argued that the analysis of the robustness of the model used here with the ISI cost function shows the model to be quite robust (Fig. 5A). However, an analysis using the Trajectory-density function alone may lead one to decide that the model is not robust (Fig. 5B).
Third, using only somatic recordings of the membrane potential is not a sufficient condition for constraining compartmental models containing nonhomogeneous distributions of ion channels. Many of the minimizations carried out using only somatic target data set (Figs. 2 and 3) show considerable deviations in the values of parameters related to the axonal initiation zone in some cases and to the dendrite in others. Additionally, even under the best optimization, obtained using the Combined cost function, although a fit in substantial agreement with the somatic data was found, large deviations from the dendritic and axonal membrane potential traces, not included in the target data set, could be observed (Fig. 6A). To reiterate this point, our simulations strongly indicate that relying on visual comparison between the target and test data sets may lead the naïve investigator to select a parameter set that may not be related to the parameter set describing the activity of the neuron.
Target data set selection
In this study we used a fixed data set of six somatic current injections to investigate the effect of various cost functions on the convergence of the model. It is clear that data sets containing more somatic current injections will provide the genetic algorithm with more information about the activity of the ion channels in the model. However, the improvement to the fit will probably be sublinear to the number of additional somatic sweeps. In addition, it is possible to extend the data set further by using less-conventional current injection protocols, such as current ramps or other protocols, thus sampling a larger portion of the activation range of the voltage-gated channels in the neuron. Moreover, although initially one full run of the genetic algorithm lasted approximately 1 wk, the duration required for the algorithm to converge appeared to increase linearly with the number of recorded sweeps and recording locations.
Addition of dendritic membrane potential traces to the target data improved the visual fit (Fig. 6B) and the fit of many of the free parameters (Fig. 7A). However, this spoiled the previous good fit of the dendritic conductance density of the voltage-gated Na+ conductance (GNa,dend). This was probably caused by the fact that, much like attempting to cover two people with one blanket, where pulling on one end to cover one person leaves a smaller piece of the blanket to the other person, adding a dendritic recording reduces the weight of the somatic data traces, forcing the membrane potential at the dendritic recording site to fit the target set of data. It was also most likely that because of the underestimation of the AP amplitude at the soma (as a result of the overestimation of GKs,soma and GKf,node) that the back-propagating AP required additional boosting from voltage-gated Na+ conductances. It was only the addition of membrane potential recordings from the axon (Fig. 6C) that reduced the overestimation of the axonal conductance densities, ultimately leading to a reduction in the error in the dendritic sodium conductance density. Furthermore, it is important to note that, because visual inspections of the somatic fit do not reveal major differences between the fit obtained only with a somatic data set and those carried out also with dendritic and axonal membrane potential traces, a naïve researcher fitting only somatic data may obtain a poor estimation of dendritic and/or axonal parameters. Finally, we conclude that, at least in the case of the model presented in this work, constraining compartmental models containing nonhomogeneous distributions of ion channels requires a target data set containing membrane potential recordings from several locations along the neuron.
Experimental and practical consequences
Based on the model and simulations presented in this work it is possible to suggest an experimental scheme. By blocking voltage-gated calcium channels it is possible to appreciably reduce the number of active channels in the membrane of most neurons. Such an experiment will eliminate from the parameter space the conductance densities of voltage-gated calcium channels and that of the calcium-activated potassium channels. The remaining channels are the voltage-gated sodium and potassium channels and the hyperpolarization activated cationic conductance Ih. Our simulations indicate that under these conditions it may be possible to constrain a model for such neurons, providing that the membrane potential will be measured at several locations along the dendritic and axonal trees. Dendritic and axonal measurements of the back-propagating AP have been performed for over a decade (Stuart et al. 1999). However, simultaneous recordings of the membrane potential from five locations, including axonal recordings, using patch electrodes have yet to be performed. Alternatively, recording the membrane potential from many locations may be performed using voltage-sensitive dyes (Antic 2003; Antic et al. 1999). It is also possible to envision even more complex experimental scenarios in which, after the recording of the membrane potential, some of the channel blockers will be omitted from the bathing medium and new recordings of the APs will be made. This snapshot of neuronal activity could then be used to constrain a more complex compartmental model using the model constrained for fewer channels with fixed parameters. Thus it may be possible to progressively constrain more complex models, eventually reaching the desired target of a full model of the investigated neuron.
|
GRANTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
ACKNOWLEDGMENTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
FOOTNOTES |
|---|
Address for reprint requests and other correspondence: A. Korngreen, Faculty of Life Sciences, Bar-Ilan University, Ramat-Gan 52900, Israel (E-mail: korngra{at}mail.biu.ac.il)
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
Antic SD, Major G, and Zecevic D. Fast optical recordings of membrane potential changes from dendrites of pyramidal neurons. J Neurophysiol 82: 1615–1621, 1999.
Antic SD. Action potentials in basal and oblique dendrites of rat neocortical pyramidal neurons. J Physiol 550: 35–50, 2003.
Ariav G, Polsky A, and Schiller J. Submillisecond precision of the input–output transformation function mediated by fast sodium dendritic spikes in basal dendrites of CA1 pyramidal neurons. J Neurosci 23: 7750–7758, 2003.
Bekkers JM. Properties of voltage-gated potassium currents in nucleated patches from large layer 5 cortical pyramidal neurons of the rat. J Physiol 525: 593–609, 2000a.
Bekkers JM. Distribution and activation of voltage-gated potassium channels in cell-attached and outside-out patches from large layer 5 cortical pyramidal neurons of the rat. J Physiol 523: 611–620, 2000b.
Bekkers JM and Delaney AJ. Modulation of excitability by alpha-dendrotoxin-sensitive potassium channels in neocortical pyramidal neurons. J Neurosci 21: 6553–6560, 2001.
Berger T, Larkum ME, and Lüscher HR. High Ih channel density in the distal apical dendrite of layer V pyramidal cells increases bidirectional attenuation of EPSPs. J Neurophysiol 85: 855–868, 2001.
Bhalla US and Bower JM. Exploring parameter space in detailed single neuron models: simulations of the mitral and granule cells of the olfactory bulb. J Neurophysiol 69: 1948–1965, 1993.
Bischofberger J and Jonas P. Action potential propagation into the presynaptic dendrites of rat mitral cells. J Physiol 504: 359–365, 1997.[CrossRef][ISI][Medline]
Brown AM, Sayer RJ, Schwindt PC, and Crill WE. P-type calcium channels in rat neocortical neurones. J Physiol 475: 197–205, 1994a.
Brown AM, Schwindt PC, and Crill WE. Voltage dependence and activation kinetics of pharmacologically defined components of the high-threshold calcium current in rat neocortical neurons. J Neurophysiol 70: 1530–1543, 1993.
Brown AM, Schwindt PC, and Crill WE. Different voltage dependence of transient and persistent Na+ currents is compatible with modal-gating hypothesis for sodium channels. J Neurophysiol 71: 2562–2565, 1994b.
Chen WR, Midtgaard J, and Shepherd GM. Forward and backward propagation of dendritic impulses and their synaptic control in mitral cells. Science 278: 463–467, 1997.
Chen WR, Shen GY, Shepherd GM, Hines ML, and Midtgaard J. Multiple modes of action potential initiation and propagation in mitral cell primary dendrite. J Neurophysiol 88: 2755–2764, 2002.
Colbert CM and Pan E. Ion channel properties underlying axonal action potential initiation in pyramidal neurons. Nat Neurosci 5: 533–538, 2002.
