Abstract
Postsynaptic integration is a complex function of passive membrane properties and nonlinear activation of voltage-gated channels. Some cortical neurons express many voltage-gated channels, with each displaying heterogeneous dendritic conductance gradients. This complexity has hindered the construction of experimentally based mechanistic models of cortical neurons. Here we show that it is possible to overcome this obstacle. We recorded the membrane potential from the soma and apical dendrite of layer 5 (L5) pyramidal neurons of the rat somatosensory cortex. A combined experimental and numerical parameter peeling procedure was implemented to optimize a detailed ionic mechanism for the generation and propagation of dendritic spikes in neocortical L5 pyramidal neurons. In the optimized model, the density of voltage-gated Ca2+ channels decreased linearly from the soma, and leveled at the distal apical dendrite. The density of the small-conductance Ca2+-activated channel decreased along the apical dendrite, whereas the density of the large-conductance Ca2+-gated K+ channel was uniform throughout the apical dendrite. The model predicted an ionic mechanism for the generation of a dendritic spike, the interaction of this spike with the backpropagating action potential, the mechanism responsible for the ability of the proximal apical dendrite to control the coupling between the axon and the dendrite, and the generation of NMDA spikes in the distal apical tuft. Moreover, in addition to faithfully predicting many experimental results recorded from the apical dendrite of L5 pyramidal neurons, the model validates a new methodology for mechanistic modeling of neurons in the CNS.
Introduction
Neurons in the CNS receive tens of thousands of synaptic inputs on their dendrites and integrate them to create an output in the form of action potentials (APs) that are relayed to the next neuron in the network. This synaptic information is shaped by the active and passive properties of the dendrites. APs initiated at or near the soma actively backpropagate into the dendritic tree (Stuart and Sakmann, 1994). Additionally, dendrites generate complex regenerative Ca2+ and Na+ spikes (Schiller et al., 1997; Magee et al., 1998; Martina et al., 2000; Migliore and Shepherd, 2002; Johnston et al., 2003), and contain electrically and chemically defined compartments (Larkum et al., 1999b, 2001; 2009; Magee, 1999; Korngreen and Sakmann, 2000; Schiller et al., 2000). The exciting finding of regenerative activity in the apical dendrite of neocortical layer 5 (L5) pyramidal neurons (Schiller et al., 1997; Larkum et al., 1999a,b, 2001) has sparked renewed interest in the computational properties of single neurons (London and Häusser, 2005).
Current interpretations of neuronal physiology are based on the watershed investigation by Hodgkin and Huxley (1952) of APs in the giant squid axon. Thus, the golden path to the description of neuronal excitability passes through kinetic modeling of the voltage-gated channels expressed by the neuron, followed by the application of these models to a full numerical model of the neuron. This methodology can be hard to apply, since many CNS neurons express >10 types of ion channels that can display heterogeneous somatodendritic conductance gradients. Given this vast parameter space, compartmental models of single neurons in the CNS are notoriously difficult to tune, which has considerably restricted our ability to understand their integrative properties. Many models simply match specific experimental observations and are thus only crude approximations of neuronal physiology. Recently, numerical simulations have shown that multiple recordings of the membrane potential can be used to accurately constrain compartmental models. This can be achieved by using a simple gradient descent algorithm, if enough recording sites are available in the training dataset (Keren et al., 2005; Huys et al., 2006; Gold et al., 2007). Conversely, constraining a compartmental model solely by using somatic recordings of the membrane potential or reduced features of the dendritic membrane potential leads to the generation of a large solution manifold (Keren et al., 2005; Druckmann et al., 2007, 2008, 2011; Hay et al., 2011). Here we show that a dataset recorded from the soma and apical dendrite combined with a pharmacological and numerical parameter peeling procedure is sufficient to constrain a compartmental model for the apical dendrite of L5 pyramidal neurons containing eight channels. This model faithfully predicts a wide range of physiological activity recorded from L5 pyramidal neurons; it also provides a quantitative mechanism for assessing regenerative dendritic spiking in these neurons.
Materials and Methods
Slice preparation.
Acute brain slices (sagittal, 300 μm thick) were prepared from the somatosensory cortex of 30- to 35-d-old Wistar rats of either sex. The animals were killed by rapid decapitation after shallow anesthesia with isoflurane (Stuart et al., 1993) following the guidelines of the Bar-Ilan University animal welfare committees. Slices were perfused throughout the experiment with oxygenated artificial CSF (ACSF) containing the following (in mm): 125 NaCl, 15 NaHCO3, 2.5 KCl, 1.25 NaH2PO4, 1 MgCl2, 2 CaCl2, and 25 glucose, pH 7.4 with 5% CO2 at 32–34°C. The bath solution for the current-clamp experiments in which all Ca2+-gated K+ channels were blocked contained the following (in mm): 135 NaCl, 15 HEPES, 2.5 KCl, 1 MgCl2, 2 CaCl2, 25 glucose, 1 mm tetraethylammonium (TEA; Sigma), pH 7.4 with NaOH. Apamin (Alomone Labs) and iberiotoxin (Alomone Labs) were stored at −20°C as stock solutions in double-distilled water and added directly to the bath solution to form final concentrations of 200 and 30 nm, respectively. Before adding iberiotoxin to the bath solution, the perfusion tubing was coated with Sigmacote (Sigma) to prevent binding of the toxin. In addition, in experiments with iberiotoxin, 0.1 mg/ml bovine serum albumin (Sigma) was added to the application solution to prevent nonspecific binding by the iberiotoxin.
Electrophysiology.
Pyramidal neurons from L5 were visually identified using infrared differential interference contrast videomicroscopy (Stuart et al., 1993). Patch pipettes (4–7 MΩ somatic, 7–12 MΩ dendritic) were pulled from thick-walled borosilicate glass capillaries (2.0 mm outer diameter, 0.5 mm wall thickness; Hilgenberg). The standard pipette solution contained the following (in mm): 125 K-gluconate, 20 KCl, 10 HEPES, 4 MgATP, 10 Na-phosphocreatin, 0.5 EGTA, 0.3 GTP, and 0.2% biocytin, pH 7.2 with KOH. A liquid junction potential of 10 mV was not corrected for. Whole-cell recordings from the soma and apical dendrite were performed at 36–38°C using a Multiclamp-700B amplifier (Molecular Devices). Voltage was filtered at 10 kHz and sampled at 50 or 20 kHz using AxoClamp9 (Molecular Devices), digitized by a Digidata-1320 interface (Molecular Devices), and stored on the hard disk of a personal computer. At the end of each experiment, slices were fixed in cold 100 mm PBS, pH 7.4, containing 4% paraformaldehyde. After fixation, the slices were incubated for 2 h in avidin-biotinylated horseradish peroxidase (ABC-Elite, Vector Laboratories), and the stain was developed using 0.015% diaminobenzidine. The stained neurons were digitally traced using a Neurolucida system (MicroBrightField), and the tracings were converted to NEURON readable code. Care was taken that, in NEURON, no compartment exceeded 50 μm in length.
Numerical simulations.
The compartmental model, genetic algorithm, and cost function were programmed using NEURON version 5.9 (Carnevale and Hines, 2006). Simulations were run on a custom-made Linux cluster with 168 central processing units sharing the same network file system. One of the machines functioned as a master, submitting and managing the jobs using the ParallelContext class of NEURON over a MPICH2 ring spanning the other computers. Ion channel models were implemented using the NMODL extension of NEURON (Carnevale and Hines, 2006). All simulations were performed with an integration time step of 50 μs.
Compartmental model.
We slightly modified our previous model and applied the assumptions used in the study by Keren et al. (2009), where full details are given. The assumptions were as follows: the morphology of the neuron is known; the passive membrane parameters are spatially uniform; and intracellular Ca2+ dynamics were simulated as a Ca2+ extrusion by a CaATPase (Destexhe et al., 1993), in which the decay time constant could be viewed as a simple intracellular Ca2+ buffer:
where [Ca2+]i is in millimoles, the unit conversion constant is k = 10,000 for ICa (in microamperes per square centimeter), F = 96489 C/mol is the Faraday constant, d = 0.1 μm is the depth of the shell beneath the membrane, [Ca2+]rest = 10−5 mm is the intracellular Ca2+ concentration at rest, [Ca2+] is the change in the Ca2+ concentration during an event, and τ = 80 ms is the rate of Ca2+ removal that was adapted to cortical neurons (Schaefer et al., 2003); the model contained eight voltage- and Ca2+-gated ion channels. These included an Na+ channel, a slow inactivating K+ channel, a fast inactivating K+ channel, and a hyperpolarization-activated cation current (Ih) channel (Keren et al., 2005, 2009). Two voltage-gated Ca2+ channels, a high voltage-activated (HVA) Ca2+ channel and a medium voltage-activated (MVA) Ca2+ channel were based on nucleated patches (Almog and Korngreen, 2009) and whole-cell recordings (Foehring et al., 2000; Magistretti et al., 2000). The models of the voltage-gated Ca2+ channels were based on the Goldman–Hodgkin–Katz equation, which appeared to be more sensitive to the intracellular Ca2+ changes within the neuron. An MVA Ca2+ model with a slow inactivation time constant (in milliseconds) was chosen. Finally, two Ca2+-gated K+ channels, a small-conductance Ca2+-gated K+ (KSK) channel, and a large-conductance Ca2+-gated K+ channel (KBK) were based on whole-cell recordings (Khaliq et al., 2003; Sun et al., 2003; Akemann and Knöpfel, 2006; Mercer et al., 2007; Deister et al., 2009). All the channels were modeled using models based on the study by Hodgkin and Huxley (1952). The units used in reporting the parameters are as follows: the permeability of the Ca2+ channels is given in micrometers per second; the conductance of the ion channels is given in picosiemens per square micrometer; the time constant is given in milliseconds; the membrane potential is given in millivolts; the permeability of the HVA Ca2+ channel (CaHVA) is given as:
permeability of the MVA Ca2+ channel (CaMVA) is given as:
conductance of the KSK channel is given as:
conductance of the KBK channel is given as:
and the kinetics of the channels are known and were based on experimental recordings.
Dendritic channel gradients were a continuous function of the distance from the soma. Reconstructed axons were removed. Where stated, the axon was replaced by an artificial axon (Mainen and Sejnowski, 1996; Schaefer et al., 2003) to allow AP initiation. In the myelinated axon, the membrane capacitance (Cm) was reduced to 0.04 μF/cm2, and the passive membrane conductance (Gpassive) was changed to 0.02 pS/μm2 in the nodes of Ranvier. Channel densities were as follows in the axon hillock initial segment in the node of Ranvier: GNa = 30,000, GKs = 1500, and GKf = 1000 (where GNa, GKs, and GKf are the maximal conductances per unit area of the sodium, slow inactivating potassium, and fast inactivating potassium conductances, respectively). Channel densities were as follows in the myelinated axon: GNa, GKs, and GKf as the soma (Table 1, cell 5). The full details of the basic model are given in (Keren et al., 2005, 2009). The reversal potentials of the voltage- and Ca2+-gated ion channels were set at 60, −80, 130, and −30 mV, respectively, for Na+, K+, Ca2+, and Ih. The NEURON code of the model and the genetic algorithm were donated to the ModelDB database.
Results of fitting the model to the data using a linear gradient for the MVA Ca2+ channel
Genetic algorithm.
A genetic algorithm is an optimization algorithm based on the mechanisms of Darwinian evolution. It uses random mutation, crossover, and selection operators to breed better models or solutions (individuals) from an originally random starting population (Mitchell, 1996). The current study used a genetic algorithm similar to those used previously (Keren et al., 2005, 2009; Gurkiewicz and Korngreen, 2007). Briefly, the population was sorted according to the value of the cost function of each individual (Eq. 1), and a new generation was created using selection, crossover, and mutation as operators. Selection used a tournament in which two pairs of individuals were randomly selected and the individual with the better score from each pair was transferred to the next generation. This procedure was repeated N/2 times (where N is the size of the population) until the new population was full. The one exception to this selection process (and later to the crossover and mutation operators) was that the best individual was transferred unchanged to the next generation to prevent genetic drift.
Each pair selected for transfer to the new population was subjected to a one-point crossover operator with a probability of 0.5. After the new population was created, each parameter value in the new population was subjected to mutation with a probability of 0.1. The mutation operator performed a substitution of the parameter value with a random value drawn from a flat random number distribution spanning the entire search space of the parameter. On average, a typical run of the genetic algorithm lasted 2–3 d on our Linux cluster.
The cost function calculated the difference between the target and the test membrane potential traces, and summed its mean squares, yielding a cost function value expressing the distance of the test from the target dataset, as follows:
where V represents the target dataset, and v represents the membrane potential changes in the test dataset. N is the total number of points in each membrane potential trace, and M is the number of membrane potential sweeps simulated in the model or recorded in the experiment.
Data analysis.
Data were analyzed using IgorPro and Excel software. The influence of the blockers on the membrane potential was examined by fitting a linear function to the repolarization phase of a single AP. This fit gave the repolarizing velocity for the three different blockers and each concentration used. For each blocker, the repolarizing velocity was normalized to the control value. The normalized repolarizing velocity was plotted versus the log concentration of each blocker. The IC50 values were calculated using the following equation:
where Δ is the relative repolarizing velocity change, base is the normalized repolarizing velocity value of the bottom plateau of the curve itself, blocker is the blocker concentration, and IC50 is for the blocker. The voltage time integral was calculated for simulated dendritic EPSPs and spikes to compare the response to simultaneous activation of NMDA and AMPA synapses (Larkum et al., 2009; ModelDB accession no. 124043) to the expected response (summation of individual activation of the synapses).
Statistical analysis.
The statistical analysis used a one-tailed Student's t test (*, **p < 0.00005). All values are displayed as the mean ± SE, unless otherwise stated. The coefficient of variation (CV) was calculated for each parameter in the model. The confidence limits for the mean gradients were calculated with a significance level of 0.1 (α).
Results
Blocking all voltage-gated Ca2+ channels by substituting the Ca2+ bath with Co2+ (Keren et al., 2009) previously enabled us to reduce the parameter space in a model of L5 neocortical pyramidal neurons. Using somatodendritic recordings obtained under these conditions allowed a genetic optimization algorithm to constrain a simple model for the apical dendrite of L5 pyramidal neurons in several optimization steps (Keren et al., 2009). Here we extended this pharmacological parameter peeling procedure to constrain a full model for the apical dendrite of L5 pyramidal neurons. We designed the current-clamp recordings required to obtain a new dataset isolating the effect of the voltage-gated Ca2+ channels during a Ca2+ spike by blocking all the subtypes of Ca2+-gated K+ channels. First, pharmacological experiments determined the required concentration of blockers of the Ca2+-gated K+ channels. Figure 1a shows the effect of 100 nm apamin, a blocker of KSK channels, on the slow after-hyperpolarization of an AP recorded at the soma (Fig. 1a, red trace). Figure 1b presents the apamin concentration dependence of the normalized repolarizing velocity of a somatic AP [i.e., the velocity (y-axis) as a function of apamin concentration (x-axis)] with an IC50 of 31 ± 6 nm (n ≥ 4).
Pharmacological block of Ca2+-gated K+ channels. a, A somatic action potential measured from a L5 pyramidal neuron before (black trace) and after (red trace) adding apamin (100 nm) to the bath solution. b, The relationship between apamin concentration and inhibition of the repolarizing velocity of the action potential plotted on a logarithmic concentration scale (n ≥ 4). The repolarizing velocity was normalized to the minimal value. The data were fitted with the Hill equation (solid line), giving an IC50 of 30.82 ± 5.89 nm. Error bars indicate normalized SEM. c, A somatic action potential measured from a L5 pyramidal neuron before (black trace) and after (red trace) adding iberiotoxin (30 nm) to the bath solution. d, A somatic action potential measured from an L5 pyramidal neuron before (black trace) and after (red trace) adding TEA (1 mm) to the bath solution. e, Comparison of the repolarizing velocity of the action potential in the different bath solutions: control solution (ACSF, n = 9), solution containing TEA (1 nm, n = 4); and solution containing iberiotoxin (30 nm, n = 4). The asterisk indicates a significant difference (p < 0.00005, one-tailed t test) between the control solution and the solution containing TEA. The double asterisk indicates a significant difference (p < 0.00005, one-tailed t test) between the control solution and the solution containing iberiotoxin. Error bars indicate SEM. f, The relationship between TEA concentration and inhibition of the repolarizing velocity of the action potential plotted on a logarithmic concentration scale (n = 4) as in b. The fit gave an IC50 of 0.58 ± 0.046 mm. Errors bars indicate normalized SEM.
Next, we tested the effect of blockers for the KBK channel. Figure 1c–f displays the effect of iberiotoxin, and the less specific blocker of voltage-gated K+ channels TEA, which in small concentrations mainly blocks BK channels (Storm, 1987; Lang and Ritchie, 1990). Iberiotoxin (Fig. 1c, 30 nm, red trace) and TEA (Fig. 1d, 1 mm, red trace) had similar effects on the shape of the AP recorded at the soma; namely, the broadening and removal of the fast afterhyperpolarization. Figure 1e compares the blockade effect of the two blockers on the repolarizing velocity of a somatic AP. We observed no significant difference in the blockade level of these two blockers. Figure 1f shows the concentration dependence of the blocker TEA as a function of the normalized repolarizing velocity of a somatic AP, with an IC50 of 0.58 ± 0.04 mm (n = 4). Thus, we decided to block the BK channels with a low concentration of TEA instead of using the toxin iberiotoxin.
Based on these pharmacological experiments, the SK channel was blocked with apamin (200 nm), and the BK channel with TEA (1 mm). Simultaneous whole-cell recordings from the soma and along the apical dendrite of L5 pyramidal neurons recorded the membrane potential before and after bath application of this blocker cocktail. Following each recording, the recorded cell was stained, the tissue fixed, and the neuron traced using Neurolucida. The morphology was converted to a code readable by NEURON. This experimental protocol provided an electrophysiological and a morphological dataset for each neuron.
We used this dataset to extend the previous model for L5 pyramidal neurons (Keren et al., 2009). Based on previous work (Foehring et al., 2000; Magistretti et al., 2000; Almog and Korngreen, 2009), we added two types of Ca2+ channel models to this basic model: an HVA Ca2+ channel and an MVA Ca2+ channel. We previously pharmacologically identified five possible voltage-gated Ca2+ channel types in somatic nucleated patches (Almog and Korngreen, 2009). Thus, the reduction to two channel types in the current model is a simplification that may have to be addressed in future work once more accurate dendritic gradients of these channels are available. Additionally, we inserted two types of Ca2+-gated K+ channels models into the model; namely, SK and BK channels (Khaliq et al., 2003; Sun et al., 2003; Akemann and Knöpfel, 2006; Mercer et al., 2007; Deister et al., 2009).
Figure 2 illustrates the application of the numerical peeling procedure to one of the five experimental datasets used here. We recorded hyperpolarizing and depolarizing membrane potential traces by negative and positive somatic current injections before (Fig. 2b, right) and after (Fig. 2b, left) blocking Ca2+-gated K+ channels. Next, a genetic algorithm constrained the model in two optimization steps. The average densities of our previously constrained four-conductance model (Keren et al., 2009) were used as an initial estimate for these conductances. We allowed the range of each of the four basic conductances to vary within the 95% confidence range of our basic model (Keren et al., 2009). The first step optimized a model without Ca2+-activated channels using the dataset recorded in the presence of blockers to these channels (Fig. 2b, left). This first step optimized the density gradients of the MVA and HVA Ca2+ channels (Fig. 2b, right). To prevent these optimized parameters from drifting away, the search range for these parameters was then set to ±25% of the obtained values during the second optimization step in which the full model was optimized using the dataset recorded in standard ACSF (“control conditions”). A linear function described the dendritic conductance gradient for the four new conductances and permeabilities [HVA Ca2+ (PHVA), MVA Ca2+ (PMVA), SK (GSK), and BK (GBK)]. This led to a set of parameters describing the passive membrane parameters, the GNa, GKf, GKs, GIh, GSK, and GBK dendritic conductance density gradients, and the PHVA and PMVA dendritic permeability density gradients (Table 1).
Application of the parameter peeling procedure to recordings from the apical dendrite of L5 pyramidal neurons. a, Reconstruction of a L5 pyramidal neuron stained with biocytin illustrating electrode placement. b, Constraining passive membrane parameters and channel density gradients of voltage-gated K+, Na+, Ca2+, and Ih channels (left), and Ca2+-gated K+ channels (right). Hyperpolarizing and depolarizing membrane potential traces were recorded at the soma (black traces) and at 415 μm along the apical dendrite (red traces). Ca2+-gated K+ channels were blocked by apamin (200 nm) and TEA (1 mm). The dashed blue traces show the membrane potential traces simulated at the dendrite using the best parameter set obtained by the genetic algorithm. Table 1 (cell 1) shows the final parameter set obtained from this optimization. c, In a simulation, the somatic membrane potential of this cell was clamped to an experimental waveform series of four APs at 41 and 157 Hz that were recorded from cell 1 (black traces; Table 1). The membrane potential traces recorded at 415 μm along the dendrite are shown in red. The dashed blue traces represent the membrane potential traces simulated at the dendrite using the best parameter set obtained by the genetic algorithm.
We used the somatic current-clamp recordings in the simulations as voltage-clamp commands and trained the genetic algorithm using only the dendritic membrane potential. This procedure provided a good fit for the membrane potential traces recorded with apamin and TEA. It also allowed adequate fitting of the traces recorded without any ion channel blockers (Fig. 2b). To test the reliability of the fit, we clamped the somatic membrane potential of the optimized model to a series of four APs at low and high frequencies. This led to a good reproduction of the measured dendritic data (Fig. 2c, dashed blue traces).
Figure 3 shows the parameter sets acquired by applying the peeling procedure to recordings of five cells (see also Table 1). The majority of the passive membrane parameters took on similar values in all five cells (Fig. 3a; Table 1), whereas the membrane resistance (Rm) showed a large variability. The gradients of the dendritic distributions of GKs GKf, GIh, and GNa (Fig. 3b; Table 1) were similar to previous experimental reports (Stuart and Sakmann, 1994; Berger et al., 2001; Schaefer et al., 2007; Keren et al., 2009). The mean linear dendritic gradient density of PHVA and PMVA decreased along the apical dendrite (Fig. 3c,d) and showed the largest variability in the five cells (Table 1). The mean dendritic gradient density of the GSK also decreased along the apical dendrite, whereas the mean gradient density of GBK was uniformly distributed over the dendrite.
Conductance and permeability gradients in the optimized model. a, Passive parameters obtained by the peeling procedure in five experiments. b, The average dendritic conductance density of GNa (pink line), GKf (blue line), GKs (red line), and GIh (gray line). c, The average dendritic permeability density of PMVA (solid purple line). The dotted purple lines are the 90% confidence limits of the average. Inset, The average dendritic permeability density of PHVA (solid purple line). The dotted orange lines are the 90% confidence limits of the average. The dendritic permeability gradient of CaMVA was calculated using the following equation: PMVA(x) = PMVA,soma + PMVA,dend(exp ((−(x − CaMVA,dist)/CaMVA,width)2)) d, The average dendritic permeability density of PHVA (solid orange line). The dotted orange lines are the 90% confidence limits of the average. e, The average dendritic conductance density of GSK (solid dark blue line). The dotted dark blue lines are the 90% confidence limits of the average. f, The average dendritic conductance density of GBK (solid green line). The dotted green lines are the 90% confidence limits of the average.
Parameters can show substantial variability after optimizing a model by a genetic algorithm (Keren et al., 2005; Druckmann et al., 2007, 2008; Gurkiewicz and Korngreen, 2007; Hay et al., 2011). Optimizing the model in several stages using the peeling procedure reduces this variability (Keren et al., 2009). We have shown that a CV of the parameter, calculated between models optimized using different datasets, can be used as an indication of the robustness of that parameter (Keren et al., 2009). Here we calculated the CV of five optimized parameter sets (Table 1). The optimized model predicted that the dendritic density of the BK channel would be uniform (Fig. 3f). As a result, the dendritic location at which the linear gradient levels off (Table 1, BKdist) was poorly constrained by the genetic algorithm, leading to a high CV for this parameter. We only included two voltage-gated channels in the model instead of the five that have been pharmacologically identified (Almog and Korngreen, 2009). To partially adjust for this reduction, we allowed the genetic algorithm to shift the activation and inactivation curves of these two pooled Ca2+ channels. The parameters describing the shift to the activation curves of these two channels displayed low CVs, indicating that the model was sensitive to changes in these parameters (Table 1, CaMVA,shift,act and CaHVA,shift,act). However, the parameters controlling the shift to the inactivation curves (Table 1, CaMVA, shift,inact and CaHVA,shift,inact) displayed very high CVs. This clearly indicated that the optimization process was insensitive to changes in the inactivation curves of both Ca2+ channels (HVA and MVA) within the search ranges that were defined for the genetic algorithm. Interestingly, the somatic value of the HVA permeability displayed a large CV, possibly indicating that the model of the somatic Ca2+ channels should be expanded (Table 1).
Previous studies have hypothesized that Ca2+ channels may form a high-density or “hot” Ca2+ zone at the distal dendrite, thus enabling the generation of dendritic Ca2+ spikes (Yuste et al., 1994; Larkum and Zhu, 2002; Schaefer et al., 2003; Larkum et al., 2009; Hay et al., 2011). Nevertheless, our optimized model successfully simulated a Ca2+ spike at the distal dendrite using a linear dendritic gradient density of two Ca2+ channels. To test whether a hot Ca2+ zone at the distal dendrite could also generate a dendritic Ca2+ spike, we optimized the model when a Gaussian-shaped hot zone was added to the spatial profile of the MVA Ca2+ permeability. Initially, the genetic algorithm freely modified the height, width, and location of the hot zone. This resulted in a hot zone that was not located in one specific place along the apical dendrite. When we constrained the location of the hot zone to the region between 400 and 700 μm along the apical dendrite, we obtained a good fit displaying a hot zone with a twofold increase in MVA Ca2+ permeability (Fig. 3c, insert). Although the genetic algorithm was able to fit the model to the measured dendritic data, it failed to generate a Ca2+ spike at distal locations on the apical dendrite (data not shown).
Next, we tested the ability of the model to predict the physiology of the apical dendrite of L5 pyramidal neurons. We aimed to determine whether the model could predict an ionic mechanism for this activity. Injecting current distally into the apical dendrite of L5 pyramidal neurons can generate a local or “isolated” dendritic spike (Schiller et al., 1997). Thus, we simulated a dendritic spike by injecting 50 ms current steps into the simulated apical dendrite (using linear gradients of MVA and HVA Ca2+ permeabilities). Consistent with experimental observations, the optimized model generated a local distal dendritic spike and a somatic subthreshold membrane potential deflection (Fig. 4a) predicting a mechanism for the dendritic spike (Fig. 4b). This mechanism had the following three phases: (1) The GNa and GKf conductances were activated and inactivated during stimulus onset (Fig. 4b1), and this initial channel activation boosted the membrane potential, thus bringing it closer to the activation threshold of the voltage-gated Ca2+ channels; (2) the two voltage-gated Ca2+ channels were activated by this initial trigger and continued to depolarize the membrane potential (Fig. 4b2); and (3) the GKs, GBK, and GSK conductances were activated, repolarizing the membrane potential (Fig. 4b3). The Ih conductance decreased during the Ca2+ spike (Fig. 4b4) and recovered slowly once the membrane potential returned to its resting value. In Figure 4b1, the GKf conductance is activated before the GNa conductance as a function of the different activation thresholds of these two channels (the midpoint activation voltage for Kf conductance was −47 mV, and for the Na conductance was −36 mV; Keren et al., 2005).
Initiation of a local dendritic regenerative spike. a, Reconstruction of an L5 pyramidal neuron stained with biocytin illustrating electrode placement (left). Right, Simulated somatic and dendritic voltage responses to 50 ms current injection (0.6–1 nA) through the dendritic pipette at 600 μm, as illustrated on the left. A local dendritic spike was evoked by a 1 nA current. The optimized model used for this simulation was cell 5 (Table 1). b, The activation of the eight ion channels inserted into the model during the dendritic regenerative potential (black trace) simulated in a. c, The simulated threshold for dendritic potential (●) and somatic AP (○) is plotted as a function of the distance from soma. The solid lines are an exponential fit. To obtain a simulated somatic AP, an artificial axon was added to the model; the dendritic Na+ conductance density was divided by 10 to simulate the passive forward propagation toward the soma.
The shallow gradients of the MVA and HVA channels suggested that it should be possible to generate a local dendritic Ca2+ spike both distal and proximal to the soma. However, experimental recordings have only reported spikes distal to the soma (Schiller et al., 1997; Zhu, 2000; Larkum et al., 2001, 2009). To investigate this difference, we generated local spikes at several locations along the apical dendrite of the optimized model (Fig. 4c). As mentioned above, the optimization was performed only for the apical dendrite, thus eliminating the need to define an AP initiation site at the axon initial segment by applying a somatically recorded membrane potential as a voltage-clamp command during the optimization process. Thus, it was possible to inject current pulses into the apical dendrite without triggering an axonal AP. As predicted by the dendritic gradients, it was possible to generate a spike both proximal and distal to the soma (Fig. 4c, ●). However, the current threshold for generating a local spike increased as the current was injected more proximally to the soma. This result can be easily explained by considering the electrotonic structure of the apical dendrite. Proximal to the soma, the dendrite is thicker than its distal part, and the proximal part also sprouts several oblique dendrites, whereas the distal part is almost bare. Thus, much of the proximally injected current flows to the current sinks generated by oblique dendrites and the soma (Schaefer et al., 2003). Furthermore, the membrane area of the proximal compartment is larger than that of the distal compartment, requiring more current to charge it to the same potential.
To qualitatively investigate the relationship between the threshold of the Ca2+ spike and that of the axonal AP, we attached a model for a cortical axon to the soma (Mainen and Sejnowski, 1996; see also Materials and Methods). This axon was not optimized by the genetic algorithm since we did not have experimental recordings from the axon's initial segment. Injecting a square pulse at various locations along the apical dendrite generated APs at the axon initial segment. As expected, the current threshold of the axonal APs increased as the injection site in the apical dendrite was moved more distally along the dendrite (Fig. 4c). While qualitative, the intersection between the curves displayed in Figure 4c provided a simple explanation for the appearance of dendritic spikes solely in the distal apical dendrite. Current injection into the proximal apical dendrite will passively propagate to the soma and will trigger an AP at the axon initial segment that fails to initiate a dendritic spike due to the electrotonic structure of the proximal dendrite. Conversely, current injection to the distal apical dendrite will attenuate to the soma, failing to generate an AP, but will be large enough to charge the local dendritic membrane to the threshold of the dendritic spike.
It has been shown that a train of low-frequency AP only backpropagates into the dendrite, whereas a high-frequency train generates a dendritic Ca2+ spike (Larkum et al., 1999b). We tested the ability of our model to predict this physiological observation. We voltage-clamped the soma of the optimized model using an experimentally recorded series of four APs (Fig. 5), and simulated the membrane potential and channel activation along the apical dendrite. At a low AP frequency, the APs backpropagated into the apical dendrite (Fig. 5), in accordance with experimental findings (Stuart and Sakmann, 1994; Larkum et al., 1999b). During backpropagation the Na+ channel and the two Ca2+ channels (MVA and HVA) were responsible for dendritic membrane depolarization up to ∼400 μm from the soma (Fig. 5, orange and green). At more distal locations, backpropagation was almost completely passive. The main channels repolarizing the dendritic membrane potential were the two Ca2+-gated K+ channels (SK and BK), whereas the somatic membrane was also repolarized by the two voltage-gated K+ channels (GKf and GKs; Fig. 5, black). Figure 6 shows a simulation similar to that in Figure 5 but with a high-frequency somatic AP train. At a high frequency, the optimized model succeeded in generating a dendritic Ca2+ spike beyond 600 μm (Fig. 6, blue and red). At this frequency, the Na+ channel contributed mainly to the first spike along the dendrite (Fig. 6). The main channels repolarizing the dendritic membrane were the two Ca2+-gated K+ channels (GSK and GBK; Fig. 6). The somatic membrane was also repolarized by the two voltage-gated K+ channels (GKf and GKs; Figs. 5, 6, black). The MVA Ca2+ channel contributed more than the HVA Ca2+ channel during the dendritic Ca2+ spike (Fig. 6), as pointed out in Figure 4b. The optimized model faithfully predicted experimental patterns of intracellular Ca2+ concentration (Larkum et al., 1999b). Following a low-frequency train of backpropagating APs, the intracellular Ca2+ concentration increased at the soma and the proximal apical dendrite (Fig. 7, left). However, following a high-frequency train of APs the intracellular Ca2+ concentrations increased from the soma all the way to the distal apical dendrite (Fig. 7, right).
Conductance and permeability activation during a low-frequency series of backpropagating action potentials. In a simulation, the somatic membrane potential was clamped to a waveform derived from a series of four APs with a firing rate of 39 Hz (black lines) recorded from cell 5 (Table 1). The changes in the membrane potential, the conductance of the six voltage-gated ion channels (Na+, Kf, Ks, Ih, SK, and BK), and the permeability of the two voltage-gated Ca2+ channels (CaHVA and CaMVA) included in the model are shown at the soma and at 200, 400, 600, and 800 μm along the apical dendrite. The traces are color coded according to the color of the simulated pipettes shown on the left.
Conductance and permeability activation during a high-frequency series of backpropagating action potentials. In a simulation the somatic membrane potential was clamped to a waveform consisting of a series of four APs with a firing rate of 147 Hz (black traces) recorded from cell 5 (Table 1). The changes in the membrane potential, the conductance of the six voltage-gated ion channels (Na, Kf, Ks, Ih, SK, and BK), and the permeability of the two voltage-gated Ca2+ channels (CaHVA and CaMVA) included in the model are shown at the soma and at 200, 400, 600, and 800 μm along the apical dendrite. The traces are color coded according to the simulated pipettes shown on the left.
Intracellular Ca2+ concentration during low- and high-frequency stimulation. In a simulation, the somatic membrane potential was clamped to a waveform comprising a series of four APs with a firing rate of 39 Hz (left) and 147 Hz (right) (black traces) recorded from cell 5 (Table 1). The parameter Cai in the optimized model was used to track the intracellular Ca2+ concentration changes that are shown at the soma and at 200, 400, 600, and 800 μm along the apical dendrite. The units of the intracellular Ca2+ concentration are millimolar due to the absence of a Ca2+ buffer in the model. The traces are color coded according to the simulated pipettes shown on the left.
It has been shown that temporal coupling of a backpropagated AP with synaptic inputs generates a backpropagation-activated Ca2+ spike (BAC) firing (Larkum et al., 1999a,b, 2001). To test whether our model can predict BAC firing, we again added an artificial axon to the cell body (Mainen and Sejnowski, 1996). Then we injected a subthreshold EPSP-like current of 1.6 nA into the distal apical dendrite (at 800 μm; Fig. 8b, gray). This potential attenuated along the dendrite (Fig. 8b, red and blue) toward the soma (Fig. 8b, black). A somatic current injection of 0.5 nA generated a somatic AP that backpropagated into the apical dendrite (Fig. 8c). Coupling the subthreshold EPSP-like current (Fig. 8b, gray) and the somatic current (Fig. 8b, black) with a time window of 5 ms caused the optimized model to generate a BAC firing (Fig. 8d, red), followed by two additional somatic APs (Fig. 8d, black). Injecting a higher EPSP-like current (3 nA) into the distal apical dendrite also generated a Ca2+ spike and somatic APs (Fig. 8e). BAC firing is highly dependent on the timing of the backpropagating AP and dendritic synaptic input (Larkum et al., 1999a, their Fig. 2). Our model was able to predict this temporal dependency. Figure 9 displays the temporal coupling of a backpropagated AP with the synaptic inputs required to simulate BAC firing. The time interval for a BAC firing was up to 7 ms (Fig. 9b). When the delay between the AP and the EPSP was changed to −10 or 10 ms, the optimized model did not generate BAC firing (Fig. 9a,c).
Reproducing BAC firing. a, Reconstruction of an L5 pyramidal neuron stained with biocytin, illustrating electrode placement. b, EPSP-like current of 1.6 nA (rising τ = 2 ms, declining τ = 10 ms) was injected through the gray pipette (800 μm, bottom) into a dendrite of one of the optimized models (Table 1, cell 5). The simulated voltage response showed the shape of an EPSP at the soma (black trace) and at the apical dendrite (400 μm, blue trace; 600 μm, red trace). c, Injection of a square current of 0.5 nA through the somatic pipette (black, bottom). The action potential generated at the soma backpropagated along the apical dendrite. d, The combination of the two stimuli used in b and c with a time interval of 5 ms generated a BAC firing at the distal apical dendrite (red trace). e, EPSP-like current injection of 3 nA (rising τ = 2 ms, declining τ = 10 ms) through the gray pipette (800 μm, bottom) was sufficient to generated a Ca2+ spike at the distal apical dendrite (red trace).
The timing of BAC firing. Top, EPSP-like current of 0.6 nA (rising τ = 1 ms, declining τ = 5 ms) was injected at the distal dendrite (800 μm) into a dendrite of one of the optimized models (Table 1, cell 5). The simulated voltage response showed the shape of an EPSP at the soma (black trace) and at the apical dendrite (600 μm, red trace). Middle, Injection of a square current of 0.5 nA through the somatic pipette (black, bottom). The action potential generated at the soma backpropagated along the apical dendrite. Bottom, Combining the two stimuli used with a different time interval (a, −10 ms; b, 7 ms; c, 10 ms). A Ca2+ spike was simulated using a time interval of 7 ms. Δt was calculated as the time between the start of the somatic injection and the start of the dendritic injection.
In our optimized model, the distal EPSP and the backpropagating AP collided at ∼470 μm from the soma. Thus, we examined the mechanism at 470 μm along the apical dendrite where the two stimuli converged, and at 600 μm from the soma, where the dendritic spike was fully developed (Fig. 8d, red). At 600 μm during the BAC firing (Fig. 10b), the multistep mechanism was similar to that observed for a dendritic spike, as shown in Figures 4 and 6. First, the GNa and GKf conductances were activated and were basically responsible for the first Na+ spike (Fig. 10c, left). Unlike Figure 4, the GNa conductance was activated before the GKf conductance, since a Na+ spike first propagated from the soma and was not generated locally. Next, the two Ca2+ channels, but mainly the PMVA, were activated and were responsible for the continued depolarization (Fig. 10b, middle). The three K+ conductances (GKs, GBK, and GSK), but mainly the GKs conductance, repolarized the dendritic membrane (Fig. 10b, right). The mechanism at 470 μm was similar, except that the PMVA permeability had less impact on the membrane depolarization (Fig. 10c, middle) and the GKs conductance also had less influence on membrane repolarization (Fig. 10c, right). These differences fit the absence of a Ca2+ spike at 470 μm. Here there were only backpropagating APs (Fig. 10c, black).
Conductance and permeability activation during BAC firing. a, Reconstruction of an L5 pyramidal neuron stained with biocytin illustrating electrode placement. b, The activation of seven ion channels inserted into the model during the BAC firing simulated in Figure 5d at 600 μm from the soma. Left, Activation of GNa (pink trace) and GKf (blue trace) conductances. Middle, Activation of PHVA (orange trace) and PMVA (purple trace) permeabilities. Right, Activation of GKs (red trace), GBK (green trace), and GSK (dark blue trace) conductances. c, The activation of seven ion channels inserted into the model during the BAC firing simulated in Figure 5d at 470 μm from the soma. Left, Activation of GNa (pink trace) and GKf (blue trace) conductances. Middle, Activation of PHVA (orange trace) and PMVA (purple trace) permeabilities. Right, Activation of GKs (red trace), GBK (green trace), and GSK (dark blue trace) conductances.
Strong dendritic spikes are known to induce AP firing at the axon initial segment (Larkum et al., 1999a,b, 2001). The model was able to faithfully predict this effect, which was simulated by injecting a large EPSP-like current into the distal apical dendrite (Fig. 8e, red). Since the Ca2+ spike generated at the distal dendrite lacked the backpropagating Na+ spike, it forward propagated toward the soma. The mechanism was essentially the same as that in Figures 4, 5, 6, and 8. First, the GNa and the GKf conductances contributed to the changes in the membrane potential (Fig. 11b,c, left), but at the distal dendrite (600 μm) these conductances did not participate in the Na+ spike, as mentioned above (Fig. 11b, left). In contrast, at the more proximal dendrite (470 μm) these conductances partially contributed to the first spike (Fig. 11c, left). Next, the two Ca2+ permeabilities continued to depolarize the dendritic membrane at both locations (Fig. 11b,c, middle); these two permeabilities had a greater impact on the membrane depolarization at the distal dendrite (600 μm), where the Ca2+ spike occurred (Fig. 11b, middle). Then, the dendritic membrane repolarized as a result of the slow activation of the GKs, GBK, and GSK conductances (Fig. 11b,c, right), but the GKs conductance had a greater effect on the repolarization at 600 μm (Fig. 11b, right).
Conductance and permeability activation during a complex Ca2+ spike. a, Reconstruction of a filled L5 pyramidal neuron stained with biocytin, illustrating electrode placement. b, The activation of seven ion channels inserted into the model during the dendritic Ca2+ spike simulated in Figure 5e at 600 μm from the soma. Left, Activation of GNa (pink trace) and GKf (blue trace) conductances. Middle, Activation of PHVA (orange trace) and PMVA (purple trace) permeabilities. Right, Activation of GKs (red trace), GBK (green trace), and GSK (dark blue trace) conductances. c, The activation of seven ion channels inserted into the model during the dendritic Ca2+ spike simulated in Figure 5e at 470 μm from the soma. Left, Activation of GNa (pink trace) and GKf (blue trace) conductances. Middle, Activation of PHVA (orange trace) and PMVA (purple trace) permeabilities. Right, Activation of GKs (red trace), GBK (green trace), and GSK (dark blue trace) conductances.
The proximal apical dendrite modulates the coupling between the axonal AP initiation zone and the Ca2+ spike initiation zone in the distal apical dendrite (Larkum et al., 2001). Depolarization of the proximal apical dendrite increases coupling, whereas hyperpolarization decreases it (Larkum et al., 2001). To test the ability of our model to predict this modulation, we injected current into the proximal apical dendrite of our optimized model to test whether this zone could affect the propagation of a distal dendritic spike. We injected an EPSP-like current of 2 nA into the distal apical dendrite (at 800 μm; Fig. 12b1, gray). The dendritic potential simulated at the distal apical dendrite (Fig. 12b1, red) failed to propagate toward the soma (Fig. 12b1, blue and black). A small depolarization (0.2 nA, 50 ms, onset 30 ms before the EPSP-like waveform) at the proximal apical dendrite (at 300 μm; Fig. 12b2, blue pipette) was able to simulate a Ca2+ spike at the distal dendrite (Fig. 12b2, red), followed by two somatic APs (Fig. 12b2, black). Alternatively, a somatic AP elicited after injection of an EPSP-like waveform at the distal dendrite (at 800 μm; Fig. 12c1, gray and black) could be blocked (Fig. 12c2, black) by the injection of a hyperpolarized current (−0.4 nA, 50 ms, onset 30 ms before the EPSP-like waveform; Fig. 12c2, blue) at the proximal dendrite.
Modulation of the dendritic spike propagation. a, Reconstruction of a filled L5 pyramidal neuron stained with biocytin, illustrating electrode placement. b, A dendritic spike changed to BAC firing. b1, An EPSP-like current of 2 nA (rising τ = 2 ms, declining τ = 10 ms) was injected through the gray pipette (800 μm) into a dendrite of one of the optimized models (Table 1, cell 5). The simulated voltage response showed the shape of an EPSP at the soma (black trace) and at the apical dendrite (300 μm, blue trace; 600 μm, red trace). b2, A small depolarizing current at the proximal dendrite (300 μm, blue electrode) was injected (0.2 nA, 50 ms, onset 30 ms before the EPSP-like current). The dendritic spike converted to a Ca2+ spike at the distal dendrite (red trace). c, A dendritic AP changed to a dendritic spike. c1, An EPSP-like current injection of 3 nA (rising τ = 2 ms, declining τ = 10 ms) through the gray pipette (800 μm, bottom) generated a dendritic potential at the distal apical dendrite (red trace) that propagated to the soma (black trace). c2, Hyperpolarizing current (−0.4 nA, 50 ms, onset 30 ms before the EPSP-like current) was injected into the proximal dendrite (300 μm, blue electrode). The hyperpolarizing current caused the dendritic AP to convert to a dendritic spike (red trace).
Previous studies have shown that within a compartment, the neuronal inputs are summed nonlinearly, whereas the output of different compartments sums linearly (Polsky et al., 2004; Larkum et al., 2009). We used the optimized model to simulate NMDA spikes at thin dendrite branches and to examine the information integration. We inserted NMDA and AMPA synapses (Larkum et al., 2009) in the model within the same dendritic branch (the distance between the pipettes was 26 μm; Fig. 13b, left) and between two branches (Fig. 13b, right). First, we activated the synapses separately and then simultaneously (Fig. 13b). Figure 13c compares the EPSP integral of the combined stimulus (Fig. 13b, red) with the expected EPSP integral (summation of the two individual stimuli; Fig. 13b, blue). The input integration within the same branch was linear for a weak stimulus and became nonlinear as the stimulus intensity increased (Fig. 13b, left, c, red). By contrast, the summation of the inputs between branches remained linear regardless of the stimulus intensity (Fig. 13b, right, c, blue).
Simulated NMDA spikes within and between dendrite branches. a, Reconstruction of a layer 5 pyramidal neuron stained with biocytin, illustrating electrode placement. b, A NMDA-like spike was evoked through two distal pipettes (green and orange) after inserting NMDA and AMPA synapses (Larkum et al., 2009) into one of the optimized models (Table 1, cell 5). Left, Integration within branches. First, the distal pipettes (green and orange) were activated separately (black traces) and then simultaneously (red traces). The EPSPs were simulated 150 μm below the orange pipette. The blue trace represents the summation of the two individual activations. The summation is shown for weak, medium, and strong stimuli. The distance between the two pipettes at the same branch was 26 μm. Right, Integration between branches. First, the distal pipettes (gray and orange) were activated separately (black traces) and then simultaneously (red traces). The EPSPs were simulated 150 μm below the orange pipette. The summation is shown for weak, medium, and strong stimuli. c, Expected versus actual simulated EPSP intervals are plotted for a range of stimulus intensities for summation within the branch (red line) and between branches (blue line).
Discussion
We applied an adaptation of the Hodgkin–Huxley methodology to derive a mechanistic model for the apical dendrite of neocortical L5 pyramidal neurons. Using experimentally known channel kinetics, a pharmacological peeling procedure (Figs. 1, 2), and a genetic optimization algorithm, we were able to reduce the parameter space to manageable segments. The model predicted an ionic mechanism for the generation of dendritic spike (Figs. 3, 4, 5), the interaction of this spike with the backpropagating AP (Figs. 4, 8, 9), the ability of the proximal apical dendrite to control the coupling between the axon and the dendrite (Fig. 12), and the generation of NMDA spikes in the distal apical tuft (Fig. 13).
Our approach allowed us to constrain a mechanistic model for L5 pyramidal neurons containing passive parameters, and eight ion conductances and permeabilities known to be expressed in the membrane of this neuron (Stuart and Sakmann, 1994; Berger et al., 2001; Larkum et al., 2001; Benhassine and Berger, 2005; Schaefer et al., 2007; Almog and Korngreen, 2009; Deister et al., 2009). The passive parameters and the gradient of the dendritic distribution of fast and slow voltage-gated K+ channels (Kf and Ks), the voltage-gated Na+ channel and the Ih evolved to resemble previous experimental reports (Stuart and Sakmann, 1994; Berger et al., 2001; Larkum et al., 2001; Schaefer et al., 2007; Keren et al., 2009). Additionally, the mean dendritic density gradient of the GSK decreased along the apical dendrite. The GBK showed a uniform distribution, as was also found along the somatodendritic axis in the few previous studies of this channel (Benhassine and Berger, 2005). Both Ca2+ channels (HVA and MVA) showed a mean dendritic density gradient that decreased linearly along the apical dendrite (Fig. 3). These ion gradients along the somatodendritic axis in the model successfully simulated many dendritic electrophysiological events observed experimentally in L5 pyramidal neurons: AP backpropagation (Stuart and Sakmann, 1994; Stuart et al., 1997); regenerative dendritic Na+ and Ca2+ spikes (Schiller et al., 1997; Larkum et al., 1999a,b); BAC firing (Larkum et al., 1999a,b, 2001); changes in intracellular Ca2+ concentration (Schiller et al., 1995; Larkum et al., 1999b; Larkum and Zhu, 2002); and synaptic integration (Polsky et al., 2004; Larkum et al., 2009).
The optimized model predicted a mechanism of regenerative dendritic spikes (Fig. 4), and showed that a regenerative dendritic spike arises through a multistep mechanism, as follows: (1) Na+ and Kf channels contribute to the Na+ spike generated during spike onset, and due to the different kinetics of the channels, the Kf channel is activated before the Na+ channel for subthreshold events and vice versa for suprathreshold events; (2) the two Ca2+ channels, mainly MVA Ca2+ channels, contribute to the continuing depolarization; and (3) the remaining K+ channels (voltage- and Ca2+-gated) are responsible for membrane repolarization.
The optimized model predicted that the mechanism for the generation of dendritic spikes would be expressed throughout the apical dendrite in an almost uniform manner. The density of most voltage-gated channels included in the model leveled to a constant value approximately half way along the apical dendrite, with some channels leveling off even proximal to the soma (Fig. 3). This suggested that it should be possible to initiate a dendritic spike in any location along the apical dendrite. This appeared at first to contradict experimental observations showing that a dendritic spike could only be triggered in the distal apical dendrite (Amitai et al., 1993; Schiller et al., 1997; Larkum et al., 1999a, 2001; Pérez-Garci et al., 2013). This can be accounted for by considering the threshold for spike initiation in the axon initial segment (Fig. 4c). Specifically, the proximal dendritic current injection passively propagates to the axon initial segment and generates an AP before the local dendritic depolarization crosses the threshold for the generation of a dendritic spike. Furthermore, the increase in the threshold of the dendritic spike proximal to the soma, probably due to the passive properties of the proximal apical dendrite, provides a computational model for the ability of this region to control the coupling between the distal dendrite and the axon. Depolarization of the proximal dendrite allows the forward-propagating dendritic spike to cross the local threshold for dendritic spike generation, thus facilitating active forward propagation of this spike in the direction of the soma (compare Fig. 12 with Larkum et al., 2001). It is important to note that the interactions of the dendritic spike with the axonally generated AP have been investigated by attaching an artificial axon to the optimized soma (Mainen and Sejnowski, 1996). This artificial axon was used since we did not have experimental recordings from the axon that could be used to optimize a model for that cellular section (Keren et al., 2005). Thus, the results presented in Figures 4c, 8, 9, and 12 should be considered as qualitative alone. Nevertheless, these simulations provide computational predictions regarding the mechanism controlling the interaction between dendritic and axonal spikes.
Our current model predicts a different mechanism for the initiation of dendritic spikes from those previously suggested by other numerical simulations (Schaefer et al., 2003; Larkum et al., 2009; Hay et al., 2011). Previous models of dendritic spiking in L5 pyramidal neurons assumed the existence of a dendritic hot zone of low-threshold activated (t-type), voltage-gated Ca2+ channels in the distal apical dendrite (Schaefer et al., 2003; Larkum et al., 2009; Hay et al., 2011). These are the two major differences between our model and these previous studies. Our model does not contain a t-type Ca2+ channel because we did not observe this type of current in our somatic nucleated patch recordings (Almog and Korngreen, 2009). Our model includes an MVA channel whose kinetics were largely based on these recordings. This was the primary reason for our lack of success in optimizing a working model when we assumed that a hot zone of MVA channels existed in the distal apical dendrite (Fig. 3). MVA channels require a larger depolarization than the t-type channels to participate in a regenerative dendritic event. Thus, a hot zone of these channels was not hot enough to generate a Ca2+ spike.
We have previously suggested an optimization scheme for a full compartmental model of a neuron using multiple voltage recordings from the soma, axon, and dendrite as a training dataset for the genetic algorithm (Keren et al., 2005). We have also shown that the point-by-point sum of squares is sufficient as a cost function to constrain a model for a single dendrite, as long as the somatic membrane potential is used as a voltage-clamp command and a parameter peeling procedure is applied (Keren et al., 2009). An intriguing alternative to the peeling procedure requires the sampling of the membrane potential from many locations along the soma, dendrites, and axon (Huys et al., 2006). Unfortunately, such massive parallel recordings of the membrane potential requires considerably better voltage-sensitive dyes than we have at present. The results presented here suggest that finding large parameter manifolds is partly due to optimizing the model by only using the averages of extracted AP and dendritic spike features as training datasets (Hay et al., 2011). It is clear that to optimize a model for a cellular region displaying regenerative activity, it is imperative to obtain membrane potential recordings from that region. Even with this expanded dataset, pharmacological parameter peeling and multistage optimization are required to avoid local minima.
The contradiction between our linear model and the previously suggested hot zone model cannot be resolved computationally. Hot zone models have predicted many physiological properties of L5 pyramidal neurons (Schaefer et al., 2003; Larkum et al., 2009; Hay et al., 2011). The best way to resolve this issue would be to experimentally determine, probably using the patch-clamp technique, the dendritic density gradient of voltage-gated Ca2+ channels. This experimental verification of our simulation will also serve to further expand the current model that only contains two voltage-gated Ca2+ channels. Sadly, although several valiant but fruitless attempts have been made, these dendritic Ca2+ gradients have yet to be determined.
The simulations illustrated in Figure 4c suggest an experiment that may be used to support the model we presented. The simulations predict that it would be possible to trigger a dendritic spike even in the proximal apical dendrite following a complete block of the voltage-gated Na+ channels in the axon initial segment. Thus, a possible experiment would entail a simultaneous recording from the soma and various distances along the apical dendrite and an application electrode puffing tetrodotoxin at the axon initial segment. Should such an experiment reproduce the predicted spike threshold profile (Fig. 4c), it would strengthen the predictions of our model. Finally, applying the peeling procedure to quantitative modeling of the axon initial segment would doubtless lead to a more comprehensive modeling of L5 pyramidal neurons and other cortical axons.
Footnotes
- Received July 8, 2013.
- Revision received October 21, 2013.
- Accepted October 28, 2013.
This work was supported by Grant 212/07 from the Israel Science Foundation.
The authors declare no competing financial interests.
- Correspondence should be addressed to Dr. Alon Korngreen, Mina and Everard Goodman Faculty of Life Sciences, Leslie and Susan Gonda Multidisciplinary Brain Research Center, Bar-Ilan University, Ramat Gan 52900, Israel. alon.korngreen{at}biu.ac.il
- Copyright © 2014 the authors 0270-6474/14/340182-15$15.00/0