The Slow Dynamics of Intracellular Sodium Concentration Increase the Time Window of Neuronal Integration: A Simulation Study

Introduction

In modeling the electrical behavior of neurons, the intracellular and extracellular concentrations of the major ions (except Ca²⁺) are typically taken as constant, due to two assumptions: (a) Existence of effective regulating mechanisms that maintain homeostatic conditions; (b) The total number of ions that flow across the membrane during electrical activity is rather small, either because the current is weak or lasts for a short time. Ionic concentration changes have been taken into account in sporadic cases such as [Cl⁻]_i fluctuations, which may result in the reversal of GABA_A receptor-mediated currents (Wagner et al., 1997; Ben-Ari, 2002), or [K⁺]_o, which can increase as a result of intense population activity, leading to initiation of seizures (Moody et al., 1974; Yarom and Spira, 1982; Cressman et al., 2009).

In the few studies where changes in neuronal [Na⁺]_i were examined, it was used to detect and analyze Na⁺ conductance hot-spots (Fleidervish et al., 2010) or as part of point neuron modeling in the context of seizures or bursting (Cressman et al., 2009; Barreto and Cressman, 2011). Forrest et al. (2012) provided a simplified model of changes in [Na⁺]_i in a cerebellar Purkinje cell in order to examine the effect of Na⁺-K⁺ pump current. The study of [Na⁺]_i changes and their consequences is significantly more common in the context of cardiac myocyte modeling, where Na⁺ is treated as an intracellular messenger (Swietach et al., 2015).

[Na⁺]_i dynamics are governed by Na⁺ influx, pumping, and longitudinal diffusion. In the present study, the primary sources of Na⁺ influx are voltage-dependent Na⁺ conductance (either transient or persistent) and glutamatergic synaptic conductance. Contrary to Ca²⁺, longitudinal diffusion is an important factor in [Na⁺]_i dynamics, since Na⁺ is free to diffuse in the cell without buffering. Moreover, at short time scales, Na⁺ diffusion seems to be the major contributing factor to changes in [Na⁺]_i (Fleidervish et al., 2010). Due to the prominent role of longitudinal diffusion, [Na⁺]_i dynamics are highly dependent on the geometry of the compartment in question and the availability of diffusion sinks in the adjacent compartments. In some cases (e.g., thin processes and/or lack of diffusion sinks), [Na⁺]_i dynamics may be several orders of magnitude slower than the dynamics of membrane potential or [Ca²⁺]_i, thus allowing prolonged temporal integration, or “memory,” of past inputs or activity (Forrest et al., 2012; Zylbertal et al., 2015). Conversely, the relative rapid diffusion of Na⁺ entails that its concentrations in the sub-membrane area and the core of the compartment are in equilibrium.

The main Na⁺ extrusion mechanism is the electrogenic Na⁺-K⁺ pump, which transports two K⁺ ions into the cell for every three Na⁺ ions it transports out of the cell, thus producing net outward current using energy from ATP hydrolysis (Carafoli, 1991). For every unit of charge introduced into the cell as Na⁺ inward current, a net 1/3 unit of charge will contribute to an outward current carried by the Na⁺-K⁺ pump. This outward current exerts a variable effect on the membrane potential and thus on neuronal information processing (Forrest, 2014). Another source for Na⁺-dependent outward current may be Slick and Slack K⁺ channels (Kaczmarek, 2013), which will not be discussed here.

Another important consequence of changes in [Na⁺]_i are variations in Na⁺ Nernst potential (E_Na). Similarly to [Na⁺]_i, this value is usually regarded as constant, although as we demonstrate below, bouts of spiking activity or synaptic inputs are sufficient to shift it enough to cause significant spike amplitude adaptation and apparent synaptic depression.

Lastly, [Na⁺]_i affects Ca²⁺ dynamics, since the Na⁺-Ca²⁺ exchanger, an important Ca²⁺ extrusion mechanism, is powered by the trans-membrane Na⁺ concentration gradient (Blaustein and Lederer, 1999). Changes in this gradient reduce or even reverse the activity of Ca²⁺ extrusion through the exchanger (Schäfer, 2001; Scheuss et al., 2006; Zylbertal et al., 2015). Since Ca²⁺ dynamics are generally faster than Na⁺ dynamics, [Na⁺]_i dictates the stable-state [Ca²⁺]_i, or the level to which [Ca²⁺]_i decays. As we demonstrate below, this may have surprising consequences, such as tagging of active synapses by prolonged elevation in [Ca²⁺]_i.

Here, we present a detailed analysis of [Na⁺]_i dynamics and their consequences in three detailed conductance-based models representing three neuronal types: an accessory olfactory bulb (AOB) mitral cell, a cortical layer V pyramidal cell, and a cerebellar Purkinje cell. We do not take into account changes in [K⁺]_o, which are usually attributed to population activity rather than the behavior of a single cell. Such changes may enhance or reduce the effects we describe, depending on the context.

This modeling was inspired by our experimental and computational analysis of persistent activity and infra-slow bursting in AOB mitral cells (Zylbertal et al., 2015), where we showed that [Na⁺]_i dynamics have a critical role in generating both phenomena. Here, we further examine the role of [Na⁺]_i dynamics in the AOB mitral cell model, and extend it to other neuronal types using a revised version of previously published pyramidal cell and Purkinje cell models (Hay et al., 2011; Masoli et al., 2015).

These models were chosen because they combine detailed morphology and biophysics. Contrary to the AOB mitral cell model, their geometry was not reduced and represents the original reconstruction of cell morphology. In addition to representing two different neuronal types, these models also diverge in several design choices, stemming from the different goals they were meant to fulfill (see section Methods). For example, the Purkinje cell model contains a detailed axon, Na⁺ channel mechanism, and Ca²⁺ dynamics, while the pyramidal cell model has a realistic and detailed distribution of conductances along its dendritic tree. These differences allowed us to explore a broad range of possibilities for [Na⁺]_i dynamics and their consequences, without adhering to any specific design. More broadly, it should be noted that our goal is not to provide accurate predictions for the behavior of each neuronal type we examine, but rather to explore the landscape of possibilities arising from modeling [Na⁺]_i dynamics in neurons.

Our results demonstrate three primary pathways by which [Na⁺]_i dynamics may influence neuronal computation: (1) prolonged modulation of local EPSP amplitude due to reduction of E_Na; (2) prolonged activity-dependent hyperpolarization due to increased activity of the Na⁺-K⁺ pump; (3) synaptic tagging by prolonged Ca²⁺ elevation following Na⁺ loading through repeated inputs.

Methods

All models were implemented in the NEURON simulation environment with Python (Hines and Carnevale, 1997; Hines et al., 2009). All current injections and voltage recordings were done at the soma.

AOB Mitral Cell Model

For a full description of this model, see Zylbertal et al. (2015, Supporting Information). Briefly, the model was based on experimental measurements and morphological reconstruction of an AOB mitral cell (Myatt et al., 2012; Zylbertal et al., 2015), and included influx, diffusion, and extrusion of Na⁺ and Ca²⁺ (the latter was restricted to the dendritic tuft). Sodium concentration change in each compartment by influx and efflux through the membrane is given by:

Where $I_{\text{N}{\text{a}}^{+}}$ is the net Na⁺ current (the product of the local current density and the membrane area), F is Faraday constant, and V is the compartment volume (fully accessible to Na⁺). Na⁺ is free to diffuse longitudinally, and is pumped out of the cell by the Na⁺-K⁺ pump, modeled using a simple kinetic scheme (see below). Ca²⁺ accumulation was modeled similarly: the entire volume of each compartment was accessible to Ca²⁺, yet, it was not free to diffuse. Ca²⁺ buffering and pumping was also modeled using simple akinetic schemes, while the Na⁺-Ca²⁺ exchanger current followed:

(Courtemanche et al., 1998)

Where I_NCX(max) is the maximal current, γ is the voltage dependence parameter, k_m(Na) and k_m(Ca) are the Na⁺ and Ca²⁺ dissociation constants, k_sat is the saturation factor, and R and F are the gas constant and Faraday constant, respectively.

The model assumes the presence of active Na⁺ channels in the apical dendrites and tufts (Ma and Lowe, 2004), as well as non-uniform channel properties across different compartments (Colbert and Pan, 2002). Evolutionary multi objective optimization algorithm (EMOO, Deb, 2001; Bahl et al., 2012) was used first to find a simplified (“lumped”) geometry that would reproduce the passive electrical properties of the detailed geometry. A second EMOO step was used to fit the model parameters, based on recorded electrophysiological and imaging data. Some membrane mechanisms were based upon published models hosted by ModelDB (Mainen and Sejnowski, 1996; Courtemanche et al., 1998; Lazarewicz et al., 2002; Hines et al., 2004; Korngreen et al., 2005). The model code is available online at: https://senselab.med.yale.edu/ModelDB/ShowModel.cshtml?model=185332.

Layer V Pyramidal Cell Model

This is an adaptation of a detailed model developed by Hay et al. (2011). Among the models presented in this paper, we selected the one that included an axon. The original model is based on reconstructed cortical layer V pyramidal cells (Le Bé et al., 2007), and was fitted using the EMOO algorithm based on experimental results derived from step current injection (Le Bé et al., 2007), Ca²⁺ spike statistics (Larkum et al., 1999), and back-propagating action potential properties (Larkum et al., 2001).

We modified the model by first introducing Na⁺ accumulation, diffusion, and a pumping mechanism to all compartments, similarly to the AOB mitral cell model (see above). We used the value of 0.3 μm²/ms for the Na⁺ diffusion coefficient, approximately the one measured experimentally in dendrites (Mondragão et al., 2016). In order to account for the apparent effect of dendritic spines, we changed this value to 0.03 μm²/ms in some simulations (see Supplementary Information). The Na⁺-K⁺ pump was modeled as in the AOB mitral cell, and was distributed using a similar order of magnitude (in mol/cm²: soma - 1·10⁻¹¹; axon - 5·10⁻¹²; dendrites - 1·10⁻¹⁵).

We next replaced the Ca²⁺ dynamics of the original model (simple exponential decay) with buffering, pumping, and Na⁺-Ca²⁺ exchange, as in the AOB mitral cell model. The electrogenic effect of the Na⁺-Ca²⁺ exchanger was removed, since its effect on the membrane potential is already taken into account in the fitting of the original model. We used the Ca²⁺ pump density value from the mitral cell model in the pyramidal cell dendrites, and a density six times higher in its soma and axon. We used the maximal current value of the Na⁺-Ca²⁺ exchanger from the mitral cell model throughout the pyramidal cell, and preserved its parameters. Additionally, we updated the Ca²⁺ channel models, so that the Goldman-Hodgkin-Katz equation, rather than the Nernst equation, is used to infer the Ca²⁺ electromotive force. The spatial resolution of the model was tripled to improve simulation accuracy of diffusional elements. This revised model is available online at: https://senselab.med.yale.edu/ModelDB/showModel.cshtml?model=230326.

Cerebellar Purkinje Cell Model

This is an adaptation of a detailed model developed by Masoli et al. (2015). The original model is based on a reconstructed dendritic tree and soma of a Guinea pig Purkinje cell (Rapp et al., 1994), to which a detailed description of the axon was later added, including the AIS and three nodes of Ranvier. It incorporates available knowledge on localization and gating of 15 Purkinje cell ionic channels, modeled using either Hodgkin-Huxley or Markovian formalism, as well as detailed Ca²⁺ dynamics with pumping, buffering, and approximation of radial diffusion. Free model parameters (maximal conductance values) were initially set according to estimates taken from the literature (Swensen and Bean, 2003), and later manually fine-tuned to match simulation results with somatically-recorded voltage traces elicited by step-current injections.

Since this model already included detailed Ca²⁺ dynamics, its adaptation was straightforward: We adapted it by introducing Na⁺ accumulation and pumping (as in the pyramidal cell model, using identical pump spatial distribution), and replacing some of the Ca²⁺ pumps with Na⁺-Ca²⁺ exchangers, which were not part of the original model. This was done by setting the maximal Na⁺-Ca²⁺ exchanger current in each compartment according to the compartment pump density (multiplied by 10⁸ mA/mol), and subsequently reducing the pump density by a factor of 10⁵.

Calculation of Na⁺ Entry Ratio

We calculated the Na⁺ entry ratio (Carter et al., 2009) in each model by dividing the integral of the global Na⁺ current during an action potential by the minimal charge needed to account for the action potential amplitude (i.e., the product of that amplitude and the membrane capacitance). Effective membrane capacitance was assessed by clamping the model cells soma to two different holding potentials, and measuring the capacitance transient. In the Purkinje cell model, this was done after removing the dendrites, due to the high load they exert on the soma, which results in overestimation of the apparent capacitance.

Simulation of Synaptic Currents

Local synaptic conductance was simulated as a sum of onset and offset exponentials, whose time constants represent a compromise between fast AMPA receptor-mediated current kinetics and slow NMDA receptor-mediated current kinetics (Spruston et al., 1995, τ_onset = 0.5 ms, τ_offset = 10 ms, g_max = 5 nS). The Na⁺ and K⁺ components of the conductance were inferred from its reversal potential (E_rev = 10 mV) and the resting ionic concentrations. Climbing fiber input was simulated as a concurrent input to all dendritic branches whose diameter is >2.2 μm, with τ_onset = 0.5 ms, τ_offset = 5 ms, Σg_max = 250 nS (Wadiche et al., 2001).

Realistic Distributed Input

The pyramidal cell model was used to simulate the effects of stochastic distributed inputs (Hay and Segev, 2015). Thousand pre-synaptic “cortical” spike trains followed Poisson distribution with an average firing rate of 1 Hz. Each of these pre-synaptic sources made five contacts in the pyramidal cell dendrites (for a total of 5,000 synaptic contacts), their location randomly chosen from a uniform distribution. Additional 1,000 pre-synaptic “thalamic” spike trains (having 1 Hz background firing rate and three 2–3 s long episodes of 5 Hz firing rate, representing sensory stimuli) made five contacts at sites chosen from a gamma distribution (k = 2.0, θ = 60 μm), resulting in higher density of contacts in proximal regions (Rah et al., 2013). These five contacts were either clustered (all in the same site), or each chosen individually. Excitatory synaptic conductance was modeled as detailed above, with reduced peak conductance (g_max = 0.5 nS) and slightly longer decay time (τ_offset = 15 ms) to account for NMDA receptor connections.

Since the original Na⁺-K⁺ pump model compensates for Na⁺ leak and background synaptic activity, it had to be adjusted to prevent gradual increase in [Na⁺]_i. This was done by lowering the [Na⁺]_i steady-state value of the pump to a quarter of its original value, and increasing its dendritic density to 50 ·10⁻¹⁵ mol/cm².

Inferring Na⁺-K⁺ Pump Current as a Function of [Na⁺]_i

The Na⁺-K⁺ pump is simulated according to the following reactions:

Therefore:

Where F is the Faraday constant. Assuming that [pumpNa] is in equilibrium ($\frac{d\left[pump\text{Na}\right]}{dt}=0$), [pump] and [pumpNa] can be inferred and substituted into (1), yielding:

Inferring Stable-State [Ca²⁺]_i as a Function of [Na⁺]_i

The plasma membrane Ca²⁺ pump is simulated like the Na⁺-K⁺ pump (with different stoichiometry). Summing the rate of the Ca²⁺ pump with the rate of the Na⁺-Ca²⁺ exchanger gives the total rate of Ca²⁺ change:

Solving for [Ca²⁺]_i when $\frac{d{\left[\text{C}{\text{a}}^{2+}\right]}_i}{dt}=0$ (equilibrium) yields a complex expression for stable-state [Ca²⁺]_i. SymPy Python package (Meurer et al., 2017) was used to generate and evaluate this expression for different values of [Na⁺]_i.

Results

Simulated Na⁺ Dynamics in Multi-Compartmental Model Cell

To examine the Na⁺ dynamics in the three model cells, we stimulated each of them with a 30 Hz, 2-s long pulse train. The pulse amplitude and duration were adjusted in each cell, so that each pulse would evoke a single action potential, resulting in a 30 Hz firing frequency. In the Purkinje cell model, hyperpolarizing DC injection was used to prevent spontaneous firing.

First, in order to test how realistic the models are in terms of Na⁺ intrusion dynamics, we quantified the sodium entry ratio (SER), i.e., the ratio between overall Na⁺ entry in each spike and the total charge needed to account for the spike amplitude (see Methods, Carter et al., 2009). The cortical pyramidal cells exhibited the most efficient dynamics (SER = 1.76), followed by the AOB mitral cell (SER = 3.6) and the Purkinje cell (SER = 5.6). These values are similar in magnitude and relations to experimentally measured values (Carter et al., 2009), and are slightly higher, due to the underestimation of the Na⁺ current in experiments compared to models.

In each of the model cells (Figure 1A), the axon initial segment (AIS) bears (by definition) the highest density of voltage-gated Na⁺ channels, therefore, the largest build-up in [Na⁺]_i occurs in this compartment (Figure 1B). The extent and duration of this build-up is determined primarily by the magnitude of Na⁺ influx in each spike and the efficiency of Na⁺ diffusion out of the AIS. Compared to these processes, Na⁺ active extrusion (pumping) is too slow to have a significant effect on Na⁺ dynamics in the AIS (Fleidervish et al., 2010). Action potentials in the mitral cell model have the largest duration (Figure 1B, inset, blue), followed by the spikes generated by the pyramidal cell (green) and the Purkinje cell (red). Furthermore, the diffusion of Na⁺ out of the AIS in the mitral cell model is the slowest among the three. As a consequence, AIS [Na⁺]_i in the mitral cell model is increased during the spike train by 40 mM, as compared to the pyramidal cell (18 mM) and Purkinje cell (6 mM) models. Moreover, contrary to the pyramidal and Purkinje cell models, in the mitral cell model Na⁺ concentration does not reach equilibrium during the spike train, and thus longer trains would result in even larger [Na⁺]_i build-up.

The large volume of the soma prevents extreme changes in somatic [Na⁺]_i (Figure 1C, note the scale). The lack of good diffusion sinks, however, results in a slow return to baseline [Na⁺]_i, governed mostly by active Na⁺ extrusion. Among the models tested, the pyramidal cell model has the highest density of voltage-gated Na⁺ channels in its soma, leading to the largest somatic [Na⁺]_i elevation (6 mM) of all. The only major somatic diffusion sink is the proximal dendrite, where [Na⁺]_i follows a similar trajectory as the soma (Figure 1D).

Similarly to the soma, no strong diffusion sinks are available near the distal dendrite (Figure 1E), hence Na⁺ concentration dynamics are determined by the density of active Na⁺ conductance and extrusion mechanisms. In the pyramidal and Purkinje cell models the level of dendritic voltage-gated Na⁺ conductance is very low, resulting in nearly no [Na⁺]_i build-up in the distal dendrite. Conversely, the mitral cell model has a high level of dendritic voltage-gated Na⁺ conductance (Ma and Lowe, 2004), causing a modest (4 mM) elevation in dendritic [Na⁺]_i. Nevertheless, this elevation is extremely prolonged, since Na⁺ extrusion from the distal dendrites is governed solely by pumping, as the lack of sinks prevents effective diffusion. Such prolonged dendritic elevation in [Na⁺]_i has been indeed observed experimentally (Fleidervish et al., 2010; Zylbertal et al., 2015) and can potentially act as a non-leaky integrator of neuronal activity (see below).

The Purkinje cell model consists of a detailed axon with myelinated sections and nodes of Ranvier, enabling the simulation of [Na⁺]_i build-up along the axon. Due to the strong diffusion sinks created by the nearby myelinated sections and the limited Na⁺ conductance, [Na⁺]_i build-up in the third node of Ranvier (similarly to the other nodes) is relatively small and transient (Figure 1F).

Overall, a large variety of temporal profiles of [Na⁺]_i were observed in the different cell models as well as at different compartments of the same model. The differences among the temporal profiles reflect differences in Na⁺ conductance, Na⁺ extrusion mechanisms, and geometrical relations between neighboring compartments.

Synaptic Activity Results in Prolonged Local Elevation in [Na⁺]_i

Since the pyramidal and Purkinje cell models lack significant dendritic voltage-dependent Na⁺ conductance, firing activity does not alter their dendritic [Na⁺]_i (Figure 1E, red and green lines). We therefore examined the effect of another possible source of intracellular Na⁺, the synaptic current, on dendritic [Na⁺]_i. Since the local dendritic properties in the context of Na⁺ are similar in the two models, we used the pyramidal cell model (Figure 2A) to demonstrate this effect.

A 5 Hz, 4-s long train of synaptic inputs, followed by two isolated subsequent inputs given to examine prolonged influences (Figure 2B, triangles) were simulated at a distal site in the apical dendritic tree (Figures 2A,B, arrow). The Na⁺ component of the synaptic conductance was inferred from the synaptic reversal potential (+10 mV), assuming that only Na⁺ and K⁺ conductances contribute to the total synaptic current. As apparent by both the temporal (Figure 2B) and the spatial (Figure 2C) profiles of the dendritic [Na⁺]_i, the synaptic activation caused a moderate (8 mM) increase in [Na⁺]_i at the synaptic site, which decayed by diffusion to neighboring segments (color coded by proximity as in A) in a few seconds.

The presence of dendritic spines results in non-linear (“anomalous”) dendritic diffusion (Santamaria et al., 2011). To linearly approximate the effect of dendritic spines, which drastically reduce Na⁺ spread (resulting in a smaller apparent Na⁺ diffusion coefficient), we repeated the simulation with a 10-fold reduction in the dendritic Na⁺ diffusion coefficient. Compared to a spiny dendrite, this reduction results in a conservative estimation of [Na⁺]_i in the stimulated site and a good approximation of [Na⁺]_i in neighboring sites, as demonstrated by simulating the two cases (Supplementary Figure). Indeed, under these conditions (Figures 2D,E), the change in [Na⁺]_i is much larger (25 mM), more localized, and persists for a longer period of time.

Compared with the changes in dendritic [Na⁺]_i, alterations in membrane potential due to synaptic activity (Figures 2F,G) are considerably faster and have a larger spatial extent. Synapse-driven [Na⁺]_i changes are therefore more compartmentalized than electrical changes and persist for longer periods of time.

Elevated [Na⁺]_i Causes a Reduction in Na⁺ Nernst Potential, Affecting the Amplitude of Action Potentials and EPSPs

One potential outcome of an increase in [Na⁺]_i is a reduction in the Na⁺ Nernst potential (E_Na). The relationship between [Na⁺]_i and E_Na, assuming physiological temperature and [Na⁺]_o = 150 mM, is shown in Figure 3A. As apparent, starting from the physiological levels of baseline [Na⁺]_i (about 10–20 mM), each 10 mM increase in [Na⁺]_i results in a reduction in E_Na of about15 mV. Therefore, an increase in [Na⁺]_i due to firing (Figure 1) or synaptic (Figure 2) activity should significantly change E_Na. Naturally, such changes would affect the membrane potential only during episodes where the Na⁺ conductance is large, such as the peaks of action potentials and EPSPs.

To assess these effects, we first examined the somatic spike amplitudes during the 30 Hz train evoked by current pulses. In the AOB mitral cell (Figure 3B), a considerable reduction in spike amplitude is apparent (blue line) compared to a model that does not take [Na⁺]_i changes into account (light blue line). Such amplitude adaptation is widely observed experimentally (Le Bé et al., 2007), and could not be reproduced using simple Hodgkin-Huxley dynamics (Zylbertal et al., 2015). It is explained by a reduction in E_Na at the AIS (Figure 3B, solid black line), rather than the somatic E_Na (dashed line), which does not change along the train. Recovery from this amplitude reduction would follow the [Na⁺]_i profile in the AIS (Figure 1B), meaning that subsequent action potentials generated during several seconds following the train would also have reduced amplitude. It should be noted that a change in E_Na is only one of many possible underlying causes for spike amplitude adaptation, and experimental work will be needed to assess its relative contribution to the phenomenon.

A similar effect is observed when simulating the spike train in the pyramidal cell (Figure 3C), where the amplitude adaptation stems from changes in both the somatic and AIS E_Na. It should be noted that the original pyramidal cell model diverged from the experimental results by not reproducing a spike amplitude adaptation during a test pulse (Hay et al., 2011). As shown here, introducing Na⁺ dynamics to this model bridges this gap between the model and the experimental results.

In the Purkinje cell model (Figure 3D), little change in the somatic or AIS E_Na is inferred from [Na⁺]_i, and the peaks of the fast action potentials are far below E_Na, causing a near-constant spike amplitude, as previously observed experimentally (Llinás and Sugimori, 1980).

In order to examine the effect of E_Na changes on synaptic inputs, we reverted to the synaptic stimulation pattern we described above in the pyramidal cell model (Figure 2), this time using both distal and proximal dendritic sites (Figure 3E, arrow and arrowhead). Figure 3F presents the EPSPs recorded in the soma in the case of a distal stimulation site, when a low Na⁺ diffusion constant was used (approximating dendritic spines). As apparent, the EPSP amplitudes are reduced along the input train, in a manner that resembles weak short-term synaptic depression. Subsequent inputs, at 2 and 9 s after the train, demonstrate the slow recovery of the EPSP amplitude.

The normalized somatic EPSP amplitude (Figure 3G, green line) shows a 17% reduction in amplitude by the end of the train. A smaller reduction was found when the effect of spines was not taken into account (magenta line) or when the proximal stimulation site, which has a larger dendritic diameter, was examined (orange line). In a model without Na⁺ dynamics, no EPSP amplitude reduction is observed (dashed green line).

As in the case of action potentials, the EPSP amplitude reduction is explained by a local change in E_Na. This change is maximal in the distal stimulation site when the effect of dendritic spines was included (Figure 3H, green line), intermediate when the effect of spines was excluded (magenta line), or small when the synaptic input was located at a proximal segment (orange line).

In conclusion, changes in [Na⁺]_i result in significant and possibly prolonged changes in both action potential and EPSP amplitude, especially following bouts of activity. Changes in action potential amplitude may not propagate down the axon, therefore, the degree of their functional significance is not clear. Changes in EPSP amplitude are likely to influence neuronal information processing, reminiscent of synaptic depression.

Elevated [Na⁺]_i Causes an Increase in Na⁺-K⁺ Pump Outward Current, Which May Result in Spike Frequency Adaptation and Prolonged After-Hyperpolarization

Similarly to any enzyme, the electrogenic Na⁺-K⁺ pump rate depends upon the concentration of its substrates, [Na⁺]_i and [K⁺]_o. Since our models assume constant [K⁺]_o, we will demonstrate the effect of changing [Na⁺]_i on the outward current generated by the pump. The pump properties (forward and reverse time constants) were taken from the fitted AOB mitral cell model. Figure 4A shows the relationship between [Na⁺]_i and the pump current per cm² (assuming the pump density is equal to the one used in the axon hillock of the mitral cell). Large changes in the pump current are expected, especially when the [Na⁺]_i level is above 20 mM. The pump can therefore be thought of as a Na⁺-dependent outward current mechanism, similar to the Slick and Slack Na⁺-dependent K⁺ channels (Kaczmarek, 2013). Like these mechanisms, it may act as a homeostatic mechanism by preventing runaway activity.

Such variations in the pump current can have significant effects, particularly when large changes in [Na⁺]_i, occur. We therefore examined the current generated by the AIS pumps during and after an evoked train of action potentials (Figure 4B). This current is a function of both AIS [Na⁺]_i and the density of Na⁺-K⁺ pumps in the AIS, hence it does not perfectly mirror the AIS [Na⁺]_i build-up (Figure 1B). In the AOB mitral cell model, this current keeps intensifying during the train, and reaches a peak value of 400 pA (Figure 4B, blue), while in the Purkinje cell model, a steady-state current with a similar amplitude is quickly established (Figure 4B, red). The strongest current is observed in the pyramidal cell model (Figure 4B, green) due to the combination of high [Na⁺]_i and a high density of pumps.

During a prolonged firing burst, such a current may cause spike frequency adaptation (not evident here, since each spike is evoked by an individual current pulse), but its effect may be masked by multiple conductances, which are active during a firing episode (such as voltage- and Ca²⁺-dependent K⁺ conductances). Yet, due to the relatively slow Na⁺ dynamics, the elevated Na⁺-K⁺ pump-mediated outward current may persist for a long period following the burst (Figure 4C), a time in which most active conductances return to their resting state. This effect is especially evident in the mitral cell and pyramidal cell models (blue and green lines), whereas in the Purkinje cell model, the somatic sink results in rapid restoration of AIS [Na⁺]_i (red line, Figure 1B) and a fast reduction in the pump current.

The immediate result of this modulation of the Na⁺-K⁺ current is a prolonged hyperpolarizing potential that follows each burst of action potentials (Figure 4D, solid lines, vs. dashed lines representing results without modeling Na⁺ mechanisms). Notably, the initial part of such hyperpolarization may be related to Ca²⁺-dependent K⁺ conductance. In the case of the mitral cell model (blue line), the latter part of the response is masked by a prolonged Ca²⁺-dependent inward current. The interplay between these two currents underlies the generation of infra-slow oscillations (Cressman et al., 2009; Zylbertal et al., 2017), where the Na⁺-K⁺ pump current acts as an adaptation current that terminates each firing burst and determines the inter-burst interval duration.

The Coupling between Na⁺ and Ca²⁺ Dynamics Leads to Prolonged Ca²⁺ Transients That May Tag Active Synapses

Another interesting consequence of the changes in [Na⁺]_i is their influence on Ca²⁺ extrusion, particularly via the Na⁺-Ca²⁺ exchanger. This is the main mechanism for removal of a large excess of Ca²⁺, such as the one that develops during a burst of spikes by the activity of voltage-activated Ca²⁺ channels (Fierro et al., 1998). Contrary to the Ca²⁺ pump, the Na⁺-Ca²⁺ exchanger does not rely on ATP hydrolysis for Ca²⁺ extrusion but rather uses the trans-membrane Na⁺ gradient to remove calcium from the cell. As a consequence, elevation in [Na⁺]_i leads to reduction in the rate of Ca²⁺ extrusion by the exchanger, and sometimes even to calcium influx (Schäfer, 2001; Jeon et al., 2003). Assuming that the Na⁺-Ca²⁺ exchanger and the Ca²⁺ pump are the only mechanisms that transfer Ca²⁺ across the membrane (i.e., when voltage-gated calcium channels are closed), [Ca²⁺]_i stable-state is a dynamic equilibrium, where Ca²⁺ flows out through the pump and in through the exchanger. The level of this stable-state [Ca²⁺]_i is therefore a function of [Na⁺]_i that may be inferred analytically (see Methods) when treating the latter as fixed (an assumption justified by its slow dynamics compared to [Ca²⁺]_i). This function is plotted in Figure 5A, with the pump and exchanger parameters fitted for the mitral cell model (solid line), or in cases of a 10-fold increase (dashed line) or decrease (dotted line) in density of the Na⁺-Ca²⁺ exchanger. The [Ca²⁺]_i stable-state is thus determined by the exchanger to pump ratio: a larger ratio results in stronger influence of [Na⁺]_i.

Since [Na⁺]_i in the distal dendrite of the mitral cell model remains elevated for a long period of time following a spike train (Figure 1E), so does the dendritic stable-state [Ca²⁺]_i (Figure 5B, dashed line). Following the spike train, the dendritic [Ca²⁺]_i quickly drops to the vicinity of the stable-state [Ca²⁺]_i, and subsequently follows its slow decay (Figure 5B, blue line). This results in a prolonged [Ca²⁺]_i elevation that, in the case of the AOB mitral cell that expresses Ca²⁺-dependent non-selective cation channels, may cause a prolonged inward current underlying persistent activity (Zylbertal et al., 2015).

In the pyramidal and Purkinje cell models, the dendritic active Na⁺ conductance is too low to account for a similar effect. Synaptic activity may, however, produce comparable, albeit local, elevation in [Na⁺]_i (Figure 2). We therefore tested the result of synaptic input on [Ca²⁺]_i dynamics in two input sites: one is the distal site in the dendritic tuft (Figure 5C, arrow), and the other is a distal site located within the Ca²⁺ “hot zone” of the model (700 μm from the soma, Figure 5C, arrowhead). As before, the stimulation pattern included a 4-s, 5 Hz train of inputs and two subsequent isolated inputs (Figure 5D, black triangles), using the low Na⁺ diffusion coefficient to approximate the effect of dendritic spines.

Stimulating in the original distal site (Figure 5D) results in a prolonged [Ca²⁺]_i elevation (120 nM peak, light green line) generated by the change in local stable-state [Ca²⁺] (black dotted line). Note that the density of Ca²⁺ pump and Na⁺-Ca²⁺ exchanger is too low for the actual [Ca²⁺]_i to closely follow the stable-state [Ca²⁺]_i during the train of inputs (i.e., the light green line is initially below the dotted black line). Later, when local [Na⁺]_i slowly decays and so does the stable-state [Ca²⁺]_i, the actual [Ca²⁺]_i closely follows. Thus, the active synapses elicited a prolonged local increase in Ca²⁺ concentration that may lead to local tagging of the active synapse (Martin and Kosik, 2002). It should be noted that such a tagging mechanism operates without activation of Ca²⁺ conductance.

To examine the effect of additional Ca²⁺ influx through voltage-gated channels, we used somatic current injections concurrent with the initial synaptic inputs to simulate eight calcium spikes, activated by back-propagating action potentials (Figure 5D, red triangles, Larkum et al., 2001; Hay et al., 2011). Since the density of Ca²⁺ channels in the examined compartment is low, the calcium spikes did not significantly alter the [Ca²⁺]_i trajectory (Figure 5D, thick green line). The effect of the calcium spikes alone, without synaptic inputs, on [Ca²⁺]_i is shown by the magenta line.

In contrast, when repeating a similar simulation in the “hot zone” site (Figure 5E), the calcium spikes evoked by somatic current injection had a significant effect on the [Ca²⁺]_i trajectory (thick green line vs. light green line, peaks at 110 and 75 nM, effect of calcium spikes alone shown by magenta line). Notably, since the actual [Ca²⁺]_i is below stable-point [Ca²⁺]_i during the train of inputs, additional [Ca²⁺]_i influx at this time contributes to the elevation in [Ca²⁺]_i for a prolonged time. In this case, the system acts as a coincidence detector, where back-propagating spikes, concurrent with the train of synaptic inputs, result in stronger Ca²⁺ tagging of the active synapses.

A similar effect was observed in a Purkinje cell simulation, where climbing fiber inputs, evoking complex spikes that are associated with global [Ca⁺]_i influx, were used instead of back-propagating Ca²⁺ spikes. A 2-s 5 Hz train of parallel fiber inputs were simulated in a single site of the dendritic tree (Figure 5F, black triangles, location indicated by arrowhead). As in the case of the pyramidal cell, the train of synaptic stimulation alone resulted in a prolonged increase in local [Ca²⁺]_i (Figure 5F, light red line), which initially remained below the stable-state [Ca²⁺]_i, and later followed it (Figure 5F, black dotted line). Climbing fiber stimulation was simulated by stimulating all segments of the dendritic tree thicker than 2.2 μm (Figure 5F, yellow segments of the dendritic tree). These stimulations were given at 1 Hz, concurrent with the parallel fiber input train (red triangles). The contribution of Ca²⁺ influx due to the ensuing complex spikes caused a larger overall [Ca²⁺]_i elevation (52 nM vs. 45 nM peak), a difference that persisted for several seconds following the inputs. The effect of climbing fiber stimulation alone is shown using an orange line. Again, this case represents coincidence detection, where active synapses are more strongly tagged by Ca²⁺ in case of simultaneous climbing and parallel fiber inputs.

In conclusion, the coupling of Na⁺ and Ca²⁺ dynamics by the Na⁺-Ca²⁺ exchanger causes a prolonged global dendritic [Ca²⁺]_i elevation in the mitral cell model, and a local [Ca²⁺]_i elevation in the vicinity of active synapses in the Purkinje and pyramidal cell models. The latter effect could result in tagging of active synapses, intensified in some cases by concurrent back-propagating calcium spikes (in the pyramidal cell model) or complex spikes (in the Purkinje cell model).

The Described Effects on Neuronal Computation Are Preserved under Realistic Distributed Input

To test the validity of the effects shown above under a more realistic setting, we simulated the pyramidal cell under a complex set of inputs mimicking realistic cortico-cortical and thalamo-cortical connections. This was done by simulating 1,000 pre-synaptic “cortical” spike trains that follow a Poisson distribution (average firing rate of 1 Hz). Each train stimulated five excitatory connections uniformly distributed along the dendritic tree. Additionally, we employed 1,000 pre-synaptic “thalamic” spike trains. These were set to background average firing rate of 1 Hz with three epochs representing sensory stimulus, where the average firing rate was elevated to 5 Hz. The thalamic connections were clustered (all five connections were co-localized), and were more densely distributed in proximal dendrites (see Methods).

Figure 6A shows the resulting activity when not taking changes in [Na⁺]_i into account (i.e., conventional model behavior). The cell exhibits constant background firing rate of ~7 Hz, and elevated firing rate during epochs of sensory stimuli (red bars). Introducing [Na⁺]_i dynamics (with a low diffusion coefficient, representing the effect of dendritic spines) under identical pre-synaptic spike trains and contact location (Figure 6B) causes several changes in the cell output. These changes are more pronounced when comparing the firing rate histogram of the two cases, averaged over three runs of the simulation (Figure 6C, red vs. black lines). First, the average background firing rate is lowered due to dendritic Na⁺ accumulation that reduces the synaptic reversal potential. More importantly, epochs of elevated firing rate cause a subsequent prolonged reduction in background firing rate (arrow), as well as in sensory evoked firing rate during the course of the same epoch (arrowhead) and subsequent epochs (asterisk). These effects are slightly reduced in a simulation where the local [Na⁺]_i changes are lower, for example, when the thalamic inputs are not clustered and the diffusion coefficient is not adjusted to account for dendritic spines (Figure 6C, green line).

Another consequence of [Na⁺]_i dynamics is integration of synaptic activity by the local Ca²⁺ level. This is shown in Figure 6D for a basal dendrite with high density of thalamic contacts (Figure 6D, blue line) and an apical dendrite with less thalamic inputs (green line). Somatic and AIS [Ca²⁺]_i changes are smaller and transient (Figure 6D, red and cyan lines respectively). For comparison, [Ca²⁺]_i dynamics under fixed [Na⁺]_i conditions are shown with dashed lines. The [Na⁺]_i dynamics underlying these effects are shown for the same four compartments in Figure 6E. Thus, segments having strong synaptic inputs are more effectively tagged by [Ca²⁺]_i.

To conclude, the effects of [Na⁺]_i dynamics demonstrated above using a “clean” stimulus (isolated spatially and temporally) are also applicable when driving the cell using stochastic and distributed inputs.

Discussion

In this study, we used three realistic neuronal models to explore the predicted dynamics of [Na⁺]_i in various cellular compartments and the possible effects these dynamics might have on neuronal computation. Our results show that [Na⁺]_i can affect neuronal information processing by changes in E_Na, the Na⁺-K⁺ pump current, and the Na⁺-Ca²⁺ exchanger rate. The main effects we encountered are spike amplitude adaptation, apparent synaptic depression, prolonged after-hyperpolarization, and prolonged and localized elevation in [Ca²⁺]_i, which can be used as a tagging mechanism of active synapses.

Functional Significance of the Sodium Accumulation

The effects we describe may have several possible implications on neuronal computation:

EPSP Amplitude Modulation

We show that repeated activation of an individual synapse results in gradual and sustained reduction of the EPSP amplitude due to local Na⁺ accumulation and E_Na reduction (Figures 3E–H). This effect is similar to the effect of synaptic saturation (Bush and Sejnowski, 1994), only acting at much longer time scales. It can be seen as a post-synaptic mechanism of use-dependent synaptic depression, which operates primarily at distal synaptic locations, and reduces the impact of redundant vs. novel inputs. A similar effect was observed when action potential amplitude was measured, but its implication on neuronal computation is not clear. In some cases, however, it might be used as an experimental proxy to measure AIS and/or somatic [Na⁺]_i changes.

Prolonged Activity-Dependent Hyperpolarization

Na⁺ accumulation during activity results in increased activity of the Na⁺-K⁺ pump, leading to hyperpolarization of the membrane potential. Since the bulk of this current is in the AIS, its effect is mainly on the output of the neuron. This outcome may protect neurons from runaway positive feedback by limiting the duration of action potential bursts, and acting as a short-term memory of past spiking activity (Pulver and Griffith, 2010; Zhang and Sillar, 2012).

Prolonged Dendritic Ca²⁺ Elevation

According to our model, the coupling of Na⁺ and Ca²⁺ dynamics will result in prolonged local Ca²⁺ elevations whenever activity evokes a prolonged Na⁺ elevation (Figure 5). This will happen in cases of sufficient density of dendritic voltage-gated Na⁺ channels (as in the AOB mitral cell, Figure 5B), or due to repeated synaptic stimulation (Figures 5C–F). The latter is a possible mechanism for specific tagging of active synapses, where [Ca²⁺]_i elevation could activate cellular pathways related to gene expression and synaptic plasticity (Martin and Kosik, 2002). Such tagging would be gradual, since repeated synaptic activation will result in a further increase in [Ca²⁺]_i, due to the wide dynamic range of [Na⁺]_i. Since Ca²⁺ is an important second messenger, prolonged elevation in its concentration could have several other end results. For example, it could directly gate a K⁺ or non-specific ion currents, thus causing prolonged membrane potential changes.

Our results further show that, in some cases, the amplitude of the long-lasting post-stimulus [Ca²⁺]_i reflects the coincidence of specific synaptic activity and backpropagating action potentials or complex spikes. In the pyramidal cell model, this effect occurs in synapses located near the Ca²⁺ conductance hot-zone (700 μm from the soma in the apical dendrite), triggered by a short train of Ca²⁺ spikes evoked by backpropagating action potentials. A similar effect was shown for the concurrent activation of parallel fiber synapses and the climbing fiber synapse in Purkinje cells, triggered by the global [Ca²⁺]_i influx caused by a complex spike. This effect may contribute to spike timing-dependent plasticity over extended time-scales, since it is evident only when backpropagating or complex spikes follow the synaptic inputs.

The common trait among these possible outcomes is the long time-scale, resulting from the relatively slow Na⁺ dynamics. This reflects the potential role of Na⁺ as a slow chemical integrator of past activity, as intense bouts of activity are needed to effectively change its level, yet its decay to baseline is relatively slow. This contrasts with Ca²⁺, the concentration of which rapidly rises as a function of neuronal activity, and rapidly declines as Ca²⁺ is buffered, stored, and actively extruded from the cytoplasm.

Sensitivity to Chosen Parameters and Model Design

Since our goal in this study was to explore the possible effects of changes in [Na⁺]_i rather than to accurately reproduce experimental findings, some of the parameter values we chose for the Na⁺-related mechanisms were not highly constrained. In particular, we often applied values fitted using experimental results in the AOB mitral cell model to the pyramidal and Purkinje cell models and had to guess the densities of pumps and exchangers in several compartments (see Methods). Some aspects of the original models (e.g., the axonal structure and the distribution of several conductances) were themselves poorly constrained, increasing the uncertainty in the revised models. Most of our conclusions, however, hold with varying levels of confidence despite these uncertainties.

First, most of our conclusions regarding the dynamics of [Na⁺]_i themselves (Figures 1–2) are well constrained, since they rely on the density of voltage-gated or synaptic Na⁺ conductance, cellular morphology, and diffusion coefficient—properties which are adequately constrained (except for the way dendritic spines should affect the apparent diffusion coefficient—hence we present simulation results with and without this correction). The agreement between the SERs calculated using the models and the ones experimentally measured validates the active Na⁺ current dynamics used in the models. Therefore, conclusions regarding changes in E_Na and their effects (Figure 3) are also constrained as a direct consequence of the morphology and the synaptic or voltage-gated current.

Regarding Na⁺-K⁺ pumps, Ca²⁺ pumps, and Na⁺-Ca²⁺ exchangers (Figures 4–5), the specific profile of their effects indeed depends on parameters such as ionic affinities and mechanism densities. Nevertheless, their outcome should at least resemble the one we describe, as long as two assumptions regarding Na⁺ and Ca²⁺ extrusion in neurons are accepted: (a) The bulk of Na⁺ extrusion is carried out by the electrogenic Na⁺-K⁺ pump, which is less significant for the [Na⁺]_i dynamics than diffusion, at least for short time scales (Fleidervish et al., 2010). (b) The bulk of Ca²⁺ extrusion (particularly in cases of large Ca²⁺ excess) is carried out by the Na⁺-Ca²⁺ exchanger (Fierro et al., 1998). These assumptions, combined with our results concerning [Na⁺]_i dynamics, necessarily lead to conclusions similar to ours.

Experimental Findings

This modeling work was inspired by our experimental findings in AOB mitral cells, where we observed several of the effects discussed above by measuring [Na⁺]_i, [Ca²⁺]_i, and membrane potential (Zylbertal et al., 2015). We demonstrated the involvement of the Na⁺-K⁺ pump current in the responses of these cells by observing its voltage independence and its sensitivity to ouabain (Na⁺-K⁺ pump blocker) and replacement of Na⁺ with Li⁺. The effect mediated by the Na⁺-Ca²⁺ exchanger was also inferred using Na⁺ replacement by Li⁺, which does not bind to the exchanger. Moreover, [Na⁺]_i dynamics were directly measured using Na⁺ imaging (SBFI, Rose et al., 1999).

Rose and Konnerth (2001) observed strong and prolonged Na⁺ transients in dendrites and spines of hippocampal pyramidal cells following synaptic stimulation, which were qualitatively similar to the results of our simulation (Figure 2). In their case, NMDA receptors mediated the bulk of Na⁺ influx, hence, contrary to our model, it was only observed when suprathreshold stimulation was applied. As in our model, dendritic voltage-gated Na⁺ conductance hardly contributed to the observed Na⁺ accumulation.

Forrest et al. (2012) used experimental and computational methods to demonstrate the involvement of the Na⁺-K⁺ pump current in shaping the firing patterns of Purkinje cells, including prolonged after-hyperpolarization. Scheuss et al. (2006) showed how activity-dependent [Na⁺]_i elevation modulates Ca²⁺ extrusion by the Na⁺-Ca²⁺ exchanger in hippocampal pyramidal cells, and Kiedrowski et al. (1994) demonstrated that glutamate-evoked Na⁺ elevations impair Ca²⁺ extrusion in cerebellar granule cells.

Eilers et al. (1995) observed AMPA receptor-dependent local fast Ca²⁺ transients at sites of parallel fiber synapse activation in cerebellar Purkinje cells, but showed that the bulk of Ca²⁺ influx originates from voltage-dependent Ca²⁺ channels. Konnerth et al. (1998) applied a train stimulus to the same synapse and observed a prolonged Ca²⁺ transient caused by Ca²⁺ release from internal stores. These two observations are in disagreement with our proposed mechanism, where modulation of stable-state [Ca²⁺]_i (Figures 5C–F) at the site of synaptic stimulation is independent of voltage-gated Ca²⁺ conductance or release from internal Ca²⁺ stores. This discrepancy may be explained by the authors' use of a single pulse stimulation or a single short burst of pulses, producing local change in [Na⁺]_i that is insufficient to produce the effect we propose. A short burst of inputs may mimic the real granule cell discharge pattern in case of an isolated stimulus (Chadderton et al., 2004), but is too short to account for granule cell activity during locomotion (Powell et al., 2015).

Some of the predictions we present tend to elude experimental observation using conventional methods, possibly due to two major reasons: (a) time scale bias: these dynamics are typically longer than those usually considered in electrophysiological or imaging experiments, thus they might end up filtered out or not being noticed. (b) weak stimuli: significant effects require Na⁺ build-up by sufficiently high-frequency or prolonged stimulation, rather than a single action potential or brief synaptic input. Thus, effective experimental examination of the phenomena described by us requires strong stimuli and prolonged recording epochs.

Overall, our results demonstrate that while the assumption of constant [Na⁺]_i holds in short time scales and during transient or low frequency activity, changes in [Na⁺]_i are need to be taken into account in multicompartmental neuronal modeling when considering long time scales, synaptic plasticity, and extensive synaptic input or bursting.

Author Contributions

AZ designed and wrote the simulations and analyzed the data. AZ, YY, and SW wrote the manuscript.

Funding

This work was supported by the Israel Science Foundation (Grant 1350/12 to SW) and the Gatsby Charitable Foundation.

Conflict of Interest Statement

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.

The reviewer MS and handling Editor declared their shared affiliation.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fncom.2017.00085/full#supplementary-material

References