Automated optimization of a reduced layer 5 pyramidal cell model based on experimental data

The construction of compartmental models of neurons involves tuning a set of parameters to make the model neuron behave as realistically as possible. While the parameter space of single-compartment models or other simple models can be exhaustively searched, the introduction of dendritic geometry causes the number of parameters to balloon. As parameter tuning is a daunting and time-consuming task when performed manually, reliable methods for automatically optimizing compartmental models are desperately needed, as only optimized models can capture the behavior of real neurons. Here we present a three-step strategy to automatically build reduced models of layer 5 pyramidal neurons that closely reproduce experimental data. First, we reduce the pattern of dendritic branches of a detailed model to a set of equivalent primary dendrites. Second, the ion channel densities are estimated using a multi-objective optimization strategy to fit the voltage trace recorded under two conditions - with and without the apical dendrite occluded by pinching. Finally, we tune dendritic calcium channel parameters to model the initiation of dendritic calcium spikes and the coupling between soma and dendrite. More generally, this new method can be applied to construct families of models of different neuron types, with applications ranging from the study of information processing in single neurons to realistic simulations of large-scale network dynamics.


Introduction
To incorporate realism into large-scale simulations of cortical and other networks (Traub et al., 2005;Markram, 2006), one needs to construct biophysically realistic compartmental models of the individual neurons in the circuit. Many parameters of these models have not been directly measured experimentally; therefore, these parameters must be tuned to match the experimentally observed input-output relation of the neuron. Solving the resulting nonlinear optimization problem is difficult and requires extensive computing resources, especially for models comprising a detailed neuronal morphology and a large number of compartments (Traub et al., 2005;Achard and De Schutter, 2006;Markram, 2006;Druckmann et al., 2007;Hay et al., 2011).
Here we develop reduced models of neocortical layer 5 pyramidal cells with a small number of compartments to represent the dendritic geometry. Compared to fully detailed compartmental models, reduced models confer significant speed advantages both for the optimization of the single neuron model and for simulation of networks of such neurons. The reduced model's geometry, even though simplified, should still incorporate the fact that synaptic inputs arriving at different layers in the dendritic tree are integrated differently (Larkum et al., 1999(Larkum et al., , 2004Schaefer et al., 2003;Branco et al., 2010). This implies that a sufficient number of compartments must be used for the reduced dendritic morphology. The ability of a neuron to lock onto fast fluctuations in the input depends critically on how sharp the action potential onset is in its voltage time-course (Naundorf et al., 2005;Palmer and Stuart, 2006;Kole et al., 2007;Popovic et al., 2011). As the initiation site of the action potential, located in the axon (Stuart and Sakmann, 1994;Palmer and Stuart, 2006;Kole et al., 2007;Popovic et al., 2011), determines how steep the action potential onset is (Yu et al., 2008), a reduced model must also have a minimum number of compartments for the axon.
We construct the reduced pyramidal cell models step by step, applying methods adapted to the problem (Roth and Bahl, 2009). The first step, which determines the reduced cell morphology and passive membrane properties, follows the tradition of defining certain passive electrical properties of the full model that are preserved in the reduced model (Stratford et al., 1989;Bush and Sejnowski, 1993;Destexhe, 2001). Second, we apply a powerful optimization strategy, Evolutionary Multi-Objective Optimization (EMOO) (Deb et al., 2002;Druckmann et al., 2007). This method, starting from a family of models characterized by multiple features, generates a suite of new models at each step, without making an a priori determination as to which one of the multiple desired objectives is most important.
In previous approaches to evolutionary model optimization, up to twenty different features are defined (Hay et al., 2011), from the action potential width to the depth of the after-hyperpolarization potential (AHP). Each feature is associated with a single number. In contrast, we return to the classical approach of minimizing the least square difference between the recorded trace and the model's response. We define four least square difference objective functions, one for each of the different time scales present in the dynamics -from the AP onset dynamics on the microsecond time scale to the slow sub-threshold charging phase lasting several tens of milliseconds.
The experimental data we fit include the responses of a layer 5 pyramidal neuron to somatic current injection when the apical dendrite is occluded or pinched (Bekkers and Häusser, 2007). Such a procedure yields essential information about the electrical properties of the cell and its apical dendrite. This information allows parameter optimization to narrow down the possible combinations of channel densities and properties along the dendritic geometry that could explain the neuron's voltage response to somatic current injection.
Finally, we show how the evolutionary approach can be extended to ensure that a neuronal model captures features that do not take on continuous values, but are discrete. For this purpose, we took another data set from a different experiment on older layer 5 pyramidal neurons. We optimized the dendritic calcium channel parameters and adjusted these values such that the model reproduces the shape of the dendritic calcium AP as well as somato-dendritic coupling factors found in experiments (Schaefer et al., 2003).
After pursuing these three steps in optimizing neuronal models, we present a family of 10 reduced models of layer 5 pyramidal neurons whose input-output relation matches a range of experimental data. These models could be used in large-scale network simulations of the neocortex.

The cell model
Our aim is to create a reduced pyramidal cell model that is simple and fast but detailed enough to show complex somato-dendritic interactions. The model is based on standard techniques from compartmental and ion channel modeling (Hodgkin and Huxley, 1952;Rall, 1962), implemented in NEURON 7.1 (Carnevale and Hines, 2005) and controlled via the NEURON-Python interface (Hines et al., 2009).
To obtain a simplified geometry of the dendrites, we model the functional neuronal sections (soma, basal dendrites, apical dendrite and the apical dendritic tuft) each by a single cylinder whose length and diameter will be later determined by the optimization algorithm we describe. The axonal geometry is based on a detailed reconstruction (Zhu, 2000) and consists of a conical axon hillock (l = 20 m) which has a diameter of 3.5 m at the soma connection and tapers to 2.0 m. The conical axon initial segment (iseg; l = 25 m) is connected to the hillock and its diameter tapers from 2.0 m to 1.5 m. The actual axon (l = 500 m) is connected to the initial segment and has a uniform diameter of 1.5 m. We did not model nodes of Ranvier or myelination. As the reduced model should be fast, the number of compartments ought to be as small as possible. We chose the following compartment numbers for the functional sections: soma = 1; basal dendrite = 1; apical dendrite = 5; apical dendritic tuft = 2; axon hillock = 5; initial segment = 5; axon = 1. Hence, the model has a total of 20 compartments.
Ion channels were selected and distributed based on recent experimental findings and modeling studies and were downloaded from ModelDB (Hines et al., 2004): A hyperpolarization-activated cation channel (HCN) (Kole et al., 2006) was inserted into the basal dendrite, the apical dendrite and the dendritic tuft. A transient sodium channel (Nat) (Kole et al., 2006) was placed into the soma, the axon hillock, the initial segment, the apical dendrite and the dendritic tuft. The voltage dependency of the channel kinetics was shifted to higher values (vshift Nat = +10 mV) in all compartments to yield higher thresholds (Mainen and Sejnowski, 1996). We introduced a second voltage shift (vshift2 Nat ) for the Nat channel in the initial segment to account for different channel properties in this area (Colbert and Pan, 2002). Nat channel density decayed linearly in the apical dendrite with distance from the soma (Mainen et al., 1995;Keren et al., 2009). A fast potassium channel (Kfast) (Kole et al., 2006) was inserted into the soma, the apical dendrite and the tuft. Its density decayed exponentially from the soma towards the tuft (Keren et al., 2009). A slow potassium channel (Kslow) was inserted into the soma, the apical dendrite and tuft and channel densities decayed exponentially with distance from the soma (Korngreen and Sakmann, 2000). A persistent sodium channel (Nap) was inserted into the soma to adjust the neuron's excitability (Traub et al., 2003). A muscarinic potassium channel (Km) was inserted into the soma. The Km channel is a non-inactivating voltage-dependent slow potassium channel which is thought to play a role in spike frequency adaption (Winograd et al., 2008). A slow calcium channel (Cas) was inserted into the tuft. The voltage dependency of its kinetics could be shifted (vshift Cas ) to adjust activation thresholds. Finally, a calcium dependent potassium channel (KCa) (Mainen and Sejnowski, 1996) as well as a calcium pump (CP) (Kole et al., 2006) was inserted into the tuft. The 10 optimized reduced model neurons and the ion channel models are available for download at http://senselab.med.yale.edu/ModelDB.

Free model parameters
Reliable data on the biophysical properties of ion channels in pyramidal neurons are rare and measurements mostly stem from different neurons from different animals or even species. Moreover, due to experimental limitations many parameters just cannot be measured. In particular, information about ion channel densities and kinetics in distal dendritic branches are not available. It is, therefore, not sufficient to put together all existing information on pyramidal neurons and build a working model; a long list of uncertainties remains. We made a selection of the most uncertain parameters that we thought could be estimated by means of an optimization strategy: The lengths (l), the diameters (d) and the axial resistances (Ra) of the functional sections of the reduced morphology were free parameters, as were the values of the membrane resistance (R m ) and capacitance (C m ). We introduced a dendritic scaling factor (dendscaling) that allowed the adjustment of dendritic membrane resistance and capacitance relative to the soma. This should account for dendritic spines, systematic errors in dendritic reconstructions or different ratios between the soma size and the dendritic tree size in different cells. Further free parameters were the ion channel densities in different compartments, which were described by the channel's maximal ionic conductance per membrane area (e.g. somaḡ Nat ). It was shown that certain ion channel densities along the apical dendrite can be described by simple Table 1 Full list of all optimized parameters with their lower and upper search bounds (LSB, USB) together with a summary of model properties for all 10 reduced models. (a) The morphological parameters found by reducing the detailed model (first step) are shown. The remaining parameters are given because the method maintains the total surface area of the functional sections during the optimization (Asoma = 1682 m 2 , A basal = 7060 m 2 , A apical = 9312 m 2 and A tuft = 9434 m 2 and therefore d soma = 23 m, L soma = 23 m, d basal = 8.7 m, d apical = 5.9 m, d tuft = 6.0 m). (b) The model parameters are shown that were found after the optimization of ionic conductances (second step). Models 1-6 were optimized using target data before and during pinching, while models 7-10 were optimized only with target data from the intact neuron. (c) Model parameters found after optimizing the dendritic calcium dynamics (third step) are shown. In models 4, 5, 9 and 10 no parameters could by found satisfying the requirements for proper calcium dynamics. (d) We summarized a few properties of the 10 resulting models: the ratio of the somatic sodium and potassium channel density, the ratio of the sodium channel density in the initial segment and soma as well as the ratio of the tuft HCN and basal HCN channel density. It is also summarized that all models show back-propagating APs and which models show a decay of BAP amplitude with somatic distance (illustrated in Fig. 7). Finally we summarize the modulation of the resting potential along the apical dendrite (distal voltage-somatic voltage at rest, see Fig. 6). If proper calcium dynamics could be found in the third step we show the resulting coupling factor.  (Berger et al., 2001;Kole et al., 2006). We assumed an exponential decay of the potassium channel density along the apical dendrite but allowed the space constant ( Kfast , Kslow ) to be a free parameter. Nat and HCN channel densities were changed linearly along the apical dendrite with slopes calculated based on the channels' densities in the soma and tuft. Finally, some kinetic parameters of the same ion channel type are distinct in different parts of the neuron (Colbert and Pan, 2002). Therefore we chose the voltage dependency of the Nat channel in the axon initial segment (iseg vshift2 Nat ) and the voltage dependency of the Cas-channel in the tuft (tuft vshift Cas ) as further free parameters. The free model parameters used for optimization as well as lower and upper search bounds are listed in Table 1.

The optimization algorithm
In all three optimization steps we used the Evolutionary Multi-Objective Optimization (EMOO) algorithm (Deb, 2001;Deb et al., 2002). The EMOO-algorithm allows one to simultaneously minimize multiple and possibly conflicting error functions and is therefore especially well suited for the optimization of spiking neuron models to experimental data (Druckmann et al., 2007(Druckmann et al., , 2008. EMOO is a genetic algorithm and uses mechanisms inspired by biological evolution, such as selection, crossover and mutation, to grow a population of size N to a certain capacity C and to transfer good individuals into the next generation. We used a Simulated Binary Crossover operator (Deb and Agrawal, 1995) and a Polynomial Mutation operator (Deb, 2001). The efficacy of an operator is controlled by Á c and Á m respectively. Large values Á c and Á m mean that the operator's effect is weak. The EMOO-algorithm was implemented in Python and run parallelized on an AMD x86 64 cluster with 80 processors running Linux. The Python-Framework for the evolutionary multi-objective optimization is available for download under www.g-node.org/emoo.

Stepwise fitting strategy
An important aspect of a good model is that it closely reproduces data used to construct the model and, even more importantly, predicts key features of new data that were not used to constrain the model (Druckmann et al., 2011). Tuning a set of parameters by hand such that the model fulfills these requirements is a daunting, timeconsuming or even impossible task and a new set of data or minor model modifications might require a repetition of that process. Therefore the process of creating a model and tuning its parameters should be as automated as possible. To divide the parameter space into smaller units and to optimize each parameter subset independently we have developed a three-step fitting strategy: In the first step we estimated the geometry of the model, including the axial resistances of the sections. The geometrical parameters of the model were fixed after this step. In the second step all ion channel and membrane parameters affecting somatic spiking but not calcium spiking dynamics were optimized. Once this step was carried out, these parameters were also fixed. Finally, in the third step we estimated the parameters needed for dendritic calcium spike dynamics. In an optimal scenario all these data would come from the same neuron from the same animal. This would be an experimental challenge and currently such a set of data does not exist. Therefore we have to combine data from different experiments in each of these steps.

First step: optimizing the reduced geometry
To obtain a representative geometry for the reduced model we started with a detailed model reconstruction of a layer 5 pyramidal neuron from a young (≈P21) rat (Stuart and Spruston, 1998). To do so we adopted a model simplification strategy (Destexhe, 2001): As we are only optimizing the neuronal geometry we removed all active conductances and globally set the membrane resistance (R m ) to 15,000 cm 2 , the membrane capacitance (C m ) to 1 pF/m 2 , the reversal potential (e pas ) to −70 mV in both the detailed and in the reduced model. The specific membrane parameters (R m , C m , e pas ) do not change with geometry and therefore do not need to be optimized in the first step. However they will need further tuning when the model is matched to spiking data in the second step. Next, we assigned each compartment in the detailed model to one of the functional sections (basal dendrite, soma, apical dendrite, tuft), determined their surface areas and set the size of the functional sections in the reduced model to the same values (A soma = 1682 m 2 , A basal = 7060 m 2 , A apical = 9312 m 2 and A tuft = 9434 m 2 ). Then, in the detailed model, we set the axial resistance (Ra) globally to the commonly used value of 100 cm. As it is not clear what the axial resistance and the length of the functional sections in the reduced model should be, we needed to estimate these values. Once we know the length of a functional section, we can also compute its diameter, as this should lead to the same surface area as measured in the reconstruction. Destexhe (2001) only used the steady state voltage in response to constant current injections to fit these parameters across all compartments. We used this target function, as well, but extended the method to reproduce the detailed neuron's somatic input impedance and phase shift functions for oscillatory somatic input currents (0-1000 Hz). These impedance functions were shown to contain information about the sub-threshold neuronal dynamics of the whole morphology (Fox, 1985;Borst and Haag, 1996). We determined these three functions in the detailed model and, for a given parameter combination, in the reduced model and calculated the sum of squared differences. This gave us three features to minimize by EMOO. We used a population size of N = 350 and a capacity of C = 700 individuals and evaluated 100 generations. The crossover parameter started at Á c = 5 and increased linearly to Á c = 50, while the mutation parameter increased from Á m = 10 to Á m = 500, thereby reducing the strength of these operators during evolution. The mutation probability per parameter was constant at 10%. We dropped the axon in this step as the detailed reconstruction lacks an axon. The axon was appended again afterwards and its membrane properties were chosen to equal those found for the soma. Appending the axon reduced the somatic input resistance from 69 M to 63 M .
After optimization the passive response properties of the reduced model matched those of the detailed model ( Fig. 1). One set of optimal parameters is given in Table 1; further optimization trials have led to different parameter combinations (not shown) that reproduced the passive response properties of the detailed model similarly well. To challenge the optimization method we injected the same Gaussian noise into the soma of the reduced and complex model and measured the resulting voltage responses in the soma and in the dendrite. As demonstrated in Fig. 2, the passive voltage responses in both models are almost indistinguishable.

Second step: optimizing the ion channel parameters
In the next step we estimated the ion channel parameters affecting somatic spiking. This was done using experimental data on the somatic spiking dynamics under two different conditions, first in the intact neuron, and second while the apical dendrite has been occluded using a method called Pinching (Bekkers and Häusser, 2007). Dendritic calcium spikes develop more fully in older pyramidal neurons  and are activated by strong dendritic local depolarization or by somatic input combined with dendritic input (Larkum et al., 1999. This means that the experimental data does not contain information about calcium dynamics. Therefore this feature could not be optimized here and we only tuned the parameters affecting somatic spiking in the second step. As the recordings were made from young rat (P17-25) pyramidal cells we can use the optimized reduced geometry that we obtained after the first step as it is based on a detailed geometry of a pyramidal neuron from a rat of similar age.
Pinching has a number of effects on neuronal dynamics (compare black and blue traces in Fig. 3): (1) the input resistance increases (from ≈82 M to 131 M ) and hence the spike frequency. By occluding the apical dendrite, one path for current loss is blocked, allowing the somatic membrane to be charged more effectively (Bekkers and Häusser, 2007). (2) The somatic resting potential becomes more hyperpolarized (≈4 mV). It has been shown that the density of I h increases with distance from the soma in the apical dendrite of layer 5 pyramidal neurons (Berger et al., 2001;Kole et al., 2006). When the apical dendrite is blocked during pinching, the depolarizing influence of dendritic HCN channels is reduced, which might explain the hyperpolarization of the somatic resting potential. (3) The threshold for AP initiation becomes lower (≈3 mV) and the AP peak voltage increases (≈0.8 mV), which, at least in a model (Bekkers and Häusser, 2007) could be explained by the decrease in the dendritic capacitance, reducing the electrical  (Stuart and Spruston, 1998) that we used as the starting point to create the reduced model. (b) An illustration of the reduced model (same scale as detailed model, geometrical parameters can be found in Table 1). We divided the complex morphology into four functional sections: The soma, the basal dendrites, the apical dendrites and the tuft. The oblique dendrites are considered to be part of the apical dendrites. (c) We injected a constant current (−1 nA) into the somata of both model neurons and measured the steady-state voltage at different locations. (d) We used a low-amplitude oscillatory somatic input current and measured the resulting membrane potential oscillation to determine the somatic frequency-impedance curve for both models. (e) We calculated the somatic phase-shift between the oscillatory input current and resulting membrane potential oscillation for both models. The black dots and curves in (c-e) describe the passive response properties of the detailed model and served as target functions for the optimization procedure in the first step. The red dots and curves show the corresponding response of the reduced model after optimization.

Fig. 2.
Comparison of the voltage traces in the complex and in the reduced model in response to noisy input current. To test whether the reduced model is a good approximation of the complex model, we analyzed responses to Gaussian noise current injections. The same random current was injected into the somata of the models (green electrodes in a and b and green trace in e). For both models, the somatic voltage as well as the voltage distally (≈280 m and ≈425 m from the soma) was recorded (black and red electrodes in a and b) and the recordings overlaid (traces in c-e).

Fig. 3.
Fitting results for the reduced model 2 after evolving the ion channel conductances using EMOO in the second step of optimization. Four sub-threshold step currents and one supra-threshold step current were delivered, and the neuron's and model's voltage response was recorded in the soma before and during pinching of the apical dendrite. These voltage traces were used for the optimization. The left part of the figure compares the experimental recordings (black) with the model responses (red) for the intact neuron before pinching. The traces in the right part (blue and orange) show the corresponding responses during pinching. A comparison of the black and blue traces (a-c, left and right parts) reveals the effects pinching has on the neuronal response properties. We used the following four objectives for the optimization: (1) the shape of the AP onset (a, left traces); (2) the shape of the AP offset (a, right traces); (3) The interspike interval times (ISIs) of the spike train (b, only the first 600 ms are shown for better visualization). (4) The four sub-threshold traces (c). The five step current injections were Iamp = −0.1 nA, −0.05 nA, 0 nA, 0.05 nA and 0.4 nA (d, green traces). The four objectives (a-c) were determined before and during pinching and compared with the experimental data. The resulting distances were summed up yielding four distance functions that we minimized by EMOO. No calcium channels were present in this step of the optimization. load on the AP initiation site in the axon. (4) The spike afterhyperpolarizing potential (AHP) becomes stronger during pinching (≈3 mV). It is unknown what causes the increase in the AHP. Under normal conditions in the intact cell, back-propagating action potentials (BAPs) (Stuart and Sakmann, 1994) are carried by dendritic voltage-dependent currents, and this in turn leads to dendritic currents flowing back to the soma that counteract the AHP. Upon pinching, this source of depolarization is removed, which could result in an apparent increase of the AHP.
For our target data, we chose a pyramidal layer 5 cell's voltage traces in response to four weak current injections (−0.1 nA, −0.05 nA, 0 nA, 0.05 nA) and one stronger current injection (0.4 nA) that led to spiking responses under both conditions, before as well as during pinching. The goal of parameter optimization using EMOO was to reproduce these data and the effects of pinching. EMOO automatically distributed the active conductances in the axon, soma and dendrites.
A spike was defined as a voltage excursion above a threshold (Â = −20 mV). The spike time was given by the time at which the voltage reaches its maximum, whereas the spike width is measured by voltage crossing the threshold before and after the spike.
In response to supra-threshold stimulation of 0.4 nA, the model had to respond with at least six spikes. Further prerequisites for the model's response at this current were: no spike width could exceed 3 ms; the absolute spike heights (voltage peaks) from the third to the penultimate spike could not change by more than 20%; the voltage minimum between the third and the fourth spike compared to the voltage minimum between the penultimate and the last spike should not change by more than 10%; and, finally, there should not be any interspike interval (ISI) below 15 ms. For all other currents, the model was required not to spike.
If these prerequisites were not met, for both the intact neuron as well as for the model neuron in which the apical dendrite was pinched, the EMOO algorithm severely punished the model by assigning it an extremely high error value. With the prerequisites fulfilled, four squared distance values were measured between the model response and the experimental data: (1) We determined the distances between model and data for each of the four sub-threshold traces between t = −50 ms and 300 ms (relative to the current injection onset, t 0 = 100 ms) and summed these distances up: (2) We determined the squared distance between the average AP onset in the model and the AP onset in the data. For this purpose, the time segment between −0.5 and −0.1 ms before the AP peak was considered, and both the voltage and its time derivative were used: (3) We also determined the distance between the average AP offsets, including their first derivatives (time window ranged from 0.1 ms to 14 ms after the spike peaks): As we included the time derivative for feature extraction our approach resembles the idea of taking the phase-plane of a spike train for the optimization (LeMasson and Maex, 2001). That approach, however, puts more weight on the sub-threshold responses than on the spike shape. This can lead to imperfect fitting results, especially when experimental data is used as a target (Druckmann et al., 2008). We instead focus on the average spike rather than on the whole spike train and hence overcome that limitation. The experimental data we used for the optimization was recorded at 50 kHz ( t = 20 s) and does not provide sufficient time resolution for proper spike alignment and hence for calculating the average spike and its time derivatives. Therefore we applied a cubic spline interpolation (new t = 5 s) to each spike and then aligned and averaged these interpolated spikes (Wheeler and Smith, 1988).
(4) Finally we also determined a distance for each ISI and summed these differences up. This distance function has already proven to be well suited for the optimization of conductance-based models (Keren et al., 2005): For each of the four distances, the value was determined before and during pinching and summed, yielding four final distance values. If all four distances are minimized, then the resulting model reproduces experimental sub-threshold responses, spiking responses as well as detailed AP shape before and during pinching. The effect of pinching was introduced into the model by increasing the axial resistance of the apical dendrite to a high value (10 6 cm).
To minimize the four distance values, we used EMOO with a population size of N = 1000 and a capacity of C = 2000 individuals. 1000 generations were evaluated. The crossover and mutation parameters remained constant at Á c = 10 and Á m = 20 respectively. The probability of a mutation per parameter was 20%. By design, EMOO ends with a population containing models that perform well on one of the distance measures alone, but could be far off the mark in the others. But the final population also contains models that have intermediate fitness in each of the distance functions. It is up to the modeler to choose an individual from the population. In order to do this step automatically as well, we normalized each distance function by the lowest value found in the respective distance function of the last generation. This natural normalization then allowed us to go through all generations and pick the individual for which the sum (=total error, E T ) of these normalized distance values was minimal.
We performed six independent optimization trials, with each trial requiring approximately two days of runtime. The relative improvement of the total error during evolution was 86 ± 5% (measured as (E T (0) − E T (i))/E T (0); E T (0) is the minimal total error in the initial population and E T (i) is the minimal total error in the final Fig. 5. Model prediction of firing frequency. We compare the experimentally measured with the predicted firing frequencies before and during pinching (black and red dots and blue and orange dots respectively) for current injections ranging from 0 nA to 1.1 nA. The current injection' amplitude taken for model optimization was 0.4 nA which is illustrated with the black dashed line. Error bars for the experimental data show the standard deviation of the two measurements for each data point. population) indicating a significant improvement of the fit. We obtained six models that closely reproduce the experimental subthreshold responses, the spike train as well as the AP shape before and during pinching, but have different sets of parameters (models 1-6, Table 1). Fig. 3 shows the fitting results for model 2. The other five models reproduce the data similarly well (not shown).

Optimization without the pinching data
The data set we have used for optimization is unusual, given that the neuron's response with and without the apical dendrite was measured. We wished to quantify how much is gained by using such data. Therefore, we repeated the optimization strategy described above but excluded the distances obtained during pinching and only used model responses and data from the intact neuron. This was done for four independent runs (models 7-10, Table 1). The quality of the fit improved significantly during evolution as shown by a relative improvement of the total error of 82 ± 10%. These models reproduce the experimental recordings for the intact neuron well, including the AP shape. However the prediction of the recordings during pinching is poor in all these models and also variable between models (see Figs. S1 and S2). This shows that the difference in neuronal responses before and during pinching contains useful information about dendritic parameters.

Model evaluation
We checked how well the optimized reduced models generalize to input currents that were not used in the optimization process. For this we compared the responses of model 2 with the experimental data for another current injection (0.6 nA), which the model predicts well (Fig. 4, the other models' predictions were similarly good). Moreover, the model's spike frequency before and during pinching is qualitatively correct for a broad set of current injections (between 0 and 1.1 nA), with a slightly higher firing frequency for stronger input currents (Fig. 5).
We were also interested in how the optimized models predict other experimental findings in pyramidal neurons. Models 1-6 showed a resting potential modulation with distance from the soma of approximately +2 mV ( Fig. 6A and Table 1) which is qualitatively also seen in experiments . Following previous studies (Keren et al., 2009) we used model 2 to evaluate the conductances active at rest along the apical dendrite. It can be seen that the enhanced dendritic depolarization at rest is due to an increased dendritic depolarizing HCN current (Fig. 6B) and due to a lack of hyperpolarizing dendritic potassium currents ( Fig. 6D and E). The modulation of the resting potential in models 7-10 was variable, and in model 10 the neuron's voltage became even more hyperpolarized the farther away from the soma (Table 1). This shows that the pinching data is beneficial to properly constrain dendritic parameters. Next, we tested if the model predicts realistic BAPs (Stuart and Sakmann, 1994;Stuart et al., 1997). In all 10 models somatic APs actively propagated into the apical dendrite while the AP half-width (width at halfway from −60 mV to the AP peak) increased. In models 1-3 and 6-8 the AP amplitude decreased with somatic distance, while in models 4, 5, 9 and 10 the AP amplitude remained nearly constant along the apical dendrite (not shown). Fig. 7 illustrates BAPs for model 2.
Finally we tested how well an average of the optimized models 1-6 would fit the experimental data. The resulting model presented a surprisingly good fit to the data (Fig. S4) and the prediction of the experimental IF -curve appeared even better than in any of the 6 optimized models (Fig. S5). It might be possible that the parameter range leading to good fitting results is broad and that the average model still lies within this range. To test for this we replaced only a single parameter (like R m ) per model with its average, which produced a model that failed reproducing the experimental data. This shows that the quality of the averaged model is rather an indication that certain parameter ratios (that are maintained during averaging) determine the quality of the fit.

Third step: optimizing the calcium spike dynamics
Older pyramidal neurons show elaborate calcium spike dynamics, i.e. a strong dendritic current input induces a local dendritic calcium spike, which can in turn depolarize the soma sufficiently to evoke axosomatic APs Larkum et al., 1999). Moreover, axosomatically initiated spikes can actively back-propagate along the apical dendrite and reduce the current threshold for dendritic calcium spike initiation (Larkum et al., 1999. This threshold reduction translates into a coupling factor between soma and dendrite, estimated to be around 0.5 for pyramidal neurons (Schaefer et al., 2003).
We wanted our reduced model to reproduce these calcium spike dynamics and show how the method can be extended to match rather different data sets. In general, four parameters in the model dominate the calcium dynamics. The strength and initiation threshold of a dendritic calcium spike are determined by the calcium channel density in the tuft (ḡ Cas ) and the voltage shift (vshift Cas ) of this channel. The calcium-dependent potassium channel curtails the length of the calcium plateau and thereby the number of somatic action potentials in a burst. The density of this channel (ḡ KCa ) was also optimized. The apical resistivity (apical Ra), even though it had been previously optimized in the first step, had to be left as a free parameter as well. All remaining parameters that had been previously optimized were not modified.
To illustrate the power of the evolutionary optimization approach, we do not use quantitative experimental data to construct the distance function. Instead, we took general experimental observations of how dendritic calcium dynamics depend on the coupling of the soma and the dendrite (see experimental traces in Fig. 8) and constructed a discrete step-like distance function describing which of these interactions the model qualitatively reproduced. Genetic algorithms for the optimization, as opposed to gradient descent, for instance, can handle such a step-like distance function. For a given parameter combination we set the distance function (or error) E as follows: The HCN channel density increases along the apical dendrite, adding to the depolarization at rest (b). The sodium conductance at rest is negligible and decays linearly with somatic distance (c). Both potassium channel conductances decay rapidly with somatic distance (d, e). At rest, these conductances thus only hyperpolarize the soma.
(1) We determined the threshold somatic current pulse required for a somatic spike. Observing more than a single somatic spike at threshold implies that a single back-propagating AP already elicits a dendritic calcium spike and thereby further somatic spikes. This is not seen in experiments, however. Hence, a model exhibiting multiple spikes was penalized and associated with the highest error value (E = 5000) during optimization.
(2) Now, if at threshold only a single somatic spike was initiated, we next tested whether a strong EPSP shaped dendritic current injection alone could induce a local dendritic calcium spike. We also checked whether this also resulted in a quickly forward spreading Ca spike and eventually multiple somatic spikes. We searched for the dendritic current amplitude threshold (=ICA) that led to this behavior. Based on experimental observations, such a current should elicit a burst of 2-4 spikes in the soma. Yet a somatic spike can occur simply due to depolarization of the soma, without a dendritic spike. So we double-checked whether the first somatic spike was, in fact, due to a somatic depolarization resulting from a dendritic Ca spike. For this purpose, we integrated the voltage in the tuft from −10 ms to 0 ms before the first somatic spike peak. A large value for this integral indicates that the calcium spike was triggered locally; as long as the tuft voltage integral was larger than 500 mV ms we considered the calcium spike to be local.
If no current amplitude was found to produce a locally initiated and forward propagating calcium spike we set the error value to E = 4000.
(3) The next test was to determine how a back-propagating Naspike from the soma influences that threshold. To measure this, we injected a brief somatic current pulse into the soma to initiate a back-propagating AP and searched for the minimal dendritic current threshold (=IBAC) needed to initiate a calcium spike. The resulting calcium spike had to fulfill the following requirements: It should produce a somatic burst of 2-4 spikes and it should consist of prolonged dendritic depolarization with a fast shut off. To quantify these requirements we determined 2 voltage integrals in the tuft (I 1 from 0 ms to 50 ms and I 2 from 100 ms to 150 ms after the first somatic spike respectively): I 1 had to be larger than 500 mV ms while I 2 had to be smaller than 1000 mV ms. If no dendritic current could be found that initiated a calcium spike we set the error value to E = 3000.
(4) If the model passed all the preceding tests, the error was determined completely by the somato-dendritic degree of coupling (C) (Schaefer et al., 2003): Fig. 7. Shape transformation of the back-propagating AP and the underlying ionic conductances in the reduced model. A short current pulse ( t = 5 ms, Iamp = 1 nA) was injected into the soma to elicit an AP. We recorded the voltage of the resulting back-propagating AP at the initial segment (dark red line), in the soma (red line), and in the apical dendrite and tuft (light red lines) (a and b). The initiation of the AP occurs in the initial segment. The amplitude of the back-propagating AP decays with distance from the soma (c), while the half-width (width at halfway from −60 mV to the AP peak) increases with distance (d). The sodium and the fast and slow potassium conductance at the AP peak decay with somatic distance (e-g). These results (black lines) are compared with a passive spread of the somatic AP, for which dendritic Nat channels were removed (dashed black lines). Note the rapid decay in amplitude in (c), implying that that back-propagating APs were not elicited (dashed black line). The illustrated axon in (a) is not to scale.
The target value for C was set to 0.5, and the error was measured as the squared difference between the true coupling and the target.
We used the EMOO-Framework for the optimization with a population size of N = 100 and a capacity of C = 200. The crossover and mutation parameters remained constant at Á c = 10 and Á m = 20, respectively. Parameters mutated with a probability of 20%. As the distance function is step-like, the algorithm accumulates models that fulfill the general experimental observations in the first phase of evolution; once all requirements are met, the algorithm selects those models that come close to the desired somato-dendritic coupling factor. The evolution was stopped once a model was found with a coupling factor of 0.5 or if, after 5 generations, no further error minimization could be observed.
We optimized the calcium dynamics for the 10 models found in the second step. For models 1-3 and 6-8 we obtained parameter combinations satisfying the experimental constraints while for model 4, 5, 9 and 10 we did not. In the latter set, the backpropagating APs do not decay with distance from the soma, which consequently increases the coupling between soma and dendrite such that a single somatic spike always elicits dendritic calcium spikes. Fig. 8 compares experimental recordings with the responses from model 2 (models 1, 3, 6-8 have similarly good waveforms).
We also tested whether the modification of apical axial resistance and the activation of dendritic calcium channels would have an influence on the fitting results from the second step. We only observed minor differences in model behavior (Fig. S3).

Discussion
In this study, we have presented a three-step strategy to automatically generate accurate models of neurons with simplified morphologies, without resorting to tuning parameters by hand. In each step, we defined objective functions for Evolutionary Multi-Objective Optimization, a genetic algorithm (Deb, 2001;Deb et al., 2002;Druckmann et al., 2007).
We fitted the somatic voltage recordings of a layer 5 pyramidal neuron with and without its primary apical dendrite, as this dendrite was reversibly squeezed or pinched to block current flow (Bekkers and Häusser, 2007). These data allowed us to study the (nonlinear) electric load of the dendrite, as seen from the soma.
Some prior knowledge about the distribution of ion channels along the dendrites (Berger et al., 2001) and their kinetics (Hille, 2001) exists and was built into the models. Even with this knowledge and even though the entire geometry of dendritic branches is reduced to a single cable, many free model parameters remained that were subject to optimization (Table 1).
To our surprise, after randomly initializing the algorithm six times, different local densities of ion channels were found each time, but the model in each case fit the voltage traces well (Fig. 3). As in single compartment models (Goldman et al., 2001;Prinz et al., 2003;Schulz et al., 2006) no unique set of optimal model parameters was found; certain parameters can be offset by others (Taylor et al., 2009) which could enable a cell to homeostatically regulate its firing pattern (Prinz et al., 2004). Certain parameters, however, are consistent across optimizations. For instance, the ratio of the sodium channel density in the initial segment and that in the soma was around 40-70, which is in agreement with experimental studies . Moreover, the shift in voltage activation of the transient Na channel, a free parameter in the model, was consistently found to be around −10 mV, which is close to what has been found in experiments (Colbert and Pan, 2002; and has been ascribed to the Na v 1.6 channel (Rush et al., 2005).  (Larkum et al., 1999;Schaefer et al., 2003). The experimental data illustrated in the left part of the figure were kindly provided by Matthew Larkum. (b) A small EPSP shaped dendritic current injection does not lead to spiking. (c) A somatic current pulse elicits a single somatic spike that back-propagates into the dendrite without crossing the threshold for the initiation of dendritic calcium APs. (d) However, if a small dendritic current injection is added the threshold can be crossed and dendritic calcium spikes are elicited and eventually somatic bursting. (e) A calcium spike can also be initiated locally in the dendrite when a strong dendritic current injection is used. The somato-dendritic coupling factor describes the relative reduction of the strong dendritic current injection threshold when a back-propagating AP is present. Black, blue and red electrodes indicate the current injections and voltage recordings in the soma, the apical dendrite and tuft, respectively (a). We used these general experimental observations for the optimization in the third step. The responses of model 2 (right part of the figure) are compared with the respective experimental traces. We used t = 5 ms, t0 = 100 ms for the somatic current pulse and found that a threshold current of Iamp = 0.9 nA was needed to elicit a somatic spike (c and d). The dendritic current injections were modeled as I = Iamp · (1 − e −t/2 ms ) · e −t/8 ms , t0 = 106 ms. Current threshold amplitudes to elicit a dendritic calcium AP were Iamp = 3.6 nA when no somatic spike was present (e) and Iamp = 1.8 nA otherwise (d) resulting in a somato-dendritic coupling of 0.5.
Evolving a population of compartmental models allows one to explore an entire subspace of ion channel properties and densities that give rise to the same response. Golowasch et al. (2002) showed that this parameter space can be non-convex, so that averaging conductance parameters over different models fitted to the data produces a new model that fails. We present here six different models for the same cell as optimized by the evolutionary algorithm, with six very different sets of parameters. Averaging the parameter values, however, across models 1-6 yielded a representative cell model that mimicked the experimental data similarly well as the 6 optimized models on its own (compare Fig. 3 and Fig. S4).
The 'training data' consisted of the somatic voltage response to a fixed current pulse, in the presence or absence of the apical dendrite, whereas the 'test data' were the responses to all other current amplitudes, which were measured experimentally, but not used for optimization. All models generalized well to other amplitudes of injected current (Figs. 4 and 5).
Fitting the model in the second step to both the data with the intact apical dendrite and the data in which the apical dendrite was pinched proved to be essential in narrowing down the degrees of freedom in parameter space. For example, the relative density of HCN channels varied greatly across regions when we only used data from the intact neuron (models 7-10), indicating that a low HCN density in the tuft, for instance, can be compensated by a higher density in the basal dendrites to yield the same somatic voltage response. The density ratio, ranging from 0.15 to 13.03 for optimizations for which only data from the intact cell were taken (models 7-10), narrowed to the range from 1.15 to 2.38 when both data before and after pinching were used (models 1-6). Data from pinching also narrowed down the ratio of the somatic sodium channel density to the sum of somatic potassium channel densities. Without pinching data, this ratio varied by a factor of four across optimizations, but with the pinching data, the ratio varied by a factor of less than two. Moreover, models that were only optimized to predict the voltage response in the case of intact dendritic geometry failed to generalize to the case in which the dendrite was removed.
Of course, observing only the voltage trace at a single point in the neuron means that much of the neuronal dynamics remains hidden. Embedding the voltage data in higher dimensions (Takens, 1981;Eckmann and Ruelle, 1985;Kennel et al., 1992) or examining the relationship between the voltage V and its time derivative dV/dt, known as phase plane analysis (LeMasson and Maex, 2001;Achard and De Schutter, 2006), can uncover some of the hidden dynamics. Several authors have used the phase plane as the basis for fitting neuronal data. Yet many divergent time scales govern the voltage dynamics at the soma, from the slow charging to threshold to the rapid discharge during the action potential. The algorithm used here treats these multiple time scales separately, so that action potential onset, offset, sub-threshold behavior each define an objective function. With multiple objective functions, some parameter changes may improve one objective, while adversely affecting another. Yet the evolutionary algorithm does not enforce hard trade-offs between the objective functions, which is one of the algorithm's core advantages: at each generation, many models compete, and a model that may perform poorly on one objective, but well on another, is not discarded.
Understanding the function of dendrites may well require a more detailed model of dendritic geometry and how ion channels are distributed along the dendrites. Yet even the reduced model presented here is under-constrained, as many configurations of channel densities and properties will give rise to the same observed response in the soma. Some authors have used multiple simultaneous somatic and dendritic recordings to refine computational models (Keren et al., 2005;Metz et al., 2007;Keren et al., 2009); ideally, one would want to observe each compartment's response in isolation.
While the original data used for fitting did not include the dendritic response, we were able to use a new data set from a different experiment to fine-tune the calcium spike in the apical dendrite (Fig. 8). Instead of seeking an exact quantitative fit to disparate experiments, we used a qualitative measure to create a discrete, as opposed to continuous, objective function. The value of this objective function was conditional on the dendritic and back-propagating action potentials occurring in the right temporal sequence, based on the stimulation of the neuron. The success of such an approach demonstrates that evolutionary multi-objective optimization can go beyond finding the least square difference between model and experiment.
Reproducing realistic calcium action potentials and incorporating active properties in the dendrites is crucial to a neuron's ability to discriminate different sequences in time (Larkum et al., 1999;Branco et al., 2010). The computational models' simplified geometry captures both the succession of dendritic events and how the back-propagating action potentials decay along the apical dendrite. In a commonly used model of pyramidal neurons (Mainen and Sejnowski, 1996), back-propagating action potentials do not decay. Indeed, the second step of optimization in our case is indifferent to whether the back-propagating action potentials decay or not. If they do not decay, then the third step of optimization fails. In contrast, if all three steps succeed, the resulting models, as presented here, more accurately reflect the qualitative and quantitative behavior of layer 5 pyramidal neurons. The intrinsic simplicity of such models makes them ideally suited for simulating networks of neurons.
Our method could be further extended. If for example, local Na-APs are required in the basal dendrite (Nevian et al., 2007) this feature could be added to the optimization algorithm in the same fashion as the apical Ca-AP's were tuned. This, however, would require experimental data from the basal dendrite, such as dendritic patch-clamp recordings. It might also be suitable to build a distance function based on knowledge about the decay of the back propagating AP in the basal dendrite.
The method we have presented should directly relate to fitting more complex models using a detailed morphology, and thereby creating models that incorporate features like more fine-grained dendro-dendritic interactions or the influence of synaptic clustering on the firing dynamics (Poirazi et al., 2003), but the use of a detailed morphology would require even larger computing facilities than we had available in this study. One possible strategy might be to start with the optimized reduced model and then to reverse the dendritic simplification process while maintaining the suband supra-threshold properties. Still, a detailed model also requires more computing resources just to simulate it, especially in the case of large networks of such neurons. Which degree of complexity is truly required is a matter of great debate (Herz et al., 2006;Markram, 2006).