The cable equation (Rall, 1977) describes how the temporal evolution of the electrical potential \(\Phi\) along a passive neuronal process with length l and varying radius a depends on the membrane currents, Ohmic axial currents, and on the membrane capacitance. In this section it is demonstrated that these equation can be extended to include currents induced by the diffusion of intracellular ions. For the resulting equations a explicit algorithm based on finite differences is derived.

Similar conceptual derivations can be found in works such as (Henry et al., 2008), which derived fractional cable equations for dendrites. Previous studies have also employed similar equations to model calcium signals in dendrites (Li, 2023). Electrodiffusion models have been used to describe electric signals in spines (Qian & Sejnowski, 1989). However, the formalisms used in these studies still require complex numerical methods (Breit & Queisser, 2021) or make strong simplifications to the composition of the electrolyte (Lagache et al., 2019). In contrast, the equations and the finite difference scheme presented here allow for simple numerical solvers and accurate solutions independent of the composition of the electrolyte.

The subsequent derivation relies on an analogy between the electrical potential \(\Phi\) and the chemical potential \(\mu _k\) (where k denotes the ion species), incorporating the Nernst-Planck equation (Cohen & Cooley, 1965). The Nernst-Planck equation reads as

$$\vec = -D_k (\vec n_k - \dfrac\vec\Phi).$$

(1)

\(\vec\) describes the directed particle current density (number of particles per unit area per unit time) of an ion-species k with density \(n_k\) (number of particles per unit volume), as a sum of diffusion induced by concentration gradients and drift of ions induced by the electric field \(-\nabla \Phi\). \(D_k\) denotes the diffusion constant, \(z_k\) the charge number, e the elementary charge, \(k_B\) the Boltzmann constant and T the temperature.

In the following a 1D-cable with radial symmetry and without radial dependence of the electrical field and the chemical potential gradient is considered. As indicated by Eq. (1), the total particle current \(j_k\) of ion species k is a result of the diffusion current \(j_\) (electric current induced by the diffusion of ions) and the drift current \(j_\) (electric current induced by a gradient in the electric potential). To transform a particle current into an electrical current, one can simply multiply the particle currents \(j_k\) by \(z_k e\). The axial drift current (charge per unit time induced by the electric field) of an ion species k along a cylinder with a cross-section \(A=\pi a^2\) (and radius a) is therefore:

$$\begin j_ \pi a^2 e z_k = -\dfracn_k\dfrac. \end$$

(2)

Now define

$$\begin I_e =\sum _k j_ \pi a^2 e z_k \end$$

as the total axial electric drift current (charge per unit time induced by the electric field) and

$$\begin r_e=\dfrac \end$$

(3)

as the electrical resistivity for electric drift currents (unit ohm metre). Together with Eq. (2) this leads to

$$\begin I_e = - \dfrac \dfrac. \end$$

(4)

The chemical potential \(\mu _k\) of an ion-species k if given by

$$\begin \mu _k = k_B T \ln \left( \dfrac} \right) . \end$$

(5)

\(n_\) is a reference concentration. The gradient of the chemical potential is simply

$$\begin \dfrac = \dfrac\dfrac. \end$$

(6)

According to the Nernst-Planck equation, the axial electric diffusion current (charge flux per unit time) induced by the diffusion of ion species k can be computed as

$$\begin I_ = - j_ \pi a^2 e z_k = - D_k e z_k \pi a^2 \dfrac, \end$$

(7)

where \(j_\) is the charge flux of ion species k per unit area per unit time, induced by diffusion.

Then one defines the diffusional resistivity for electric diffusion currents (unit ohm metre coulomb)

$$\begin r_ = r_e z_k = \dfrac, \end$$

(8)

in analogy to the electrical resistance \(r_\) and insert it into Eq. (7). Together with Eq. (6) this leads to

$$\begin I_ = - \dfrac} \dfrac. \end$$

(9)

Next, following the derivation of the cable equation, as described in (Dayan & Abbott, 2001), on subdivides the cable along the main axis into N discrete segments of identical size \(\Delta x\). \(x_i\) denote the central points of the segments (Fig. 1B). The electric potential \(\Phi\) along the cable results from a very slight imbalance of positive and negative ions, charging the membrane. The electric potential \(\Phi (x_i)\) of the cable segments at location \(x_i\) can be computed based on the membrane capacitance

$$\begin 2 \pi a \Delta x c_m \Phi = \pi a^2 \Delta x \sum _k z_k e n_k, \end$$

(10)

Here, \(c_m\) represents the specific capacitance of the membrane (capacitance per unit area). To ensure exact electro-neutrality at 0 mV membrane potential, a fixed constant negative background charge with density \(n_0\) is incorporated into the model.

To derive the cable equation, one balances the capacitive current and all currents that enter (or leave) a cable segment of finite length \(\Delta x\) (Fig. 1A). For simplicity, a system without any membrane currents is assumed here.

$$\begin 2 \pi a \Delta x c_m \dfrac =& - \left( \dfrac \dfrac\right) |_ + \left( \dfrac \dfrac\right) |_ \\& +\sum _k\left[ - \left( \dfrac} \dfrac\right) |_ + \left( \dfrac} \dfrac\right) |_ \right] \end$$

(11)

This equation is divided by \(2 \pi a \Delta x\) and the limit \(\Delta x \rightarrow 0\) is taken. Further, the gradient of the chemical potential is replaced, using Eq. (6). Compared to its standard form (Dayan & Abbott, 2001), this leads to a cable equation with an additional diffusion term.

$$\begin c_m \dfrac = \dfrac \dfrac \left( \dfrac \dfrac \right) +\sum _k\left[ \dfrac \dfrac \left( a^2 D_k e z_k \dfrac \right) \right] \end$$

(12)

The same derivation can be carried out for the chemical potential to describe temporal changes in the ionic concentrations along the cable. First one sums over all particle currents entering or leaving the volume (Fig. 1A):

$$\begin \pi a^2 \Delta x \dfrac =& - \left( \dfracz_k e} \dfrac\right) |_ + \left( \dfracz_k e} \dfrac\right) |_\\& - \left( \dfrac z_k e} \dfrac\right) |_ + \left( \dfrac z_k e} \dfrac\right) |_ \end$$

(13)

Then the above equation is divided by \(\pi a^2 \Delta x\) and the limit \(\Delta x \rightarrow 0\) is taken again, to arrive at

$$\begin \dfrac = \dfrac \dfrac \left( \dfrac} \dfrac \right) + \dfrac \dfrac \left( a^2 D_k \dfrac \right) . \end$$

(14)

Equations (12) and (14), together with adequate boundary conditions and initial conditions, uniquely determine the solution. By setting \(\dfrac=0\) one recovers the well known standard form of the cable equation.

2.1 Numerical implementation

For K different ion-species, excluding the constant background charge, e.g. K=3 for \(Na^+\), \(K^+\) and \(Cl^-\), the system of equations derived above (14), together with (12), consist of \(K + 1\) equations. Considering, that \(\Phi\) is a function of \(n_k\) (Eq. (10)), and inserting its spatial derivative

$$\begin \dfrac = \dfrac \sum _k \dfrac. \end$$

(15)

into (14), one has to solve a system of only K equations, effectively.

For numerical treatment, an explicit algorithm based on finite differences can be employed to solve this system of equations. To derive the update rule, one rewrites Eq. (14) with coefficients α and σ as

$$\begin \dfrac = \alpha _ \dfrac \left( \sigma _ \dfrac v_e \right) + \alpha _ \dfrac \left( \sigma _ \dfrac v_ \right) . \end$$

(16)

The two terms \(\dfrac \left( \sigma \dfrac v \right)\) are then approximated based on a central difference scheme with \(\Delta x = h\). After approximating the inner derivative \(\left( \sigma \frac v \right)\), one arrives at

$$\begin \dfrac \left( \sigma \dfrac v \right) \approx \dfrac\left( \dfrac\right) \end$$

(17)

The positions \(x_i\) are indicated by the black dots in Fig. 1B and represent the central points in each segment. h is the length of a segment and the distance between two consecutive points. The interfaces between two segments i and \(j=i+1\) are located at positions \(x_i+h/2\).

In the above equation, \(\sigma (x_i+h/2)\) can be considered as the conductivity of the interface, that connects the segments i and \(i+1\). The conductivity values of the interfaces are the average conductivity between \(x_i\) and \(x_i+h\), and can be computed as the harmonic mean of the conductivities of the connected segments i and \(i+1\)

$$\begin \sigma (x_i+h/2) = \dfrac. \end$$

(18)

As \(\sigma (x_i+h/2)\) is the value on the interface, it can be assumed that \(\sigma (x_i+h/2)\) is constant between \(x_i\) and \(x_i+h\), and one can approximate the outer derivative based on a central difference scheme:

$$\begin \dfrac \left( \sigma \dfrac v \right) &\approx \sigma (x_i+h/2)\dfrac \\&- \sigma (x_i-h/2)\dfrac \end$$

(19)

Next, the left side of Eq. (14) can be approximated based on a forward difference scheme:

$$\begin \dfrac = \dfrac \end$$

(20)

Combining the finite difference approximations of the left side and the right side of Eq. (14) leads to an update rule in which \(n_k\) can be stepped forward in time:

$$\begin n_k(x_i,t_k+\Delta t) = n_k(x_i,t_k) + \Delta t \left( A(x_,x_,x_,t_k) + B(x_,x_,x_,t_k) \right) \end$$

(21)

Here A and B are

$$\begin A(x_,x_,x_,t_k)=&\; \alpha _(x_i) \sigma _(x_i+h/2)\dfrac\\& - \sigma _(x_i-h/2)\dfrac \end$$

and

$$\begin B(x_,x_,x_,t_k) .=&\; \alpha _(x_i) \sigma _(x_i+h/2)\dfrac \\&- \sigma _(x_i-h/2)\dfrac, \end$$

where \(\alpha _=\dfrac\), \(\alpha _d=\dfrac\), \(\sigma _=\dfrac}\), \(\sigma _=a^2 D_k\), \(v_e=\Phi\) and \(v_=n_k\). Code for numerical implementation is available at: https://github.com/feblmu/Eberhardt2024Ion-concentration.

2.2 Boundary conditions and initial conditions

The spatial grid of the computational domain is shown in Fig. 1B. The region with a larger radius a on the right side represents the dendritic end, and on the left side, the synaptic end. The equispaced black dots are the central points \(x_1\) to \(x_N\) of the individual cylindrical segments. Each segment contains values for the electrical potential \(\Phi (x_i, t)\) and the ion concentrations \(n_k(x_i, t)\).

To apply the boundary conditions, one additional segment is appended on the left and right sides (not shown in Fig. 1B). The respective locations are \(x_0\) and \(x_\).

Dirichlet boundary conditions are applied on the dendritic end, where the concentrations \(n_k\) and the electric potential \(\Phi\) at position \(x=x_\) are held constant. The ion concentrations are always held at their initial resting values \(n_k^\). The electric potential \(\Phi\) is set to the resting value \(\Phi ^\) during most simulations. In simulations where the dendritic end gets depolarized, the value is set to a different value \(\Phi (x_)=\Phi ^\).

Neumann boundary conditions are applied on the synaptic end (left side). To prevent flux of potassium and chloride through the synaptic end, the derivatives of the electric potential and the concentrations of chloride and potassium are set to zero at the synaptic end: \(\frac = 0\), \(\frac}=0\) and \(\frac}=0\). The boundary condition \(\frac}=\frac} ez_ \pi a^2)}\) forces the influx of sodium. Applied to the finite difference scheme, sodium influx can be realized by setting

$$\begin n_(x=x_0,t)=n_(x=x_1,t)+h \frac} ez_ \pi a^2)}. \end$$

Here, \(I_\) denotes the input current on the synaptic end of the computational domain (the red line in Fig. 1B), and h represents the spatial discretization, i.e., the distance between the central points of the segments.

As initial conditions at time \(t=0\), the variables were set to the resting values \(\Phi (x,t=0)=\Phi ^\) and \(n_k(x,t=0)=n_k^\).

2.3 Stability analysis

A related finite difference scheme for a linear problem, the parabolic heat equation, indicates that the solution can achieve stability at low spatial and high temporal resolution (Landau et al., 2008). To analyze the accuracy and stability of the presented finite difference scheme, the spatial and temporal resolutions is varied (Fig. 2A-F). As stability is guaranteed by choosing a sufficiently low temporal resolution, the required temporal resolution for a fixed spatial discretization is determined by reducing temporal resolution until a stable solution is found (Fig. 2G-L).

2.4 Simulation parameters

The simulation parameters were chosen as follows:

Specific membrane capacitance: \(c_m = 0.01~F/m^2\)

Temperature: \(T=310~K\)

The resting potential was set to \(\Phi ^=-70~mV\) For the analysis of the voltage-dependent NMDAR-channel current, \(\Phi ^\) was set to \(-85\) mV.

Intracellular ion concentrations at rest: \(n_^=140~mM\), \(n_^=10~mM\), and \(n_^=10~mM\)

The concentration of the fixed negative background charge \(n_0\) was chosen according to Eq. (10) to ensure a resting potential \(\Phi ^\) of \(-70~mV\) (Milo et al., 2010).

Diffusion constants: \(D_=0.65\cdot 10^\dfrac\), \(D_=1.00\cdot 10^\dfrac\), and \(D_=1.00\cdot 10^\dfrac\). Diffusion values are chosen in agreement with (Samson et al., 2003) and (Eberhardt, 2023).

A spine with a total length of 1.4 \(\mu\)m was discretized into 5 head, 5 neck, and 4 dendritic segments of equal length.

In all simulations, a time step of 0.1 ns was used.

The radii of the cylindrical segments in the spine head and the spine neck were varied.

2.5 Analysis

To analyze changes in electrical resistance, the total resistance of the spine for drift currents is computed by summing over all finite segments along the spine:

$$\begin R_e =\sum _i\dfrac , \end$$

(22)

where \(r_e(x_i)\) denotes the electrical resistivity for drift currents in the i-th segment, and \(A=\pi a^2\) represents the cross-sectional area.

Estimations of the spine neck resistance (\(R_n\)) in experiments often rely on a voltage divider model (Cornejo et al., 2022) given by:

$$\begin R_n = \frac - V_}}, \end$$

(23)

where \(V_\) is the experimentally measured membrane potential of the spine head. In the simulations conducted in this study, \(V_\) is determined by the electric potential in the first segment of the spine head, denoted as \(\Phi (x_1)\). Meanwhile, \(V_\) corresponds to the measured membrane potential in the dendritic shaft, estimated here as \(\Phi (x_N)\). The parameter \(I_\) represents the experimentally estimated synaptic current, and its value is controlled by the strength of the injected current (\(I_\)) during the simulations.

To analyze NMDAR-channel currents, the model and all parameters of NMDAR-channel kinetics, but also the resting potential, are based on (Chiu & Carter, 2022). The other spine parameters are identical to those of the simulation shown in Fig. 1. The normalized channel conductance is

$$\begin g_(V) = \dfrac, \end$$

where V is the voltage difference across the synapse \(\Phi (x_1)-\Phi (x_0)\). The normalized NMDAR-current gets approximated by Ohm’s law.

$$\begin I_ = g_(V) (\Phi (x_1,t)-\Phi ^). \end$$

The reversal potential (\(\Phi ^\)) was set to 0 mV, and the parameters a and b were selected as \(a=0.073\) and \(b=-0.074.\)