The detailed model of spontaneous Ca²⁺ spike generation was derived from the original biophysical model (Schmidt et al., 2003; Anwar et al., 2012) developed in the NEURON simulator (Hines and Carnevale, 1997). The model contained four types of ion channels: P-type Ca²⁺ channel (P_max = 2 × 10⁻⁴ cm/s, GHK equation) (Swensen and Bean, 2005), T-type Ca²⁺ channel (P_max = 8 × 10⁻⁶ cm/s, GHK equation) (Iftinca et al., 2006), BK-type Ca²⁺-activated K⁺ channel (G_max = 7 × 10⁻² S/cm²) (Cox et al., 1997) and SK-type Ca²⁺-activated K⁺ channel (G_max = 3.1 × 10⁻⁴ S/cm²) (Hirschberg et al., 1998; Solinas et al., 2007), plus a leak channel (G_max = 1 × 10⁻⁶ S/cm² and E_rev = −61 mV).

Intracellular Ca²⁺ was modeled using the following Ca²⁺ buffering mechanisms.

All the simulations were run using a time step of 0.02 ms. Model scripts for all models used in this work are available at http://senselab.med.yale.edu/modeldb/ShowModel.asp?model=155731.

Many multi-compartment dendritic models use a single pool model of Ca²⁺ buffering, which simulates only the submembrane Ca²⁺ concentration to control K_Ca channels. These phenomenological models convert Ca²⁺ current passing through VGCC to Ca²⁺ concentration using a submembrane shell of fixed depth, d [Equation (1); Figure 3A]. The volume of a submembrane shell (as for SP_new) is defined as:

where Vol_f denotes the volume of a full compartment, Vol_c denotes the volume of the core, d is the depth of submembrane shell, diam is the diameter of compartment and L is its length. Using such a representation of submembrane shell, SVR_new equals

where SA is the surface area of the compartment. Note that SVR_new is less dependent on diameter than the SVR for the complete volume (1/diam), but, as shown in Figure 3A, SVR_new still increases for smaller diameters.

However, we noticed that most of the Ca²⁺ shell models implemented in NEURON use an incorrect volume for the submembrane shell (Figure 3B). In those models (SP_old), the volume of a submembrane shell is defined as:

Using such a representation of submembrane shell, gives SVR_old

Because the depth (d) of submembrane shells is usually taken constant for variable diameter compartments, SVR_old is constant and independent of the dendrite diameter. As a result, we observe changes in peak amplitudes of the simulated Ca²⁺ transients using SP_new in compartments with different diameters (Figure 3C) while the same compartments with SP_old always show exactly the same Ca²⁺ transient.

Since pool based models are phenomenological models, the values of depth (d) and decay time constants (β) can be tuned to approximate the desired behavior of intracellular [Ca²⁺]. In Figure 3D we show how the error of using SP_old (compared to SP_new) depends on the values of d and β used. The errors in peak Ca²⁺ were computed using β-values of 0.02 ms⁻¹ (Traub and Llinas, 1977), 6.86 ms⁻¹ (Anwar et al., 2012), and 10 ms⁻¹ (De Schutter and Bower, 1994) and using submembrane shells with depths ranging from 0.05 to 0.25 μm. The error increases with the size of depth used, as expected from SVR_new. More importantly, these errors become significantly larger for smaller diameters (diam <1 μm) and may reach up to 80% for 0.1 μm diameter compartments (inset of Figure 3D). Typically, distal dendrites have large numbers of dendritic branches with diameters less than 1 μm.

In the rest of our study we will only focus on SP_new and DM to investigate how well they can capture Ca²⁺ gradients in dendrites.

Figure 2 demonstrates in one dendritic morphology that SVR differences cause sharper Ca²⁺ gradients when it was modeled using DM compared to SP_new. We next investigated whether this is a systematic observation by simulating the dendritic Ca²⁺ spike model in 11 additional dendritic reconstructions of PCs using both methods and comparing the results (Figure 4). For each PC, the ion channels were distributed uniformly on its dendrite and each unbranched segment had a constant diameter.

Because this study mostly focuses on local differences in Ca²⁺ concentration, we summarized the data on spatial gradients by computing the ratio of integrated [Ca²⁺] in adjacent dendritic segments and plotting the distributions of these ratios in Figure 4. We observe a wide range of distributions of spatial fluctuations of Ca²⁺ levels in different neuron reconstructions, with the histograms of some neurons (e.g., PC 2 and 5) showing very large tails and other ones only small fluctuations (ratio < 2). This observation may be related to differences in the quality of the reconstructions (see Discussion). But in all cases, the DM model always results in a wider range of Ca²⁺ fluctuations as compared to SP_new.

In Figure 4 we used a common approach in compartmental modeling: we ignored small variations in diameter by taking only a single compartment for each unbranched segment. However, due to the large jumps in Ca²⁺ concentrations between neighboring compartments observed in some parts of the model (Figures 2, 4) we wondered about the realism of this assumption of uniform diameter. We investigated this issue both in the 12 PCs modeled previously as well as in 284 neocortical and 38 hippocampal pyramidal neuron reconstructions, because larger changes in diameter may be present in morphological classes where the level of branching is not as extensive as in PCs. For both neuron types we computed the coefficient of variation (CV) of diameters for every dendritic segment (between two branch points) based on all the measurements available in the morphological reconstruction (Figure 5). We observed a large variability in CV of reconstructed morphologies of neurons obtained from different laboratories for both cell classes (see Discussion), but overall the variability of diameter was much larger in pyramidal neurons where in many cells more than a quarter of the unbranched segments had CVs of 0.4 or more. In PCs more than half of the reconstructions had CVs of 0.2 or more in at least a quarter of their unbranched segments.

Our analysis suggests that for pyramidal neurons even more care should be taken when converting dendritic segments into cylindrical compartments. A good representation of dendritic segments with rapidly varying diameters is essential to model intracellular Ca²⁺ concentrations correctly (see also Figure 8).

Though the implementation of 1D diffusion in concentric cylindrical shells may seem straightforward, the NMODL language used in NEURON actually makes it difficult to do this in a flexible way and we discovered that many existing models do not implement it correctly. The standard example in the NEURON book (Carnevale and Hines, 2006) is a variable depth scheme where the volume is divided over a fixed number of concentric shells (4 in the standard example) with variable depth Equations (3) and (4) and Figure 1: DM, note that the submembrane and core shells have a smaller depth than the others). Many models using NEURON implement exactly this mechanism: 4 shells and all with variable depth. Because the volume of the submembrane shell is used to convert inward Ca²⁺ currents into a Ca²⁺ concentration that directly activates K_Ca channels, varying its depth will affect the computed value of this Ca²⁺ concentration. From a biophysical perspective there is no reason why the depth of a submembrane shell that is assumed to simulate the effective volume affecting the Ca²⁺ sensors of K_Ca channels (Fakler and Adelman, 2008; De Schutter, 2010) should vary greatly with dendrite diameter. We will therefore consider two issues: the number of shells to be modeled and a submembrane shell with variable (DM) or fixed depth (DM_FD).

If one wants to vary the number of shells modeled depending on compartment diameter, which is the correct solution, a separate NEURON mechanism has to be created for each specific number of shells that is required. Some authors have therefore decided to use a fixed number of shells with variable depth of each shell (Migliore et al., 1995; Lazarewicz et al., 2002; Gold et al., 2007; Lavzin et al., 2012) (Figure 1: DM), but this can lead to significant errors in simulated submembrane Ca²⁺ concentration in large diameter dendrites if the number of shells is taken to be small (Figure 6A: circles). These errors show both a positive and negative component depending on compartment diameter, suggesting that two types of error contribute. Indeed, when we repeated these simulations with a fixed submembrane shell depth d₁ of 0.1 μm and the rest of the volume divided over the remaining shells with equal, variable depths (FD: Figure 6B: triangles) only a positive error, increasing with diameter, remains. Because this error is quite small for a large number of shells, a model with the same large number of shells in every compartment will give accurate results in NEURON, but this may cause unacceptably slow runtimes (Anwar et al., 2012) so it is better to vary the number of shells [DM_FD mechanism, Equation (5)].

The next question is then how to compute shell depth as the depth of at least one shell has to vary to fit the total exactly to a variable compartment diameter. As already mentioned, in the standard NEURON implementation (Carnevale and Hines, 2006) the depths of all shells vary with compartment diameter [variable depth scheme; Equations (3) and (4) and Figure 1: DM], including that of the submembrane shell. In effect, the depth of the submembrane shell (d₁) may vary between d − 0.25 d and d + 0.25 d. In Figure 6C, the broken line shows the theoretical error of Ca²⁺ influx conversion to Ca²⁺ concentration using the variable depth scheme (range of d₁ due to discretization: 0.075–0.125 μm). The larger predicted errors in these cases are associated with small diameters, where small changes of diam will result in bigger changes in d₁ Equation (4) and submembrane shell volume.

These errors are large and should not be ignored. But what is the effect of these geometrical errors on actual computed Ca²⁺ concentrations? To quantify this we simulated Ca²⁺ transients using a mechanism with a variable number of shells, all with the same depth of 0.1 μm except for the core shell which has a variable diameter [DM_FD, Equation (5) and Figure 1], which is assumed to give the most accurate solution. We found that DM and DM_FD show very similar peak amplitudes and decay time constants for different diameter compartments (Figure 6B), resulting in much smaller errors in peak amplitudes of Ca²⁺ using DM compared to DM_FD than theoretically predicted (Figure 6C). The error depends on the size of Ca²⁺ influx in a nontrivial way, but for all levels of Ca²⁺ influx it was small with the largest error only about 4%. This significant difference with the theoretical prediction is due to strong buffering (especially in PCs; (Hartmann and Konnerth, 2005) and diffusion of Ca²⁺, which removes most of Ca²⁺ entering into the submembrane shell.

How do the different Ca²⁺ buffering models respond to local fluctuations of dendrite diameter in terms of Ca²⁺ levels? To estimate the effect of dendritic diameter changes on Ca²⁺ dynamics using each model, we computed the ratio of integrated Ca²⁺ transients measured in each of the pair of simulated compartments using SP_new, DM and DM_FD for many possible pairs of dendritic diameters (range: 0.1–6.0 μm with increments of 0.1 μm) (Figure 7). As explained previously, the SP_old model does not show any sensitivity to changes in diameters. For SP_new, large ratios (>2) are limited to combinations where a compartment with an extremely small diameter (≤0.3 μm) is connected to one with large diameters. For DM and DM_FD, this region expands to all compartments with diameter less than or equal to 1 μm that are connected to ones with larger diameters. Therefore, the detailed Ca²⁺ dynamics models are more sensitive to changes in dendritic diameter compared to pool based models, which explains the differences observed in Figure 3. But, although the sensitivity maps of DM are noisier than those of DM_FD, due to the use of variable depth submembrane shells in DM, overall these maps are quite similar to each other.

We conclude from Figures 6, 7 that correct simulation of radial 1D Ca²⁺ diffusion requires a variable number of concentric shells that scales with compartment diameter, but that the Ca²⁺ dynamics are less sensitive to the actual scheme used to compute the depth of these shells.

The predictions of Figure 7 are based on no (SP_new) or only radial 1D diffusion (DM and DM_FD). In this case, would the predicted large Ca²⁺ transients disappear in the presence of 3D diffusion? Also, what happens if dendritic diameter varies more smoothly than possible in a model using electrical compartmentalization? To address both issues, we used the STEPS simulator (Hepburn et al., 2012), which uses tetrahedral meshes to accurately represent detailed morphologies and 3D diffusion of molecules to simulate Ca²⁺ dynamics.

When we simulated Ca²⁺ transients with the detailed Ca²⁺ dynamics model in part of a PC dendritic arbor using STEPS, we still observed large fluctuations in Ca²⁺ levels along the different branches. Figure 8A shows the integrated Ca²⁺ levels for each tetrahedron located within 0.1 μm from the membrane. The large fluctuations of Ca²⁺ levels appear to be related to dendrite diameter. High Ca²⁺ levels are observed in dendritic regions with small diameter (Figure 8B) and at the tips of terminating branches. Higher levels at the tips are due to the higher SVR as a result of their small diameters and the reduced effective diffusion because of the closed end condition (the latter is not predicted by radial diffusion models). Overall we conclude that neither 3D diffusion nor smooth changes in dendrite diameter reduce the pronounced Ca²⁺ gradients caused by variable dendrite diameter, raising the question what level of detail is necessary to model this effect correctly?

To address this question, we compare NEURON and STEPS simulations using two different compartmentalization schemes in NEURON. First, we simulated Ca²⁺ transients in NEURON using multiple compartments per unbranched segment to capture all changes in dendritic diameters (Figure 8D; right panel). For comparison Figure 8D (left panel) shows the STEPS simulation with mean integrated Ca²⁺ concentration computed for all tetrahedrons corresponding to every NEURON compartment. Next, we made a similar comparison with NEURON simulations where every unbranched dendritic segment is considered as a single compartment (Figure 8F), which is the approach used in many compartmental models. Comparing these spatial maps (Figures 8D,F), we observe only small differences between simulations with 1D diffusion (NEURON) or with 3D diffusion (STEPS). However, the actual Ca²⁺ levels are different in the respective simulations. To quantify the difference in Ca²⁺ levels between the two approaches and how they relate to fluctuations in dendritic diameters we computed the ratios of Ca²⁺ levels and diameters for all adjacent segments. Figure 8E shows Ca²⁺ ratios in adjacent compartments for small compartment sizes (data shown in Figure 8D) and Figure 8G shows Ca²⁺ ratios in adjacent compartments with one compartment per unbranched dendritic segment (data shown in Figure 8F). The comparison of Figure 8E with Figure 8G clearly shows that the use of large compartments will result in larger jumps in Ca²⁺ levels between adjacent compartments. Using many small compartments to capture the continuous change of dendritic diameters results in much smoother and smaller changes in Ca²⁺ levels. The overall behavior of these changes in Ca²⁺ levels (Figures 8E,G) is similar, respectively, to the ratios of diameters in the original morphological reconstruction and to the ratios for adjacent compartments diameters for one compartment per unbranched segment (Figure 8C). This confirms that the simulated Ca²⁺ gradients are largely caused by the SVR effect. Finally, notice that the effect of 3D diffusion is more prominent when using small compartments (Figure 8D, bigger difference between NEURON and STEPS simulation).

The complexity of dendritic geometry and structure has been studied extensively to investigate its effects on propagation of action potentials, its role in synaptic efficacy and its effects on limiting interaction across different active dendritic regions (Mainen and Sejnowski, 1996; Vetter et al., 2001). Although variable levels of Ca²⁺ in different dendritic regions have been reported previously (Tank et al., 1988; Lev-Ram et al., 1992; Schiller et al., 1995), only a few studies specifically investigated the effect of dendritic diameters on Ca²⁺ dynamics (Holthoff et al., 2002; Rozsa et al., 2004). Therefore, those effects are often omitted while constructing biophysical models of dendrites. Due to limited quantitative information about the mechanisms controlling Ca²⁺ levels in many neurons, phenomenological models of Ca²⁺ buffering, such as the single exponential decaying pool, are commonly used for biophysical neuronal modeling. Such models when used correctly capture only some aspects of the highly complex behavior of intracellular Ca²⁺ buffering dynamics. In our previous work (Anwar et al., 2012), we showed that pool based models of Ca²⁺ buffering fail to correctly predict peak Ca²⁺ concentrations and decay time constants important for the interaction between VGCC and K_Ca channels. In this study, we investigated the effect of dendritic diameters on Ca²⁺ dynamics using a modeling approach. Our results (Figures 2, 4) show that pool based models have limited ability to capture the spatial variability of Ca²⁺ dynamics in morphologically complex dendrites as compared to a detailed Ca²⁺ dynamics model with radial 1D diffusion. The detailed Ca²⁺ dynamics model shows different peak amplitudes of Ca²⁺ levels as well as different (and multiple) decay time constants (Figure 6B). In contrast, pool based models with correct submembrane volume only show different peak amplitudes of Ca²⁺ levels (Figure 3C).

In general, many studies, including this one, ignore additional properties of real neurons that will affect Ca²⁺ dynamics. The most important of these simplifications are the assumption of constant density of Ca²⁺ channels, which is known to be not true for many neurons (for review see Johnston et al., 1996; Migliore and Shepherd, 2002), and the omission of the effect of organelles in the cytoplasm that block diffusion and have additional membrane Ca²⁺ pumps (mainly endoplasmic reticulum and mitochondria). Another important determinant of Ca²⁺ dynamics is inhomogeneous distribution of Ca²⁺ buffers in dendrites of a given neuron, as well as their properties, causing competitive binding of Ca²⁺ to available Ca²⁺ buffers and Ca²⁺ pumps (Markram et al., 1998). While the density of channels can easily be changed in compartmental models, accurate representation of intracellular organelles is possible in mesh based models only. We do not expect that inclusion of these properties would significantly change our conclusions.

Dendrites have variable diameters. Typically, the diameter of dendrites taper with increasing distance from the soma. It is generally assumed that the change in diameter of an unbranched dendrite is relatively small as compared to the change in diameter at branching, which allows representation of an unbranched dendrite segment as a single uniform diameter compartment. However, as shown in our morphological analysis in Figure 5, many reconstructions of both PCs and pyramidal neurons show great diameter variability within their unbranched segments, with CV values sometimes reaching above 0.5. This implies that a correct Ca²⁺ dynamics model should represent this diameter variability by having several compartments for each unbranched segment (Figures 8D–E), but also that the quality of the morphological reconstruction is of utmost importance. We observed great differences of the diameter variability between different neural reconstructions which often could be related to the laboratory where the reconstructions have been done, as was reported previously for pyramidal neurons (Scorcioni et al., 2004; Szilagyi and De Schutter, 2004; Holmes et al., 2006). Because it seems more likely that human error causes an undersampling of diameter changes than an exaggeration, we assume that the reconstruction with high diameter CV tend to be more reliable. Finally, one should be aware that software like CVapp (Cannon et al., 1998), which converts morphology files into formats suitable for NEURON simulation, uses a specific discretization scheme that changes diameters at branch points (Figure S1).

Although morphological reconstructions obtained using electron microscopy (EM) capture dendrites much more precisely, because of rapid fluctuations in dendritic surface those reconstructions are not suitable for compartment based models. The proper use of EM reconstructions in modeling Ca²⁺ dynamics will require more advanced simulators with support for surface or tetrahedral meshes (e.g., M-Cell, STEPS). Also, this will require more detailed description of Ca²⁺ related mechanisms (e.g., spatial distribution of VGCa channels, K_Ca channels, buffers, pumps and internal calcium stores).

Almost all biophysically detailed models have been constructed using either the GENESIS or NEURON simulators. These software packages are based on compartmentalization of dendritic structures into multiple iso-potential cylinders, where voltage, currents and concentrations are computed for each of those compartments independently. Since these compartments are based on electrical properties of dendrites, biochemical representation of intracellular mechanisms in these compartments is always an approximation of the related biophysical process. Such a simplified molecular representation may result in unrealistic behavior of the model, depending on the rationale behind the assumption and its accuracy. One such example is the commonly used single pool model to simulate intracellular Ca²⁺ in the NEURON simulator Equations (8) and (9). The conversion of influx to intracellular Ca²⁺ concentration in these models is incorrect, which is due to the use of an inaccurate volume of the submembrane shell (Figure 3B). Although this inaccuracy does not influence the results of single compartment models dramatically because it can be easily tuned by adapting the shell depth, it becomes critically important in multi-compartment models. This incorrect single pool model will always underestimate the influx (Figure 3C), which may require unrealistic distribution of dendritic VGCC and K_Ca channels during model construction and cause a mismatch in input resistance between model and actual cell, and it will not predict any spatial gradients of Ca²⁺ concentration due to fluctuations of dendrite diameter (Figure 2A).

More detailed Ca²⁺ dynamics models using radial 1D diffusion are thought to be more accurate, but again the compartmentalization of the dendrite may result in either an inaccurate or incomplete representation of model. It is a major challenge to model diffusion in the NEURON simulator correctly. NEURON allows radial (toward the center of each dendritic compartment) and longitudinal (from one compartment to neighboring compartment) diffusion. Radial diffusion requires virtual submembrane shells (Figure 1), where shells typically have a variable depth, depending on the diameter of each compartment. Furthermore, longitudinal diffusion is only allowed if the adjacent compartments have the same number of shells which will introduce a larger error (Figure 6A) unless a very large number of shells is used everywhere. Conversely, though theoretically the variable submembrane depths of the standard NEURON scheme (DM) should result in large errors, this effect was strongly filtered by the diffusion and buffering mechanisms, resulting in only small differences (Figure 6C) with a method (DM_FD) that ensures a fixed depth of the submembrane shell. It should be noted, however, that these differences may be larger in models of other neurons because the buffering capacity of PCs is exceptionally high (Hartmann and Konnerth, 2005). Control simulations showed that although the extent of the changes in Ca²⁺ levels varied in models with lower buffers concentrations, the dependence of Ca²⁺ levels on changes in diameters persisted (results not shown).

Neither of the issues just mentioned are relevant for the GENESIS simulator. Both submembrane pools (as the concpool object, De Schutter and Bower, 1994) and radial diffusion (as the difshell object, De Schutter and Smolen, 1999) are implemented correctly and are easy to set up. Conversely, it is time consuming to create multiple calcium dynamics mechanisms with radial diffusion in NEURON because a separate mechanism has to be written for each set of diameters (see Materials and Methods) and this requires a lot of extra care. And then, for every different morphology, one will have to repeat the process. We expect that multilevel declarative model description languages (Raikov and De Schutter, 2012) may allow transparent and correct compartment based assignment of molecular mechanisms in NEURON in the future.

Even in PCs, where the estimated Ca²⁺ diffusion range is only about 5 μm (Santamaria et al., 2006), we observe effects of 3D diffusion on Ca²⁺ transients compared to when only radial 1D diffusion is used, especially when Ca²⁺ concentration is averaged over short distances only (Figures 8D,E). Nevertheless, the error introduced by the 1D approach is much smaller than the errors caused by inaccurate morphologies (Figure 8C) and simulating 3D diffusion in tetrahedral meshes is quite slow. However, 3D diffusion must be included in biophysically accurate models of synaptic plasticity or models involving Ca²⁺ based signaling pathways.

For most of modeling projects with the goal of capturing excitability and integrative properties of dendrites, a correct implementation of 1D radial diffusion and buffering in NEURON (or any other compartment based simulator) will be an adequate approximation. It is then important to implement a variable number of submembrane shells, with larger number of shells in larger diameter compartments, and best using a fixed depth of the submembrane shell [DM_FD model: Figure 1, Equation (5)]. The model should be based on a high quality morphological reconstruction (Jacobs et al., 2010) and the variability of diameter along dendritic segments should be retained by having as many compartments as required to capture diameter changes (Figures 8D,E).

Finally, we do not recommend the use of simple pool models, unless good data on the properties of Ca²⁺ buffering (e.g., Schmidt et al., 2003 for PCs) in the neuron type to be modeled is completely absent. If one is forced to use a simple pool model, make sure it is implemented correctly [SP_new, Figure 3A and Equations (2) and (6)].

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.