## Abstract

Compartmental models with many nonlinearly and nonhomogeneous distributions of voltage-gated conductances are routinely used to investigate the physiology of complex neurons. However, the number of loosely constrained parameters makes manually constructing the desired model a daunting if not impossible task. Recently, progress has been made using automated parameter search methods, such as genetic algorithms (GAs). However, these methods have been applied to somatically recorded action potentials using relatively simple target functions. Using a genetic minimization algorithm and a reduced compartmental model based on a previously published model of layer 5 neocortical pyramidal neurons we compared the efficacy of five cost functions (based on the waveform of the membrane potential, the interspike interval, trajectory density, and their combinations) to constrain the model. When the model was constrained using somatic recordings only, a combined cost function was found to be the most effective. This combined cost function was then applied to investigate the contribution of dendritic and axonal recordings to the ability of the GA to constrain the model. The more recording locations from the dendrite and the axon that were added to the data set the better was the genetic minimization algorithm able to constrain the compartmental model. Based on these simulations we propose an experimental scheme that, in combination with a genetic minimization algorithm, may be used to constrain compartmental models of neurons.

## INTRODUCTION

The majority of synapses in the CNS connect to dendrites. These dendrites transform information received from synapses into a code that is then translated by the axon to action potentials (APs) that are transmitted to other neurons. The properties and functions of dendrites in the CNS have been intensively studied especially over the past decade (for reviews see Johnston 1999; 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 Ca^{2+} 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

### Simulation environment

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

##### THE PASSIVE MEMBRANE PARAMETERS ARE HOMOGENEOUS.

In this study the passive parameters were set to *R*_{m} = 21,000 Ω · cm^{2}, *R*_{a} = 100 Ω · cm, and *C*_{m} = 1.5 μF/cm^{2} with a passive reversal potential of *E*_{passive} = −65 mV. In the dendrite *C*_{m} was multiplied by a factor of two and *R*_{m} 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: *K*_{f}, *K*_{s}, *I*_{h}, 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 (*I*_{h}) 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 where the symbols have their usual meaning in the Hodgkin–Huxley formalism. Four parameters describing the spatial distribution of this conductance were defined as target parameters for the genetic algorithm: the somatic conductance density (*G*_{Na,soma}), the nodal conductance density (*G*_{Na,node}), the dendritic conductance density (*G*_{Na,dend}), and the distance along the dendrite over which the conductance density changes from its somatic value to the dendritic value (Na_{dist}). The relative distribution of this conductance is shown in Fig. 1. It was reported that the activation and inactivation curves of the voltage-gated Na^{+} conductance were shifted in the axon (Colbert and Pan 2002). Therefore two additional parameters were used as target parameters for the genetic algorithm: the shift in the voltage of half activation (Na_{shift,act}) and the shift in the voltage of half inactivation (Na_{shift,inact}).

Fast inactivating K^{+} conductance (*K*_{f}) Four parameters describing the distribution of this conductance were defined as target parameters for the genetic algorithm: the somatic conductance density (*G _{K}*

_{f,soma}), the nodal conductance density (

*G*

_{K}_{f,node}), the dendritic conductance density (

*G*

_{K}_{f,dend}), and the distance along the dendrite over which the conductance density changes from its somatic value to the dendritic value (

*K*

_{f,dist}). The relative distribution of this conductance is shown in Fig. 1.

Slow inactivating K^{+} conductance (*K*_{s}) Four parameters describing the distribution of this conductance were defined as target parameters for the genetic algorithm: the somatic conductance density (*G _{K}*

_{s,soma}), the nodal conductance density (

*G*

_{K}_{s,node}), the dendritic conductance density (

*G*

_{K}_{s,dend}) and the distance along the dendrite over which the conductance density changes from its somatic value to the dendritic value (

*K*

_{s,dist}). The relative distribution of this conductance is shown in Fig. 1.

Hyperpolarization activated cation conductance (*I*_{h}) This conductance was inserted only in the somatic and dendritic compartments, as shown in Fig. 1, assuming a homogeneous distribution. Thus only one parameter was used to represent it in the genetic minimization scheme—the conductance density (*G _{I}*

_{h}).

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 *I*_{h}, 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).

### Genetic algorithm

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) where *T* represents the target data set and *t* the spiking times in each membrane potential trace of the test data set. *Z* is the total number of APs in each trace and *M* the number of membrane potential sweeps simulated in the model.

##### 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) where *T* represents the target data set and *t* the membrane potential changes in the test data set. *N* is the total number of points in each membrane potential trace and *M* the number of membrane potential sweeps simulated in the model.

##### ISI + WAVEFORM.

This cost function was the simple sum of the above two functions (ψ^{2} + *ISI*^{2}).

##### TRAJECTORY-DENSITY.

This cost function excluded the time domain by plotting the membrane potential *V*_{(t)} versus its first time derivative *V̇ _{(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) where

*X*and

*Y*are the number of elements in the phase plane matrix and

*M*is the total number of membrane potential sweeps simulated in the model.

##### 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) where *O*_{1} is the penalty value given for the difference between the target and test AP amplitude to encourage the algorithm to favor parameter sets that produce APs with large amplitudes.

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) where *S* represents the total deviation score from the target parameters, *N* is the number of free parameters in the simulation, *P _{s}*(

*i*) is the best-fit value of the

*i*th parameter, and

*P*(

_{D}*i*) is the target value.

## RESULTS

### Cost function comparison

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. 2*Ai*) 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. 2*Ai*).

The value of the best individual in each generation calculated by the *ISI* cost function converged, in agreement with previously published investigations (Eichler West 1996; 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. 2*Aii*). Furthermore, in both runs of the genetic algorithm the *ISI* cost function caused an overestimation of the values of *G _{K}*

_{s,node},

*G*

_{K}_{s,soma}, and

*G*

_{K}_{f,node}(Fig. 2

*Aiii*). This overestimation was probably responsible for the small amplitude of the APs produced by the best parameter set (Fig. 2

*Ai*). Additionally, the value of the axial resistance (

*R*

_{a}) 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. 2*Bi*) 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. 2*Bi*). 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. 2*Bii*); 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, *G*_{Na,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 *G _{K}*

_{s,node},

*G*

_{K}_{s,}

_{soma},

*G*

_{K}_{f,node}, and

*R*

_{a}(Fig. 2

*Biii*).

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. 2*Ci*). Additionally, the convergence of the best solution obtained by the genetic algorithm (Fig. 2*Cii*) displayed what appeared to be a combination of the convergence of the *ISI* (Fig. 2*Aii*) and *Waveform* functions (Fig. 2*Bii*); 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. 2*Aiii*) and *Waveform* (Fig. 2*Biii*) cost functions, the parameter set obtained using the *ISI* + *Waveform* did not display large deviations in parameters related to the AP axonal initiation zone (*G _{K}*

_{s,node}and

*G*

_{K}_{f,node}, Fig. 2

*Ciii*). 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. 3*Ai*). 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. 3*Aii*). Additionally, like the runs using the *ISI* and *Waveform* cost functions the value of *G _{K}*

_{s,node}, was overestimated (Fig. 3

*Aiii*); and, in one of the runs, the value of

*G*

_{Na,dend}was overestimated as well. Another similarity to the minimizations using the

*ISI*cost function is the overestimation of the value of

*R*

_{a}(Fig. 3

*Aiii*). 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. 3

*Bi*). 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. 3

*Bii*). 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. 3

*Biii*).

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. 5*A*). We also found that the distribution of the results obtained using the same parameter vectors with the *Waveform* cost function (Fig. 5*B*) 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. 5*C*). 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.

### Recordings from several dendritic and axonal locations

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. 3*B*) are shown again in Fig. 6*A* (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. 6*A*). 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. 6*B*) 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. 6*B*). 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. 7*A*). 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 *G*_{Na,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. 7*B*).

When another dendritic recording of the membrane potential at 200 μm along the dendrite was added to the target data set (i.e., soma and two recording locations on the apical dendrite at 400 and 200 μm), an improvement in the result was not observed (data not shown in Fig. 6 but see analysis in Fig. 7). Still, visual improvement in the fits at the soma, dendrite, and axon could be observed (Fig. 6*C*) when a recording of the membrane potential from the axon's initial segment, 20 μm away from the soma, was added to the target data set (four recording locations in total). This addition of data from the axonal recording site reduced the error in most of the parameters and especially in the overestimated *G _{K}*

_{s,soma},

*G*

_{K}_{f,node}, and

*G*

_{Na,dend}(Fig. 7

*A*). The mean geometric distance of the parameters obtained by the genetic algorithm from the target parameters for this simulation was shorter than that obtained when the target data set contained only somatically recorded membrane potential traces (51 and 73% respectively, Fig. 7

*B*). Similar results were obtained when the target data set contained recordings from the soma, axon, and one dendritic location 400 μm away from the soma (data not shown). Nonetheless, increasing the target data set further by an additional measurement from 600 μm along the dendrite did not change the quality of the fit compared with the four locations (data not shown in Fig. 6 but see analysis in Fig. 7).

## DISCUSSION

Every optimization of a model using numerical data relies on three basic conditions: the cost function, the quality of the target data set, and the minimization algorithm. In this study we selected the genetic algorithm as the minimization algorithm and tested several cost functions and data combinations in an attempt to determine the answer to two questions: 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? We show that a cost function combining analyses of many facets of passive and active responses of the model enables the genetic algorithm to produce the best set of parameters (Figs. 2 and 3). Furthermore, addition of membrane potential recordings from the dendrite and axon allows the genetic algorithm to progressively find a better parameter set (Figs. 6 and 7). Our results thus suggest that combining multiple membrane potential recordings with genetic minimization algorithms may be used to constrain functional models for multicompartment neurons containing nonhomogeneous distributions of ion channels.

### 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. 5*A*). Indeed, the convergence of the genetic algorithm using the *ISI* cost function was rapid (Fig. 2*Aii*) 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. 2*B*). 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. 2*C*); 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. 3*Ai*). 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. 3*B*) 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. 5*B*).

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. 6*A*). 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. 6*B*) and the fit of many of the free parameters (Fig. 7*A*). However, this spoiled the previous good fit of the dendritic conductance density of the voltage-gated Na^{+} conductance (*G*_{Na,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 *G _{K}*

_{s,soma}and

*G*

_{K}_{f,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. 6

*C*) 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 *I*_{h}. 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

This research was supported by Israel Science Foundation Grant 345/04.

## Acknowledgments

The authors thank Drs. M. Helmstaedter, A. Roth, Y. Louzoun, R. Unger, M. Gurkiewicz, D. Bar-Yehuda, and A. Schaefer for commenting on the manuscript.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2005 by the American Physiological Society