Controlling Ca²⁺-Activated K⁺ Channels with Models of Ca²⁺ Buffering in Purkinje Cells

Abstract

Intracellular Ca²⁺ concentrations play a crucial role in the physiological interaction between Ca²⁺ channels and Ca²⁺-activated K⁺ channels. The commonly used model, a Ca²⁺ pool with a short relaxation time, fails to simulate interactions occurring at multiple time scales. On the other hand, detailed computational models including various Ca²⁺ buffers and pumps can result in large computational cost due to radial diffusion in large compartments, which may be undesirable when simulating morphologically detailed Purkinje cell models. We present a method using a compensating mechanism to replace radial diffusion and compared the dynamics of different Ca²⁺ buffering models during generation of a dendritic Ca²⁺ spike in a single compartment model of a PC dendritic segment with Ca²⁺ channels of P- and T-type and Ca²⁺-activated K⁺ channels of BK- and SK-type. The Ca²⁺ dynamics models used are (1) a single Ca²⁺ pool; (2) two Ca²⁺ pools, respectively, for the fast and slow transients; (3) detailed Ca²⁺ dynamics with buffers, pump, and diffusion; and (4) detailed Ca²⁺ dynamics with buffers, pump, and diffusion compensation. Our results show that detailed Ca²⁺ dynamics models have significantly better control over Ca²⁺-activated K⁺ channels and lead to physiologically more realistic simulations of Ca²⁺ spikes and bursting. Furthermore, the compensating mechanism largely eliminates the effect of removing diffusion from the model on Ca²⁺ dynamics over multiple time scales.

Introduction

Purkinje cells (PCs) are known to express large conductance (BK) [1–10] and small conductance (SK) [6, 11–14] Ca²⁺-activated K⁺ (K_Ca) channels on their dendrites. K_Ca channels together with voltage-gated Ca²⁺ channels significantly control the dendritic excitability [6, 15]. Ca²⁺ entering through the voltage-gated Ca²⁺ channels causes the free cytosolic Ca²⁺ concentration to rise, which in turn controls the activation of K_Ca channels. A significant difference between those channels is their spatial separation. It has been shown that BK channels are in closer vicinity of Ca²⁺ sources compared to SK channels, suggesting that BK channels require a brief large amount (~10–100 μM) of Ca²⁺, while SK channels are activated by lower concentrations (~0–2 μM) [16–18]. Intracellular Ca²⁺ mechanisms like diffusion, endogenous buffers, internal stores, and pumps significantly control the cytosolic spread of Ca²⁺ and shape the Ca²⁺ signal that activates K_Ca channels [19]. Therefore, it is important to model intracellular Ca²⁺ dynamics carefully in a biophysically detailed PC model.

Existing PC models use a Ca²⁺ pool with a single relaxation time [20–22], which can reasonably approximate the fast transient [23] and can activate a BK channel well [17, 18]. On the other hand, it provides inappropriate Ca²⁺ signals to activate SK channels to support their role in excitability modulation. A simple and effective solution is to use two Ca²⁺ pools [24], respectively, for the fast and slow transient. A more comprehensive and biophysically realistic solution is to use a detailed model [25] with Ca²⁺ buffers, diffusion of those buffers, and diffusion of Ca²⁺ and Ca²⁺ extrusion pumps. Large compartments also require modeling of radial diffusion of Ca²⁺ ions and buffers [26]. The presence of diffusion in these detailed models can result in a large computational cost [27], which may need to be avoided in simulating morphologically detailed PC models.

In this study, we compare the effectiveness of different Ca²⁺ buffering models in controlling Ca²⁺-activated K⁺ channels. We build a single compartment model of a PC dendritic segment, including P-type and T-type voltage-gated Ca²⁺ channels and BK-type and SK-type Ca²⁺-activated K⁺ channels, and tune this to generate Ca²⁺ spikes. Intracellular Ca²⁺ dynamics are modeled using a single Ca²⁺ pool, two Ca²⁺ pools, or detailed Ca²⁺ dynamics with calbindin (CB), parvalbumin (PV), a pump, and diffusion. Our results show that a detailed Ca²⁺ dynamics model with buffers, pumps, and diffusion has better control over Ca²⁺-activated K⁺ channels and can generate physiologically more realistic Ca²⁺ spikes. We also introduce a method to specify a diffusion compensating mechanism (DCM) that can replace the radial diffusion in the model. Our result shows that the consequences of removing diffusion from the model on the simulated Ca²⁺ dynamics can be largely eliminated by these compensating mechanisms for both short and long time scales. This allows the generation of physiologically realistic Ca²⁺ spikes that can be simulated with significantly less computational cost.

Materials and Methods

PC Model for Ca²⁺ Spikes

A dendritic compartment with a diameter of 4 μm and a length of 20 μm was used as a PC dendrite model for simulating dendritic Ca²⁺ spikes. All the simulations were run in the NEURON simulation environment [28]. The details about the currents used in our model are given below.

P-Type Ca²⁺ Channel

The P-type calcium current was based on data from Swensen and Bean [29] and included three activation (m) gates. The current density was based on the Goldman–Hodgkin–Katz equation [30]. The equations describing its kinetics are summarized below:

$$ \begin{array}{*{20}{c}} {{I_{{{\rm{C}}{{\rm{a}}^{{2 + }}}}}} = \overline {{P_{{{\rm{C}}{{\rm{a}}^{{2 + }}}}}}} \times {m^3} \times {g_{\rm{GHK}}}} \hfill \\{{m_{\infty }} = \frac{1}{{1 + {e^{{ - \frac{{V + 24.758}}{{8.429}}}}}}}} \hfill \\{{\tau_m} = \left\{ {\begin{array}{*{20}{c}} {0.2702 + 1.1622{e^{{ - \frac{{{{\left( {V + 22.098} \right)}^2}}}{{164.19}}}}}} \hfill & {{\hbox{if }}V \geqslant - {40}\,{\hbox{mV,}}} \hfill \\{0.6923{e^{{\frac{{V - 4.7}}{{1,089.372}}}}}} \hfill & {\hbox{otherwise}} \hfill \\\end{array} } \right\}} \hfill \\\end{array} $$

The factor g _GHK is a current per unit permeability.

T-Type Ca²⁺ Channel

The low threshold voltage-activated, T-type calcium current was modeled using Cav 3.1 data from Iftinca et al. [31]. It included two activation (m) gates and an inactivation (h) gate. The current density is based on the Goldman–Hodgkin–Katz equation [30]. The equations describing its kinetics are summarized below:

$$ \begin{array}{*{20}{c}} {{I_{{{\rm{C}}{{\rm{a}}^{{2 + }}}}}} = \overline {{P_{{{\rm{C}}{{\rm{a}}^{{2 + }}}}}}} \times {m^2} \times h \times {g_{\rm{GHK}}}} \hfill \\{\begin{array}{*{20}{c}} {{m_{\infty }} = \frac{1}{{1 + {e^{{ - \frac{{V + 52}}{5}}}}}},} \hfill & {{\tau_m} = \left\{ {\begin{array}{*{20}{c}} 1 \hfill & {i{\hbox{f}}\,V\, \leqslant - 90\,{\hbox{mV}},} \hfill \\{1 + \frac{1}{{{e^{{\frac{{V + 40}}{9}}}} + {e^{{ - \frac{{V + 102}}{{18}}}}}}}} \hfill & {\hbox{otherwise}} \hfill \\\end{array} } \right\}} \hfill \\\end{array} } \hfill \\{\begin{array}{*{20}{c}} {{h_{\infty }} = \frac{1}{{1 + {e^{{\frac{{V + 72}}{7}}}}}},} \hfill & {{\tau_h} = 15 + \frac{1}{{{e^{{\frac{{V + 32}}{7}}}}}}} \hfill \\\end{array} } \hfill \\\end{array} $$

The factor g _GHK is a current per unit permeability.

BK-Type K_Ca Channel

The BK-type K_Ca channel was based on a kinetic scheme proposed by Cox et al. [32]. The channel model has a single voltage-dependent gate and four binding sites for Ca²⁺ ions, and Ca²⁺ binding at each site modulates the opening/closing rate coefficients. The kinetic scheme and rate constants used in the BK model were obtained from scheme II and patch 3 data in Table III in reference [32].

SK-Type K_Ca Channel

We used the SK-type channel from a Golgi cell model [33], which was based on data from Hirschberg et al. [34].

Leak Current

The leak current was modeled as a linear voltage-independent conductance following Ohm's law:

$$ {I_{\rm{leak}}} = {G_{\rm{leak}}}\left( {V - {E_{\rm{leak}}}} \right), $$

where ${G_{\rm{leak}}} = {10^{{ - 6}}}\,{\hbox{S/c}}{{\hbox{m}}^2}$ and ${E_{\rm{leak}}} = - 61\,{\hbox{mV}}$.

Intracellular and Extracellular Ca²⁺

Intracellular Ca²⁺ was modeled using different buffering mechanisms described in the following section. Extracellular Ca²⁺ was maintained at 2 mM during all the simulations.

Ca²⁺ Buffering Models

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

Single Pool

The exponential decaying Ca²⁺ pool was modeled as

$$ \frac{{d\left[ {{\text{Ca}}_i^{{2 + }}} \right]}}{{dt}} = - \frac{{{I_{{{\rm{C}}{{\rm{a}}^{{2 + }}}}}}(t)}}{{2Fd}} - \beta \left( {\left[ {{\hbox{Ca}}_i^{{2 + }}} \right] - \left[ {{\hbox{Ca}}_o^{{2 + }}} \right]} \right), $$

where $\left[ {{\hbox{Ca}}_i^{{2 + }}} \right]$ is intracellular Ca²⁺, $\left[ {{\hbox{Ca}}_o^{{2 + }}} \right]$ is Ca²⁺ at rest and is 45 nM, ${I_{{{\rm{C}}{{\rm{a}}^{{2 + }}}}}}(t)$ is total current through P-type and T-type Ca²⁺ channels, F is Faraday's constant, d is depth of a submembrane shell to define the volume for effective Ca²⁺ concentration, and β is the decay time constant. Instead of using commonly used values for β and d, we tuned these parameters to approximate the Ca²⁺ transients given by detailed Ca²⁺ dynamics, as described in the “Tuning of Single Pool Models” section.

Two Pools

Fast and slow Ca²⁺ transients, ${\left[ {{\hbox{Ca}}_i^{{2 + }}} \right]_{\rm{f}}}$ and ${\left[ {{\hbox{Ca}}_i^{{2 + }}} \right]_{\rm{s}}}$, were modeled using two pools with different decay time constants, β _f and β _s (β _f > β _s), and depths, d _f and d _s (d _f < d _s). The total intracellular Ca²⁺ $\left[ {{\hbox{Ca}}_i^{{2 + }}} \right]$ was modeled as the weighted sum of fast and slow Ca²⁺ transients.

$$ \left[ {Ca_i^{{2 + }}} \right] = {f_{\rm{f}}}{\left[ {Ca_i^{{2 + }}} \right]_{\rm{f}}} + {f_{\rm{s}}}{\left[ {Ca_i^{{2 + }}} \right]_{\rm{s}}} $$

The parameters β _f, β _s, d _f, d _s, f _f, and f _s were tuned to approximate the Ca²⁺ transients given by detailed Ca²⁺ dynamics, as described in the “Tuning of Double Pool” section.

Detailed Ca²⁺ Dynamics

The detailed Ca²⁺ dynamics model included calbindin (CB) and parvalbumin (PV) as buffers. In addition to Ca²⁺, both PV and 80% of CB were diffusible. A single surface-based Ca²⁺ pump was modeled using Michaelis–Menten kinetics [35]. The kinetics of CB and PV was obtained from a pre-existing Ca²⁺ dynamics model [25] and is described in Table 1. The kinetics of the pump was tuned to approximately match the Ca²⁺ decay measured at very high intracellular Ca²⁺ concentration [36], where the buffers are almost saturated and the decay is solely due to pump and diffusion. The outer radial shell had a depth (d) of 0.1 μm, whereas other radial shells had a depth of 0.2 μm each.

Table 1 Detailed Ca²⁺ dynamics model parameters

Detailed Ca²⁺ Dynamics with Diffusion Compensation

The detailed Ca²⁺ dynamics model was modified by removing the diffusion of Ca²⁺ and buffers. A diffusion compensating mechanism (DCM) was included to compensate for the reduced removal of Ca²⁺ from the submembrane region due to the lack of diffusion towards the center of the compartment. The DCM follows a standard buffering scheme:

$$ \left[ {{\hbox{C}}{{\hbox{a}}^{{2 + }}}} \right] + \left[ {\hbox{DCM}} \right]\mathop {{{k_{\rm{off}}}}}{\overset{{{k_{\rm{on}}}}}{\longleftrightarrow}}\left[ {{\hbox{C}}{{\hbox{a}}^{{2 + }}} - {\hbox{DCM}}} \right] $$

The total concentration of DCM and rate constants (k _on and k _off) were tuned to compensate for the diffusion of Ca²⁺ and mobile buffers in the absence of radial diffusion. Moreover, the depth of the single submembrane shell, which was used to compute the effective Ca²⁺ concentration for K_Ca channels, was also a tunable parameter.

Tuning of Ca²⁺ Buffering Models

Generating Target Traces

A single compartment was modeled with a diameter of 4 μm and a length of 20 μm. The P-type channel was included in the model for Ca²⁺ influx, and the detailed Ca²⁺ dynamics were used as a buffering model. A voltage step protocol, shown in Fig. 1a, was used to depolarize the compartment to the voltage at which physiological dendritic Ca²⁺ spikes are generated (an example of an experimental recording is shown for comparison in Fig. 1b). Then, different conductance values for P-type Ca²⁺ channel were used to generate peak amplitude of Ca²⁺ concentrations of 0.5, 1, 2, 4, and 8 μM in the volume defined by the submembrane shell. These simulated calcium transients formed the target traces for the automated fitting.

Fig. 1

a The depolarization command used to estimate Ca²⁺ profiles during generation of Ca²⁺ spikes. The peak voltage −22 mV was computed by taking an average of membrane potential during Ca²⁺ spikes (as shown in b), and the time of peak voltage ~12 ms was estimated from their average duration. b An example of electrophysiological dendritic Ca²⁺ spikes (provided by Ede Rancz and Michael Häusser, UCL, UK). Note that the small spikelets caused by retrograde conduction of somatic action potentials are not modeled in this study

Similarly, the experimental voltage trace of Ca²⁺ spikes, shown in Fig. 1b, was used as the voltage clamp control signal to generate calcium transients with peak amplitudes of Ca²⁺ concentrations of ~1.25, 2.5, 5, and 10 μM. These simulated calcium transients were later used to tune parameters for pool-based Ca²⁺ buffering models.

Tuning of Single Pool Models

The same model as described above was used, but now with a single pool. We used Neurofitter [37], an automated parameter search method to find optimal values of d and β to match the single pool simulations with the target traces. We used the random search mode, and model traces were compared to the target using a standard root mean square (RMS) error measure. While using the step voltage protocol (shown in Fig. 1a), since each of the traces used in fitting had two distinct features: (1) a fast rise and decay and (2) a slow decay (Fig. 3), we separated those features by specifying two time periods, 500 to 550 ms and 550 to 5,000 ms, over which separate RMS values were computed and added together. Similarly, while using the experimental voltage protocol (shown in Fig. 1b), we separated each of the traces used in fitting by specifying two time periods, 10 to 25 ms and 25 to 50 ms, over which separate RMS values were computed and added together. The best values for d and β obtained to fit the data from Fig. 1b were used with the single pool model to simulate dendritic Ca²⁺ spikes.

Tuning of Double Pool

Similar procedures were used to search for the values of β _f, β _s, d _f, d _s, f _f, and f _s values.

Tuning of DCM

We replaced diffusion of Ca²⁺, of free mobile buffers, and Ca²⁺-bound mobile buffers (keeping the total buffer concentrations constant) from the detailed Ca²⁺ dynamics model with a single equation for DCM. Similar procedures were used as in the case of tuning the single pool, except that we only used target traces generated using the step voltage protocol (shown in Fig. 1a), to find values of [DCM], k _on, k _off, and d. In order to find a robust solution, Neurofitter was run with five different seeds for random number generation. The five best solutions found using each seed were selected. We assume that the parameters for the diffusion compensating mechanism should be robust for different Ca²⁺ concentration peak amplitudes, but not for different compartment diameters. Therefore, the fitting procedure was repeated for diameters of 0.8, 2, 4, 6, 8, 10, 12, 16, and 20 μm.

DCM Parameter Interpolation

In actual compartmental models of neurons [38], a compartment can have an arbitrary diameter. Therefore, the parameters fitted to specific diameters were fitted to continuous functions. The 25 best solutions for each of the four parameters were plotted against corresponding diameter values. Matlab routines were used to fit the data with a single exponential function, double exponential function, and fifth order polynomial function based on an RMS measure. Later, those functions were used to predict the diffusion compensating parameter values for any compartment diameter.

Tuning Channel Densities for Dendritic Ca²⁺ Spikes

A dendritic Ca²⁺ spike, with three spikelets, was generated using the detailed Ca²⁺ dynamics model by hand tuning the channel density values. Later, Neurofitter was used to search proper channel densities for models with single pool, double pool, or DCM to generate similar dendritic Ca²⁺ spikes. The search was based on the RMS error measurement. We used five features for the fitness function. Three features comprised the three spikelets; other features were the pre-spike depolarization and post-spike hyperpolarization.

Results

Detailed Ca²⁺ Dynamics Model

A single compartment model with 4 μm diameter and 20 μm length was simulated with a P-type Ca²⁺ channel and the detailed Ca²⁺ dynamics model. The detailed Ca²⁺ dynamics model included calbindin (CB) and parvalbumin (PV) as buffers. The four binding sites of CB were simulated by assuming a pair of fast and a pair of slow binding sites as, respectively, a single fast and a single slow binding site (described as CB_f_s). PV was simulated together with Mg²⁺. Due to medium affinity of PV for Mg²⁺, its binding kinetics to Ca²⁺ is significantly affected by Mg²⁺. Diffusion of Ca²⁺, PV, and 80% of CB was included in the model. We simulated the voltage step of Fig. 1a in our model using different conductance values for P-type Ca²⁺ channel to generate peak amplitudes of 0.5 to 8 μM Ca²⁺ concentration. As the intracellular Ca²⁺ concentration started rising due to the influx through the P-type channel, the buffers and the pump became active. The changes in concentration of Ca²⁺ and Ca²⁺-bound buffers and outflux due to the pump over time are shown in Fig. 2.

Fig. 2

Intracellular Ca²⁺ profiles simulated using the depolarization command shown in Fig. 1a, for different amplitudes of Ca²⁺ influx (indicated by color). The panels show, respectively, Ca²⁺ concentration ([Ca²⁺]), Ca²⁺ bound to PV (PV.ca), Ca²⁺ bound to slow binding site of CB (CB.f.ca), Ca²⁺ bound to fast binding site of CB (CB.ca.s), Ca²⁺ bound to fast and slow binding sites (CB.ca.ca), Ca²⁺ bound to slow binding site of immobile CB (iCB.f.ca), Ca²⁺ bound to fast binding site of immobile CB (iCB.ca.s), Ca²⁺ bound to fast and slow binding sites (iCB.ca.ca), and pump current over a period of 2,000 ms

Dynamics of Single Pool and Double Pool Models

The pool parameters were tuned with Neurofitter [37], as described in “Materials and Methods” section. The best solution for the single pool, using the step voltage protocol (Fig. 1a), was β = 1.35/ms and d = 0.891 μm. The best parameter values for the double pool were β _f = 3.77/ms, d _f = 0.351 μm, β _s = 0.00306/ms, d _s = 0.928 μm, f _f = 0.994, and f _s = 0.006.

Figure 3 compares the intracellular Ca²⁺ profiles obtained using the tuned single pool and double pool models with those in the detailed Ca²⁺ dynamics model. The comparison demonstrates that the single pool approximates the Ca²⁺ transient only during the active Ca²⁺ influx and the initial fast decay. Moreover, the approximation is valid only for a small range of Ca²⁺ influx. Therefore, the single pool cannot approximate the effects of physiological buffering. The double pool can approximate Ca²⁺ transient during the active Ca²⁺ influx and following fast and slow Ca²⁺ decays, but its validity is also limited to a small range of Ca²⁺ influx values. Although the double pool model could not approximate the full range of detailed Ca²⁺ concentrations, it approximated the major components of the detailed model reasonably well, i.e., the fast and slow transients. Therefore, a double pool model is expected to provide a better control over activation of K_Ca channels than a single pool model.

Fig. 3

Comparison of Ca²⁺ profiles generated with a voltage step protocol using single pool, double pool (parameters specified in the text), and detailed dynamics model. Different Ca²⁺ concentration peak amplitudes of a 0.5, b 1, c 2, d 4, and e 8 μM were simulated to demonstrate the problems the single pool or double pool models have in capturing the complex dynamics of the detailed model. See text for parameters of the pool models

Since these estimated parameters belong to phenomenological models, the parameter values could vary based on the applied protocol. To estimate the ability of pool-based models to capture more complex Ca²⁺ kinetics, we applied also the experimental waveform (Fig. 1b) and estimated new model parameters. The best solution for the single pool was very different from the model used in Fig. 3: β = 6.86/ms and d = 0.169 μm. The same applied to the best parameter values for the double pool, though the ratio between fast and slow pool did not change: β _f = 7.33/ms, d _f = 0.167 μm, β _s = 0.00795/ms, d _s = 0.683 μm, f _f = 0.995, and f _s = 0.005. The behavior of these models (Supplementary Fig. 1) was similar to that observed in Fig. 3: their optimal concentration range is limited, and a single pool model cannot capture the slow transients.

Diffusion Compensation Model

Modeling buffers and pumps without diffusion of Ca²⁺ and buffers resulted in an accumulation of Ca²⁺ in the submembrane volume defined by a shell of depth 0.1 μm in compartments with a diameter of 1 μm or larger (Fig. 4). Therefore, we developed a compensating mechanism, DCM, to reduce the steep rise in peak amplitudes of intracellular Ca²⁺. Because the detailed Ca²⁺ dynamics model had multiple buffers with different buffer affinities, the effects of diffusion depended on the amount of Ca²⁺ entry. We developed DCM to robustly reproduce model Ca²⁺ transients in compartments with a wide range of diameters and for physiological peak amplitudes in the submembrane shell of between 0.5 and 8 μM.

Fig. 4

Comparison of Ca²⁺ profiles of different peak amplitudes, a 0.5, b 1, c 2, d 4, and e 8 μM, simulated in a compartment with a diameter of 2 μm, using the detailed Ca²⁺ dynamics model, detailed Ca²⁺ dynamics model without diffusion, and DCM model. The removal of diffusion from the detailed model resulted in a steep rise in Ca²⁺, which could be well compensated by DCM. The inset (f) below the panel (c) highlights Ca²⁺ profiles simulated using five different sets of optimal DCM parameters found by using Neurofitter

The parameters for DCM were estimated using Neurofitter [37]. Examples for different levels of Ca²⁺ influx using the voltage step protocol demonstrate the effect of removing diffusion, and the good compensation by the DCM over the entire range is shown in Fig. 4.

Derivation of the Parameter Predictors

We showed in Fig. 4 that DCM can effectively compensate for excluding diffusion, but the parameter fitting produced values only for specific compartment diameters. To use DCM in a morphologically detailed PC model, we need a way to predict the DCM parameter values for any diameter present in the PC. Figure 5 shows the results of automated parameter fitting (see “Materials and Methods” section), for nine diameters. Using Matlab, the four parameters (concentration of DCM, [DCM]; forward rate constant, k _on; the backward rate constant, k _off; and the depth of submembrane shell, d) were fitted to functions (Fig. 5). The prediction functions for the parameters are

$$ \begin{array}{*{20}{c}} {\left[ {\hbox{DCM}} \right] = 64.2 - 57.3{e^{{ - \frac{\rm{diam}}{{1.4}}}}}} \hfill \\{{k_{\rm{on}}} = 0.162 - 0.106{e^{{ - \frac{\rm{diam}}{{2.29}}}}}} \hfill \\{{k_{\rm{off}}} = \left\{ {\begin{array}{*{20}{c}} {0.000267 + 0.0167{e^{{ - \frac{\rm{diam}}{{0.722}}}}} + 0.0028{e^{{ - \frac{\rm{diam}}{4}}}}} \hfill & {i{\hbox{f}}\,{\hbox{diam}} \geqslant {2,}} \hfill \\{0.003} \hfill & {\hbox{otherwise}} \hfill \\\end{array} } \right\}} \hfill \\{d = \frac{\text{diam}}{{4 \times \left( { - 0.674 + 1.94\left( {\hbox{diam}} \right) + 0.289{{\left( {\hbox{diam}} \right)}^2} - 3.33 \times {{10}^{{ - 2}}}{{\left( {\hbox{diam}} \right)}^3} + 1.55 \times {{10}^{{ - 3}}}{{\left( {\hbox{diam}} \right)}^4} - 2.55 \times {{10}^{{ - 5}}}{{\left( {\hbox{diam}} \right)}^5}} \right)}}} \hfill \\\end{array} $$

Fig. 5

Distribution of estimated DCM parameters ([DCM], k _on, k _off, and d; mean ± STD, n = 25: the five best solutions from five runs with different seeds for random number generation) for diffusion compensation plotted against compartment diameter. These were fitted respectively by a exponential function, b exponential function, c double exponential function, and d fifth order polynomial function

The fitted functions are very close to the mean values of the parameter distributions and reasonably capture the relationship between the DCM parameters and the diameter of the compartment. Therefore, we expect these functions to predict DCM parameters for any diameter in the PC. This is demonstrated for two unfitted diameters in Fig. 6. It shows that the fitted functions indeed produce good compensation for diffusion in the detailed model. Therefore, we can use the fitted functions with morphologically detailed PC models.

Fig. 6

Comparison of Ca²⁺ profiles of different peak amplitudes simulated using detailed Ca²⁺ dynamics (blue) and diffusion compensated Ca²⁺ dynamics (green) with predicted values of parameters for diffusion compensation mechanism from the functions in Fig. 5. a Diameter 4.8 μm. b Diameter 14 μm

Dendritic Ca²⁺ Spikes

Dendritic Ca²⁺ spikes with multiple spikelets comparable to those observed in vitro (Fig. 1b) were simulated using a single pool, a double pool, detailed Ca²⁺ dynamics, and DCM. Since the buffering models used in the simulations were different approximations of detailed Ca²⁺ dynamics, we expected the intracellular Ca²⁺ transients generated by these models to show considerable differences (Supplementary Fig. 3). Therefore, we first optimized values for the maximal conductance (G _max and P _max) of ion channels (current activation profiles are shown in Supplementary Fig. 4) in the model to generate as comparable Ca²⁺ spikes as possible. This is identical to the procedure that would be used to generate a neuron model with only one type of Ca²⁺ dynamics model [38]. The resulting values are listed in Table 2. The resultant voltage traces of the Ca²⁺ spikes are shown in Fig. 7 (see also Supplementary Fig. 2).

Table 2 Maximal conductance values used to generate dendritic Ca²⁺ spikes

Fig. 7

Dendritic Ca²⁺ spikes generated using different Ca²⁺ buffering models (aligned at the peak of the first spikelet in a, b, and f). The parameters of pool-based models are those used in Supplementary Fig. 1. The conductance values used to generate these spikes are listed in Table 2. a First burst of Ca²⁺ spikes. b Burst of Ca²⁺ spikes around 57 s. c, d Spontaneous Ca²⁺ spike bursting over 100 s: black asterisk indicates the first burst of Ca²⁺ spikes (shown in a), red asterisk indicates the burst of Ca²⁺ spikes around 57 s (shown in b). e Inter-burst interval (IBI) as a function of current injection. f Burst of Ca²⁺ spikes around 57 s with injection of 0.004 pA current

Comparing the dendritic Ca²⁺ spikes generated using different buffering models (Fig. 7a), we clearly see that the dendritic spikes generated using the pool-based models are different from the spikes generated using the other buffering models and strongly depend on the pool models used (compare Supplementary Fig. 2 with Fig. 7a). The amplitude as well as the width of second and third spikelets are different from the spikelets generated using the detailed Ca²⁺ model. Another noticeable difference exists in the after hyperpolarization of these bursts. These differences are due to the different activation of K_Ca channels, in particular of the SK channel (Supplementary Fig. 4). On the other hand, the dendritic Ca²⁺ spikes simulated using DCM approximated the Ca²⁺ spikes generated using the detailed Ca²⁺ dynamics model quite well.

A comparison of dendritic Ca²⁺ spikes over a time period of 50 ms may not reveal the effect of diffusion at longer time scales. Therefore, we ran all simulations for 200 s and compared the ability of different buffering mechanisms to capture the bursting behavior over that time. Fig. 7c, d shows the spontaneous dendritic Ca²⁺ bursting using the detailed Ca²⁺ dynamics model, DCM, single pool, and double pool models. Each of the bursts contains three Ca²⁺ spikelets (examples shown in Fig. 7a, b). The frequency of bursts in case of pool-based models is approximately five to six times higher than the frequency of bursts in the detailed Ca²⁺ dynamics model. However, the frequency of bursts for the DCM model is similar to that of the detailed model.

Comparing the burst of Ca²⁺ spikes around 57 s (marked with red asterisk in Fig. 7c, d) in Fig. 7b, we see that the shape has changed compared to the first burst (Fig. 7a), which is due to the buildup of Ca²⁺ over longer time. Pool-based models failed to capture this burst shape adaptation, while DCM reproduced it fairly well.

Next, we investigated the ability of DCM to approximate characteristics of diffusion over the longer time period when the bursting period was modified by current injection. The mean inter-burst interval plotted against the injected current is shown in Fig. 7e. The behavior of detailed Ca²⁺ dynamics model and DCM is similar at most current injection levels, except for bumps around 0.002 pA. These bumps represent the transition of the firing pattern from bursts with three spikelets to bursts with two spikelets (shown in Fig. 7f), which occurs at lower current levels in the detailed model. As the pool models do not show burst shape adaptation, they also entirely fail to capture this transition to two spikelet bursts. The spontaneous bursting period range is limited. Therefore, to compare the different buffering models at faster bursting rates, we injected current pulses at 1 Hz (Supplementary Fig. 5). Under these conditions, the Ca²⁺ buildup is more pronounced leading to a SK current mediated after hyperpolarization and a progressive decrease in the number of spikes in each burst, but the relative accuracy of the different buffering models is similar to that shown in Fig. 7.

Based on the data shown in Fig. 7, we conclude that DCM is quite effective over a wide range of conditions. While the approximation is not perfect, it is clearly much superior compared to the pool-based models. Finally, we briefly explored the robustness of DCM to changes in the Ca²⁺ buffers. Supplementary Fig. 6 shows a comparison of detailed and DCM models when the concentration of the CB and PV buffers is either halved or doubled (DCM parameters unchanged but ion channel densities retuned). The DCM approximation is very good for the half buffer concentration model, but less adequate for the double buffer concentration one. Therefore, we recommend to retune the DCM parameters whenever the detailed Ca²⁺ dynamics model is changed.

Discussion

We compared the effectiveness of different Ca²⁺ buffering mechanisms in controlling K_Ca channels expressed in PCs. The PC is known to express high amount of mobile endogenous buffers like calbindin and parvalbumin [39]. Ca²⁺ entering through voltage-gated channels is quickly removed by these buffers [19]. The internal Ca²⁺ stores and pumps together with calbindin and parvalbumin control the intracellular Ca²⁺ available to K_Ca channels. Therefore, careful modeling of intracellular Ca²⁺ dynamics is essential to simulate the physiological control of K_Ca channels.

The detailed Ca²⁺ dynamics model introduced by Schmidt et al. [25] includes calbindin, parvalbumin, pumps, and diffusion of Ca²⁺, parvalbumin, and calbindin. The model was tuned for Ca²⁺ dynamics in a single spine coupled to a dendrite. Since the spine is a small compartment, only longitudinal diffusion was modeled from the spine through the spine neck. In this study, we modeled dendrites, which are relatively large compartments; therefore, we modified the Schmidt et al. [25] model to incorporate radial diffusion towards the center of the compartment using the formalization implemented in NEURON [40]. Recent Ca²⁺ dynamics models [25, 41] used different pump rates to fit their data. We tuned the pump rate parameters to match the Ca²⁺ profiles recorded from a PC at high concentrations of Ca²⁺, at which most of the buffers were saturated, and the decay could be attributed to diffusion and pump only [36].

The single Ca²⁺ pool has been the choice for modeling intracellular Ca²⁺ in almost all compartmental models. The volume of sub-compartment to compute the effective Ca²⁺ concentration for activation of K_Ca channels is often defined by a submembrane shell with a depth of 0.2 μm [20]. However, a wide range of decay rates (β values) is used from 10/ms [20] to 0.02/ms [42]. These parameter values used in single pool models can reasonably approximate the changes in microdomains [23] and can be used to control a BK-type K_Ca channel [17, 18]. SK-type K_Ca channels, which are relatively far from Ca²⁺ channels [18], require smaller amounts of Ca²⁺ (0–2 μM) and activate at relatively long time scales. Therefore, it is impossible to control both BK and SK channels correctly using a single pool. We confirmed this prediction by comparing the Ca²⁺ profiles (Supplemental Fig. 3) and dendritic Ca²⁺ spikes (Fig. 7 and Supplementary Figs. 2 and 5) generated by single pool and detailed models. We also considered a Ca²⁺ pool model with two time constants as a phenomenological solution to the problem. While such models fitted simple Ca²⁺ spikes better (Fig. 3), their concentration range is limited, and they fail to simulate the period of repetitive bursting and burst shape adaptation (Fig. 7).

The limitation of using detailed Ca²⁺ dynamics is the computational cost of underlying diffusion [27]. We computed run times for simulations of 10 s with detailed Ca²⁺ dynamics and different compartment diameters (Fig. 8). It took 45 s to run the simulation for a 2-μm diameter compartment. The run time increase was supra-linear with increase in diameter; for 20 μm, it was 1,076 s. This suggests that the computation cost will increase tremendously for simulating a morphologically detailed PC model, which we want to avoid. In this study, we proposed a method to approximate the effect of diffusion. Since only Ca²⁺ in the submembrane region is required to control K_Ca channels, it is not essential to compute the Ca²⁺ concentration in and around the center of compartments. We replaced diffusion with a compensating mechanism, DCM, that made up for the Ca²⁺ diffusing towards the center. We confirmed that the DCM approximated the effect of diffusion robustly in the range of 0.5–8 μM Ca²⁺ for dendritic compartments with diameters of 1–20 μm and over time ranges from tens of milliseconds to hundreds of seconds. Though BK channels are reported to sense approximately 100 μM of Ca²⁺, those high concentrations are only available in nanodomains [18]. For volumes defined by ~0.1 μm shells, 8–10 μM of Ca²⁺ concentration should represent the limit of the physiological range [17]. To achieve this level of accuracy, it is important to fit the DCM parameters to the detailed Ca²⁺ dynamics model used (Figs. 4 and 5). For dendritic compartments with diameters less than 1 μm, the effect of diffusion is not significant and can be compensated by tuning only the depth (d) of the submembrane shell (results not shown in this work).

Fig. 8

Comparison of run time for NEURON simulations of 10 s, with detailed Ca²⁺ dynamics and with DCM using compartments with different diameters. All the simulations were run on an Apple MacBook Pro, Intel Core 2 Duo 2.33 GHz

The method proposed in this paper is not limited to modeling PCs. The DCM parameters can be tuned to replace diffusion in large compartmental models of other neuronal types. To formulate robust DCM parameters, it is necessary to have a faithful detailed Ca²⁺ dynamics model of the specific neuronal type. Then, for the desired physiological range of Ca²⁺ influx and range of diameters in the morphology, DCM parameter equations can be derived to replace diffusion.

The use of the DCM resulted in a decrease of the run time to 8 s with a loss of dependence on diameter (Fig. 8). We conclude that using a detailed Ca²⁺ buffering model combined with DCM is a clear improvement compared to the Ca²⁺ pool for modeling Ca²⁺ dynamics in a large neuron model.