Location Specificity of Transcranial Electrical Stimulation on Neuronal Electrodynamics: A Mathematical Model of Ion Channel Gating Dynamics and Ionic Flux Due to Neurostimulation

1. Introduction

Transcranial electrical stimulation (TES) is a group of neurostimulation therapies that deliver low doses of electric current to targeted brain regions via noninvasive electrodes placed on a patient's scalp. The most common type of TES is transcranial direct current stimulation (tDCS), which administers a constant amount of electrical energy during therapy sessions. Other forms of TES include transcranial alternating current stimulation (tACS) as well as transcranial random noise stimulation (tRNS), both of which utilize a non-constant dosage of electric current (Paulus, 2011; Antal and Paulus, 2013). Most recently, high-definition TES has been introduced as a neurostimulation approach that achieves a more focused delivery of electrical energy through the use of numerous smaller anode and cathode electrodes, as opposed to just the single larger-sized anode and cathode traditionally used in tDCS, tACS, and tRNS (Borckardt et al., 2012; Caparelli-Daquer et al., 2012).

Clinical experiments clearly show that TES is an effective intervention for treating conditions that manifest from neurodegenerative disorders. Parkinson's disease patients, for example, have demonstrated enhanced movement capabilities and memory skills from TES (Boggio et al., 2006; Johnson et al., 2008). Also, individuals suffering from Alzheimer's disease have demonstrated improved recognition and memory capabilities (Boggio et al., 2009, 2011). Further, TES has shown to improve language re-learning in dementia patients (Wang et al., 2013; Yun et al., 2016). In addition to clinical findings, biological experiments have begun to show the effects of TES on membrane polarization (Liebetanz et al., 2002; Bikson et al., 2004; Stagg and Nitsche, 2011; Das et al., 2016) and calcium homeostasis (Islam et al., 1995; Nitsche et al., 2003; Adams et al., 2017), however difficulties in capturing ion channel state, ionic flux, and intracellular calcium concentrations continuously over time with a high sampling frequency yields limited neurostimulation data at the cellular level (Adams et al., 2017). Thus, the direct influence of an applied TES electric current on voltage-gated ion channel states as well as other cellular level mechanism by which TES operates is largely unknown (Nitsche et al., 2008).

In partnership with biomedical research, mathematical modeling and computational simulation have helped to enhance the neurological communities' understanding of TES. Recent models have begun to describe the impact of electrical stimulation on electric potential around neural tissue (Mandonnet and Pantz, 2011). In addition, biodomain models have provided a means to begin to characterize the influence of electrical energy on transmembrane potential using volume averaging approaches (Sadleir, 2010; Dougherty et al., 2014). These models support the physiological conclusion that TES influences the neuron by slightly polarizing the cell membrane (Nitsche et al., 2008), however, the level of biological abstraction of their mathematical formulations inherently prohibits a quantitative description of individual ion species and their movements around and through the neuron cell wall. A mathematical model of TES that incorporates the electrodiffusion of distinct ion types throughout the intracellular and extracellular domains, as well as their mobility across the cell membrane via voltage-gated ion channels, would for the first time give researchers the capability to computationally assess the impact that a TES-based electric field has on ion channel gating and subsequent ionic flux.

In this paper, we present a novel mathematical model of TES that provides a description of its effects on cellular level neuronal electrodynamics. The model integrates the Poisson-Nernst-Planck electrodiffusion system of partial differential equations (PDEs) and Hodgkin-Huxley motivated boundary conditions for cell membrane ionic flux with extracellular boundary conditions that model TES treatments. Four ion species, namely sodium, potassium, chloride, and calcium, are incorporated in the model. We include calcium in this model as cytosolic calcium is known to be an essential member of the intracellular biochemical network that triggers proper neurotransmitter secretion, and in addition, holds an integral connection with neurodegenerative diseases (Bezprozvanny, 2009; Marambaud et al., 2009; Calì et al., 2014; Surmeier et al., 2017). The TES model is then simulated on a biologically-inspired computational domain (Sosinsky et al., 2005; Chang and Rasband, 2013; Arancibia-Cárcamo et al., 2017) that includes intracellular, extracellular, and membrane regions. Using in silico experiments, we examine the impact of TES on (i) extracellular and intracellular electric potential, (ii) resting membrane potential along the node of Ranvier, (iii) voltage-dependent ion channel gating, (iv) ionic membrane flux, and (v) extracellular and intracellular ion diffusion.

Results demonstrate that a simulated TES current does in fact instantaneously polarize the transmembrane resting potential, however, it does so in a location-dependent manner, where depolarization occurs in a portion of the node of Ranvier and hyperpolarization in other regions. In turn, there is a location-dependent effect on voltage-gated ion channel states, which directly impacts ion channel permeability. Additionally, results show a location-dependent influence on ion membrane flux, with regions along the membrane that exhibit significant increases in sodium and calcium intracellular influx. Of particular importance to applications focusing on neurodegenerative diseases, simulations of TES show calcium intracellular concentrations can increase by up to 71.65% along some regions of the node of Ranvier. In addition, the total calcium concentration in the intracellular domain increases by 63.86% due to TES.

To our knowledge, this paper presents the first model of TES that incorporates neurophysiology with individual ion species and transmembrane ionic fluxes. We hope that the models, simulations, and results presented in this work help expand the research communities' understanding of the neurological mechanisms by which TES operates, and in addition, broadens the utility of mathematical modeling and simulation for computational neuroscience research.

2. Materials and Methods

2.1. Poisson-Nernst-Planck Model

The time-dependent Poisson-Nernst-Planck (PNP) system of partial differential equations (PDEs) can be used to model ion electrodiffusion around and within a neuron (Horng et al., 2012; Dione et al., 2016). The Nernst-Planck equation, which describes particle movement due to both diffusion and electrostatic forces, is given by

$\begin{array}{ll}\frac{\partial n_i}{\partial t}+\nabla ·F_i=0,& (1)\end{array}$

where the ion flux, F_i, is given by

$\begin{array}{ll}{F}_i=-D_i(\nabla n_i+\frac{n_i}{{\alpha }_i}\nabla \varphi ),& (2)\end{array}$

where $n_i=n_i(\overrightarrow{x},t)$ and $\varphi =\varphi (\overrightarrow{x},t)$ represent the concentration of the ith ion and the electric potential energy, respectively, both of which are unknown quantities to be solved for. In addition, constant D_i is the diffusivity in water for the ith ion, and the constant α_i equals $\frac{RT}{F{z}_i}$, where R, T, and F are the gas constant, temperature of the medium, and Faraday's constant, respectively.

The Poisson equation portion of the PNP system quantifies the electric potential energy due to ion concentrations and their relative valences, and is given by

$\begin{array}{ll}\nabla ·(ϵ\nabla \varphi )=-F{\displaystyle \sum _{i=1}}{z}_i{n}_i,& (3)\end{array}$

where z_i is the valence of ion i. In addition, ϵ denotes the permittivity of the medium, equaling ϵ_c · ϵ₀ in intracellular and extracellular regions, and ϵ_memb · ϵ₀ in the cell membrane. Here, ϵ₀ is given by vacuum permittivity while ϵ_c and ϵ_memb are relative permittivities of the intra/extra-cellular and membrane domains, respectively.

In this paper, four ion species are used in the PNP model, namely sodium (Na⁺), potassium (K⁺), calcium (Ca⁺²), and chloride (Cl⁻); thus, equation 1 is realized four times, and the summation term of equation 3 contains four terms.

2.2. Computational Domain

The model is simulated on a biologically-inspired two-dimensional domain representing a portion of a neuron axon that includes a single node of Ranvier, the neuronal region rich in ion channels and transmembrane ionic transport. The domain was constructed using both the myelinated and unmyelinated regions of the membrane, and biologically accurate dimensions were incorporated (Sosinsky et al., 2005; Lopreore et al., 2008; Briegel et al., 2009; Chang and Rasband, 2013; Pods et al., 2013; Maxwell, 2014; Dione et al., 2016; Arancibia-Cárcamo et al., 2017; Rogers and Team of Encyclopedia, 2018). The three subregions of the computational domain consists of (i) intracellular space, (ii) membrane, (iii) and extracellular space.

Figure 1 presents the domain, noting the locations of the three regions as well as all domain boundaries. The length of the axon portion of the domain is 4 μm (Lopreore et al., 2008; Dione et al., 2016) with the nodal portion having a length of 1 μm (Sosinsky et al., 2005; Lopreore et al., 2008; Dione et al., 2016; Arancibia-Cárcamo et al., 2017; Rogers and Team of Encyclopedia, 2018). The radius of the myelinated and unmyelinated sections of the membrane are 0.406 μm (Lopreore et al., 2008; Dione et al., 2016) and 0.005 μm (Briegel et al., 2009; Chang and Rasband, 2013; Pods et al., 2013), respectively. The radius of the intracellular space is 0.434 μm (Lopreore et al., 2008; Dione et al., 2016), and the whole domain, i.e., intracellular, membrane, and extracellular spaces, has a radius of 2 μm (Lopreore et al., 2008; Dione et al., 2016).

FIGURE 1

Figure 1. Diagram of computational domain with intracellular (Ω_I), membrane (Ω_M), and extracellular (Ω_E) subdomains. The diagram also includes labels for each boundary in the domain. Γ_L and Γ_R are the boundaries for the left and right sides of the extracellular space, respectively. Γ₁ is the boundary for the top of the extracellular space and Γ₂ labels the exterior boundaries for the intracellular subdomain. Γ₃ is the exterior boundary of the membrane and Γ₅ labels the boundary between the membrane and intra/extra-cellular space other than in the node of Ranvier, which is labeled by Γ₄.

Figure 2 displays the discretized computational mesh used in each simulation; in this mesh, there are 725,528 elements, with 67,810 nodes in the membrane, 502,644 in the intracellular space, and 159,410 in the extracellular space. The mesh has a much finer grid resolution in the Debye layer, the extracellular space directly adjacent to the membrane, as well as its neighboring intracellular space; this finer discretization is necessary to accurately model the rapid solution changes that take place in these regions of the domain (Pods et al., 2013).

FIGURE 2

Figure 2. Computational mesh with nodes on which the PDEs are solved. Intracellular and extracellular subdomains are shown in blue and the membrane region is shown in red.

2.3. Boundary Conditions

Equation (1) is defined on the intracellular and extracellular regions of the domain, namely Ω_I ∪ Ω_E, whereas Equation (3) is defined on the entire domain Ω = Ω_I ∪ Ω_M ∪ Ω_E (Pods et al., 2013). Thus, boundary conditions for these equations must be stipulated on these respective boundaries. In addition, to appropriately model TES at the cellular level, boundary conditions for the Nernst-Planck equation and the Poisson equation must be specified to model TES administration as well as ion transport across the cell membrane. These conditions are described in the following sections.

2.3.1. TES Boundary Conditions

On the extracellular space top boundary (Γ₁), the concentrations of each ion are set to a constant bulk solution value using the non-homogeneous Dirichlet boundary condition

$\begin{array}{ll}{n}_i=n_i^0,\overrightarrow{x}\in {\Gamma }_1.& (4)\end{array}$

In this work, we focus on a constant stimulation source, i.e., tDCS, and the next two boundary conditions achieve this form of TES. The electric potential on right side of the extracellular space (Γ_R) is maintained at a value of zero using the homogeneous Dirichlet boundary condition

$\begin{array}{ll}\varphi =0,\overrightarrow{x}\in {\Gamma }_R.& (5)\end{array}$

For the first 2 ms of the simulation, the electric potential on the left side of the extracellular space (Γ_L) is set to zero; after this time, TES is simulated by changing this value to 0.1 V (Faria et al., 2011; Gasca et al., 2011; Datta et al., 2013). This TES administration is represented with the time-dependent Dirichlet boundary condition

$\begin{array}{ll}\varphi =\{\begin{array}{cc}0\hfill & :t\le 2\mathrm{\text{ms}},\overrightarrow{x}\in {\Gamma }_L\\ 0.1& :t>2\mathrm{\text{ms}},\overrightarrow{x}\in {\Gamma }_L\end{array}& (6)\end{array}$

Note that alternative forms of TES can be simulated simply by implementing a non-constant value of ϕ once stimulation is activated in Equation (6).

For Equation (1), on all boundaries except the membrane, ion flux is set to zero:

$\begin{array}{ll}{F}_i·\overrightarrow{n}=0,\overrightarrow{x}\in {\Gamma }_2\cup {\Gamma }_3\cup {\Gamma }_5\cup {\Gamma }_L\cup {\Gamma }_R,& (7)\end{array}$

and the charge density flux for Equation (3) is set to zero on all the boundaries not governed by the TES source or ground using the homogeneous Neumann condition

$\begin{array}{ll}ϵ\nabla \varphi ·\overrightarrow{n}=0,\overrightarrow{x}\in {\Gamma }_1\cup {\Gamma }_2\cup {\Gamma }_3.& (8)\end{array}$

2.3.2. Hodgkin-Huxley Gating Equations

The transport of ions across the membrane wall within the node of Ranvier is governed by the non-homogeneous Neumann boundary condition

$\begin{array}{ll}{F}_i·\overrightarrow{n}=f_i^{memb}(n_i,\varphi ,t),\overrightarrow{x}\in {\Gamma }_4,& (9)\end{array}$

where $f_i^{memb}$ is a time and position dependent function that incorporates a Hodgkin-Huxley based model (Hodgkin and Huxley, 1952; Kay and Wong, 1987; Tuckwell, 2012) to quantify sodium, potassium, chloride, and calcium ion flux (see Appendix B). The membrane flux $f_i^{memb}$ uses the transmembrane voltage V = ϕ_I - ϕ_E in its calculation, which is computed at every point along the membrane in the discretized mesh (Dione et al., 2016).

2.4. Numerical Implementation

Equations (1) and (3) are decoupled using the Gauss-Seidal method (Sundnes et al., 2006). The solution algorithm consists of the following steps:

1. Solve equation 3 for ϕ at time step k + 1 given ion concentrations at time step k, $n_i^k$, with boundary conditions given by Equations (5), (6), and (8). Let ϕ^{k + 1} denote this solution.

2. Solve for $f_i^{memb}$ given $n_i^k$ and ϕ^{k + 1} (see section 2.3.2).

3. Solve Equation (1) for n_i, for each ion type, at time step k + 1 given ϕ^{k + 1}, with boundary conditions given by Equations (4), (7), and (9). Let $n_i^{k+1}$ denote these solutions.

The result is numerical solutions of ϕ and n_i at time step k + 1. This iterative sequence is initiated using prescribed intracellular and extracellular initial concentrations of each ion type, and is repeated until the end of the simulation. Within this loop, an inner iteration is used in step 2 to solve the Hodgkin-Huxley system with a smaller time step. This approach ensures the accuracy of the ion flux at the membrane and enables a larger time step for the more computationally intensive PDE solvers in steps 1 and 3. Given that the transmembrane voltage and subsequent flux vary along the node of Ranvier, a different realization of these ordinary differential equations (ODEs) is needed to be solved at every point along the membrane. In this work, the discretized domain generates 1,700 nodes along the membrane, thus the Hodgkin-Huxley ODE system was instantiated and solved for 1,700 times at each simulation time step.

The PDE in step 1 is solved using the finite element method. The PDE system in step 3 is discretized in time using the θ-rule and space using the finite element method (Mardal et al., 2003). The value of θ was set equal to 1, which corresponds to the Backward Euler method, due to its L-stability properties (Hairer and Wanner, 1996). Resulting weak formulations for these equations are presented in Appendix A. The Hodgkin-Huxley ODEs are solved using LSODE (Hindmarsh, 1983; Hairer and Wanner, 1996). This iterative implementation approach enables numerical solvers tailored to each individual equation to be used (Langtangen and Tveito, 2003), as well as individualized time steps for the PDEs and ODEs.

2.5. Computational Tools

The computational domain (Figure 2) was constructed and discretized using GMSH (Geuzaine and Remacle, 2009). The FEniCS computing platform (Alnæs et al., 2015) was used to solve the partial differential equations. This Python based library offers packages to solve finite element weak formulations subject to all boundary and initial conditions. In addition, Python's SciPy library was used to access the LSODE method (Jones et al., 2001).

Given the complexity of the mathematical model and solution approach, an object-oriented implementation of the code was developed. This approach compartmentalizes major modeling components into “classes,” and in doing so, facilitates debugging as each class can be analyzed independently, and in addition, improves code readability. Furthermore, while object-oriented implementations often take more time to design and implement than traditional procedural implementations, a significant advantage of using a class-based structure is its inherent ability to support alternative applications. For example, changes in domain geometry, mesh resolution, TES parameters, or even in the set of ions used can be effortlessly incorporated with virtually no changes to the software (Dougherty and Turner, 2016).

A class for the Nernst-Planck equation incorporates all information needed to solve this equation. This includes its associated weak formulation, diffusivity values, boundary conditions, time steps, and domain information. There are eight instantiations of this class, one for each ion type for both the intracellular and extracellular domains. A separate class is used to solve for $f_i^{memb}$ needed in step 2 of the iterative solution algorithm. There is an instantiation of this class for each of the four ion types. These membrane current classes in turn possess an object dedicated to solving the Hodgkin-Huxley differential equations, which generates solutions for the gating variables m, n, and h (see section 2.3.2). There are 1,700 instantiations of this class, one for each discretized point on the membrane. Information in this Hodgkin-Huxley class is used by the membrane current class to resolve $f_i^{memb}$ along the membrane, which is then used by the Nernst-Planck class via access to the membrane current class.

2.6. Numerical Simulations

A 20 ms simulation of TES was performed via the boundary condition given by Equation (6). A time step of 0.01 ms was selected for the outer iteration of the solution algorithm (section 2.4) as this value is small enough to accurately model the changes in electric potential and ion concentrations (Pods et al., 2013). For solving the inner iteration of step 2, the ODE system was solved with a maximal time step of 0.0005 ms. As described by Equation (6), TES is simulated by changing the Dirichlet boundary condition value from 0 V to 0.1 V on the left boundary of the extracellular space after t = 2 ms, which was selected as this allows concentration gradients and transmembrane ionic flux to achieve equilibrium; this stimulation dosage is consistent with electric potentials achieved during TES sessions (Fregni et al., 2006; Miranda et al., 2006; Datta et al., 2013). This allows the electric potential, transmembrane voltage, ion channel gating variables, ionic flux, and ion concentrations before and after electrical stimulation to be directly compared, thus enabling a direct assessment of the specific impact of TES on neuronal electrodynamics. All simulation parameter values are presented in Table 1, and all values used in the model and simulations are taken from published biomedical literature and previous neuronal mathematical models (Lopreore et al., 2008; Pods et al., 2013; Dione et al., 2016). Simulation run time was approximately 3 days, 18 min, and 9 seconds, on a computer using a fourth generation Intel XEON processor with 3.7 GHz and eight cores.

TABLE 1

Table 1. Simulation parameters.

An iterative implementation and testing approach was used to verify the accuracy of the model implementation. First, individual solvers for the PDEs given by Equation (1) and (3) were constructed and validated against the online PDE solver DiffpackSE (Bruaset and Langtangen, 1997; Langtangen, 2003). Second, the Hodgkin-Huxley ODE model was implemented and verified independently of the PDEs, thus ensuring that changes in intracellular and extracellular electric potential and ion concentrations at the membrane correctly compute gating variable states as well as flux during membrane polarization (Hodgkin and Huxley, 1952; Kay and Wong, 1987; Tuckwell, 2012). Third, these three solvers were integrated into a single solution code using the object-oriented approach as detailed in section 2.5. Fourth, verification of the complete code came by comparing sodium and potassium membrane flux time courses and magnitudes to results from previous PNP modeling implementations (Lopreore et al., 2008; Pods et al., 2013; Dione et al., 2016). Fifth, the transmembrane voltages, intra/extra-cellular ion concentrations, ion channel gating variables, and membrane current fluxes predicted by the complete, fully-coupled model were compared to the isolated Hodgkin-Huxley code (see section 2.5) to validate the accuracy of the fully integrated, coupled implementation used in all simulations. Finally, we draw comparisons between results of the model and those of published medical studies and biological experiments when available.

3. Results

3.1. Transmembrane Voltage Polarization Exhibits Location Specificity

The electric potential energy, ϕ, throughout the neuronal domain at both the beginning and the end of the simulation is shown in Figure 3. Here, changes in both the distribution and magnitude of ϕ from TES are observed. In particular, prior to neurostimulation application, the electric potential distribution is highly symmetric (Figure 3A), however, after TES administration, the domain is highly asymmetric; the majority of high voltage areas are concentrated on the left side of the domain, juxtaposed with the stimulation source boundary condition ϕ = 0.1 V, and electric potential declines more rapidly as the ground boundary is approached (Figure 3B). In addition, the maximum extracellular electric potential value increases by 55.2% from 0.096 V at the start of the simulation to 0.149 V at the end, which due to ionic electrodiffusion, is 49.0% greater than the anode source voltage of 0.1 V. Further, intracellular values for ϕ increase themselves from a minimum and maximum of 0.020 V and 0.026 V to 0.063 V and 0.078 V, respectively.

FIGURE 3

Figure 3. Electric potential energy (ϕ) throughout the computational domain at t = 0 ms (A) and t = 20 ms (B).

Along the neuron membrane, there is a change in transmembrane voltage upon application of electrical stimulation after t = 2 ms (Figure 4). Figure 4A shows the transmembrane voltage throughout the simulation at 11 equispaced points within the node of Ranvier. These points are labeled as a percent based on their position along the node of Ranvier where, for example, 0, 50, and 100% refer to the points on the far left, middle, and far right of the node. The resting transmembrane voltage for each of these points is approximately -70.23 mV. For the point in the center of the node the transmembrane voltage does not change upon stimulation, maintaining its value of -70.23 mV throughout the simulation. For all other points, immediately at stimulation application, there is an instantaneous jump in transmembrane voltage, however, this change depends on the location along the membrane (Figure 4B).

FIGURE 4

Figure 4. (A) Shows the transmembrane voltages due to TES application at equispaced locations within the node of Ranvier. (B) Shows the percent change in the transmembrane voltage due to TES at each point along the membrane; 1.585·10^-6 μm is the far-left, 2·10^-6 μm is the center, and 2.385·10^-6 μm is the far right. A positive percent change indicates depolarization and a negative percent change indicates hyperpolarization.

These results demonstrate the location dependence of changes in transmembrane voltage due to TES. Specifically, transmembrane voltages at points left of center become hyperpolarized, whereas depolarization occurs on the right-hand side. In addition, the magnitude of the polarization from TES administration varies depending on proximity to the edges and center of the node of Ranvier; these values change to a greater degree near the edges as compared to locations near the center. Furthermore, Figure 4B shows that maximum changes in transmembrane voltage do not occur at the extreme edges of the node, but rather at locations situated at 1.64·10^-6 μm and 2.35·10^-6 μm, which correspond to approximately 9 and 91%, both well within the the edge of the node of Ranvier. Interestingly, hyperpolarization occurs for locations on the side with the 0.1 V stimulation source, whereas depolarization occurs on the side adjacent to the ground boundary condition.

In addition to these findings, it is observed that membrane voltage polarization is sustained throughout the TES application, which is consistent with clinical results that show that TES effects persist in sessions consisting of tens of minutes (Miniussi et al., 2008; Nitsche et al., 2008). This sustained increase in neural impulse sensitivity in specific regions of a node of Ranvier permits the TES treatment efficacy recognized by the medical field (Nitsche et al., 2008; Caparelli-Daquer et al., 2012). Our results are also consistent with clincal research that shows that TES has the net effect of increasing neuron excitability by depolarizing to sub-threshold potential (Liebetanz et al., 2002; Bikson et al., 2004; Stagg and Nitsche, 2011; Das et al., 2016). In addition, changes in transmembrane voltage magnitude are consistent with previous mathematical simulations of TES (Dougherty et al., 2014).

3.2. Voltage Gated Ion Channel State Variables Exhibit Location Specificity

The changes in transmembrane voltage due to TES directly impact the behavior of voltage gated ion channels due to changes in their gating variables (Figure 5). Like Figure 4A, Figure 5 displays the values of each gating variable throughout the simulation at the same 11 equispaced points within the node of Ranvier. The location specificity previously observed with transmembrane voltage is also present for the changes in all gating variables.

FIGURE 5

Figure 5. Gating variable values due to TES application at equispaced locations within the node of Ranvier. (A) Shows m, (B) shows n, and (C) shows h.

Prior to stimulation application, m, n, and h show minimal position dependence as their respective values are essentially equal throughout the membrane. For example, before TES application, m is approximately 0.0281 everywhere in the node of Ranvier. When stimulation is applied, changes in m, n, and h become location specific; points where the cell becomes hyperpolarized, i.e., locations between 0 and 50%, result in decreases in m and n as well as increases in h. On the other hand, at sites of depolarization, namely positions between 50 and 100%, m and n increase while h decreases.

Directly corresponding to the locations of maximum change in transmembrane voltage, positions of greatest change in all gating variables also occur off of the membrane edges near 10 and 90%. In addition, the curves of the gating variables are directly associated to the polarized membrane voltages at the same 11 points. In particular, the amplitudes of the gating variable curves correspond to their associated transmembrane voltages, as well as distances between the curves; more precisely, the ranking of each curve of m based on plot amplitude is identical to the ranking of the transmembrane voltage curves, and in addition, the amount of spacing between m curves (Figure 5A) is proportional to the spacing between transmembrane voltage curves (Figure 4A). The same observations apply for n (Figure 5B), and h (Figure 5C) as well with the exception that the ordering is inverted due to characteristics of h that are subsequently discussed.

While the dependence of gating variables on transmembrane voltage is not unexpected, the location specificity of the gating variables due to TES shown here is novel, and in addition, begins to explain how neurostimulation impacts ion channel gating and subsequent ionic flux. Of particular interest in this regard, a clear difference in the shapes, magnitudes and trajectories of the m, n, and h time course curves is observable; the m gating variable changes rapidly, hitting a limiting value early in the simulation, whereas n and h grow more slowly, and fail to reach an asymptotic value within 20 ms. However, m has the lowest amplitude change of the three, with a maximum change of 0.0013, which is only 26.1 and 16.25% of the changes in n and h, respectively.

Figure 6 shows the values of each gating variable at every point along the discretized membrane at seven different simulation times. At t = 2 ms, each gating variable maintains the same value along the membrane as TES application has not yet started; after administration, the value of each gating variable changes over time based upon its location in the membrane. The speed at which m reaches its limiting value is also seen here as the curves for 5, 10, and 20 ms are virtually identical. In contrast, all curves for n and h are visible and continually change throughout the 20 ms simulation. Similar to transmembrane voltage, maximum and minimum values occur approximately at the 9 and 91% locations. Furthermore, it is seen that on the left-half of the node of Ranvier, the m and n probability values are lower than those attained on the right-half, and the opposite is true for h. As will be shown in section 3.3, the time and location dependence of changes in these gating variables as a result of TES has a direct impact on transmembrane ionic current.

FIGURE 6

Figure 6. Gating variable values at each location along the node of Ranvier at simulation times t = 2.0, 2.1, 2.2, 2.5, 5, 10, and 20 ms. (A) Shows m, (B) shows n, and (C) shows h.

3.3. Membrane Ion Flux Exhibits Location Specificity

As the gating variables dictate ion channel permeability, the location specificity observed in transmembrane voltage as well as m, n, and h has a direct influence on ion flux into and out of the neuron. Figure 7 shows the ion flux for sodium, potassium, and calcium over time at the 11 equispaced points within the node of Ranvier. Given the sign convention of the boundary condition governing membrane current (Equation 9), a negative value for flux indicates current coming into the cell from the extracellular space.

FIGURE 7

Figure 7. Membrane flux for sodium (A), potassium (B), and calcium (C) over the course of the simulations for the 11 equispaced points on the node of Ranvier. A negative flux indicates ion flow into the cell from the extracellular space, and a positive value indicates an efflux out of the cell.

Due to passive electrodiffusion forces from the multi-ion environment, as well as a transmembrane voltage not precisely equal to -70 mV, a slight flux of ions across the membrane occurs prior to TES application. Upon activation after t = 2.0 ms, there are significant changes in neuronal flux. For locations on the right-half of the node of Ranvier, where the cell becomes depolarized (Figure 4B), there is an increase in sodium influx (Figure 7A). This is precisely predicted by the gating variable results (Figures 5A,C); as m represents sodium channel activation, which increases on the right-hand side, and h, sodium channel inactivation, which decreases on the right, an increase in sodium influx is this region is expected, and as shown in Figure 7A is attained. In addition, this influx is greatest at the 91% mark, which correlates with all prior results including (i) where the cell experiences its greatest depolarization, (ii) where m is maximal, and (iii) where h in minimal. On the hyperpolarized left-hand side, sodium influx still occurs, but at a decreased rate as m decreases and h increases here.

The m gating variable also controls calcium channel activation (see Appendix B), and so trends in calcium flux function similarly to sodium flux (Figure 7C). Specifically, locations where the cell becomes depolarized yield an increase in calcium influx and hyperpolarized regions experience a decreased influx. For potassium, due to its reversal potential, the opposite occurs and an efflux transpires throughout the entire node of Ranvier. In addition, as n governs potassium activation, potassium efflux increases on the left side where hyperpolarization presents and decreases on the right half of the node of Ranvier (Figure 7B).

These results are consistent with published TES studies that show an increase in calcium influx from a membrane depolarization due to an electric field applied in the extracellular medium (Nitsche et al., 2003; Adams et al., 2017). In addition, like the biological literature, our model predicts that this influx is governed by voltage gated calcium channel permeability (Islam et al., 1995). The novelty of this model is in extending this knowedge to provide a description of how the voltage gated calcium channels within the node of Ranvier operate to achieve this. First, the model allows to see the changes in flux at a greater frequency and with more spacial detail than has been captured with experiments. In addition, the model identifies the gating variable m as driving the changes in flux. Finally, these results reveal a time and spatial based dependence of the gating variable, voltage gated channel activation, and calcium flux.

3.4. TES Causes Intracellular Calcium Dyshomeostasis

As shown in section 3.3, calcium flows into the neuron from the extracellular space at different rates depending on the region within the node of Ranvier (Figure 7C). Thus, over the course of the TES simulation, an increase in intracellular calcium concentration occurs. However, the magnitude and rate of this increase is unknown. Figure 8 shows intracellular calcium concentrations at six simulation time steps. At t = 0 ms, the entire intracellular space has a constant concentration of 10^-4 mM, which is the initial condition for calcium in this domain. Over time, an increase in calcium concentrations from calcium flux due to TES is seen at all subsequent time steps. In addition, for times t > 0, a larger concentration of calcium is noticed at the membrane region, precisely due to calcium influx at the membrane, along with a diffusion throughout the intracellular domain. At the 91% membrane location calcium concentrations increase by 71.65% over the course of the simulation. Furthermore, the total amount of calcium within the intracellular space increases by 63.86% during the course of the simulation. This increase is approximately linear, as can be seen from the color gradients of the intracellular concentration plots.

FIGURE 8

Figure 8. Concentration of calcium in the intracellular space (Ω_I in Figure 1) during the simulation at time t = 0, 2, 5, 10, 15, and 20 ms.

These results are consistent with prior experiments that found an increase in calcium concentration due to an influx of calcium in the presence of electrical stimulation (Islam et al., 1995; Adams et al., 2017). In fact, the values predicted by the model are within one order of magnitude of those shown in electrical stimulation biological studies (Adams et al., 2017). Moreover, the model augments this knowledge by provding a detailed prediction of how, where, and when calcium ion flow into the neuron as described in section 3.3.

4. Discussion

Mathematical modeling and computer-based simulation has shown to be a valuable component in enhancing neurostimulation efficacy as well as providing an instrument for helping the research community learn about the mechanisms by which it operates. While both in silico and biological experimentation have facilitated a greater understanding of neuromodulation, the cellular-level electrodynamics during electrical stimulation treatments still remain highly elusive. To help address this contention, we have presented a novel mathematical model of transcranial electrical stimulation that describes the effect of TES on ion channel dynamics and transmembrane ionic flux. The model is based on the Poisson-Nernst-Planck system of partial differential equations, and to our knowledge is the first that integrates electric potential energy, individualized ion species, voltage-gated ion channels, and transmembrane flux with a medically-based TES induced electric field, and within a biologically-inspired computational domain, showcases how this treatment effects neuronal electrical functioning.

A key finding of this work is the location specificity exhibited by the cell's electrical processes due to TES. In particular, results show that TES polarizes the neuron as expected, however, the degree of voltage change is dependent on the location within a node of Ranvier, a phenomena reported by the deep brain stimulation modeling community (McIntyre et al., 2003, 2004). In turn, results show that the states of the ion channels also exhibit location-dependent changes, which directly impacts membrane flux and subsequent intracellular sodium, potassium, calcium, and chloride concentrations. While the degree and type of electrical polarization is location dependent, these results show that TES effectively elevates resting membrane potential so that ultimately neuron firing is more achievable (Nitsche et al., 2008).

It is well-known that cytosolic calcium is a key element in the intracellular signaling cascade that enables neurotransmitter secretion as well as cell viability. In addition, a disruption to calcium homeostasis is correlated with neurodegenerative disease (Bezprozvanny, 2009; Marambaud et al., 2009; Calì et al., 2014; Surmeier et al., 2017). Our results augment these findings by showing that TES directly alters calcium membrane flux and intracellular calcium concentration via voltage gated calcium channels, by almost 64% over the course of the simulation. These findings may suggest that a possible mechanism by which neurostimulation achieves therapeutic success, in addition to depolarizing the cell, is by altering calcium dyshomeostasis in diseased neurons.

By implementing the simulation software using an object-oriented approach, its utility can be seamlessly extended to other computational studies and future work. Using these tools, we have begun investigating the impact of TES on more biologically complex domains, including one that encompasses three nodes of Ranvier (Figure 9). In addition, we are starting to examine the effect of TES on three-dimensional domains (Figure 10). We have also begun to compare the influence of different forms of neurostimulation, like deep brain stimulation, on transmembrane ionic flux. Finally, we are interested in examining the effects of ionic flux and cytosolic ion concentrations on intracellular signaling pathways that have implications to neurodegenerative disorders.

FIGURE 9

Figure 9. Electric potential energy due to TES throughout a computational domain of a neuron with three nodes of Ranvier.

FIGURE 10

Figure 10. Extracellular concentration of calcium during TES in a three-dimensional computational domain.

Data Availability

The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

Author Contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Funding

ED and KL were supported by the Summer Undergraduate Research Fellowship (SURF) Program of the Rhode Island Institutional Development Award (IDeA) Network for Biomedical Research Excellence from the National Institute of General Medical Sciences of the National Institutes of Health under grant number P20GM103430.

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.

Acknowledgments

We thank Professor Christoph Börgers for helpful correspondence and information about the calcium component of the Hodgkin-Huxley model. Also, many thanks to Abigail Small for helping with GMSH and meshing related activities. We also thank Andrew DelSanto for assistance with three-dimensional result visualization.

Supplementary Material

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

References

Adams, R. D., Gupta, B., and Harkins, A. B. (2017). Validation of electrical stimulation models: intracellular calcium measurement in three-dimensional scaffolds. J. Neurophysiol. 118, 1415–1424. doi: 10.1152/jn.00223.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Alnæs, M. S., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., et al. (2015). The fenics project version 1.5. Arch. Numer. Softw. 3:100. doi: 10.11588/ans.2015.100.20553

CrossRef Full Text | Google Scholar

Antal, A., and Paulus, W. (2013). Transcranial alternating current stimulation (tACS). Front. Hum. Neurosci. 7:317. doi: 10.3389/fnhum.2013.00317

PubMed Abstract | CrossRef Full Text | Google Scholar

Arancibia-Cárcamo, I. L., Ford, M. C., Cossell, L., Ishida, K., Tohyama, K., and Attwell, D. (2017). Node of ranvier length as a potential regulator of myelinated axon conduction speed. Elife 6:e23329. doi: 10.7554/eLife.23329

PubMed Abstract | CrossRef Full Text | Google Scholar

Bezprozvanny, I. (2009). Calcium signaling and neurodegenerative diseases. Trends Mol. Med. 15, 89–100. doi: 10.1016/j.molmed.2009.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Bikson, M., Inoue, M., Akiyama, H., Deans, J. K., Fox, J. E., Miyakawa, H., et al. (2004). Effects of uniform extracellular dc electric fields on excitability in rat hippocampal slices in vitro. J. Physiol. 557(Pt 1):175–190. doi: 10.1113/jphysiol.2003.055772

PubMed Abstract | CrossRef Full Text | Google Scholar

Boggio, P. S., Ferrucci, R., Rigonatti, S. P., Covre, P., Nitsche, M., Pascual-Leone, A., et al. (2006). Effects of transcranial direct current stimulation on working memory in patients with Parkinson's disease. J. Neurol. Sci. 249, 31–38. doi: 10.1016/j.jns.2006.05.062

PubMed Abstract | CrossRef Full Text | Google Scholar

Boggio, P. S., Khoury, L. P., Martins, D. C., Martins, O. E., de Macedo, E. C., and Fregni, F. (2009). Temporal cortex direct current stimulation enhances performance on a visual recognition memory task in Alzheimer disease. J. Neurol. Neurosurg. Psychiatr. 80, 444–447. doi: 10.1136/jnnp.2007.141853

PubMed Abstract | CrossRef Full Text | Google Scholar

Boggio, P. S., Valasek, C. A., Campanha, C., Giglio, A. C., Baptista, N. I., Lapenta, O. M., et al. (2011). Non-invasive brain stimulation to assess and modulate neuroplasticity in Alzheimer's disease. Neuropsychol. Rehabil. 21, 703–716. doi: 10.1080/09602011.2011.617943

PubMed Abstract | CrossRef Full Text | Google Scholar

Borckardt, J. J., Bikson, M., Frohman, H., Reeves, S. T., Datta, A., Bansal, V., et al. (2012). A pilot study of the tolerability and effects of high-definition transcranial direct current stimulation (HD-tDCS) on pain perception. J. Pain 13, 112–120. doi: 10.1016/j.jpain.2011.07.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Briegel, A., Ortega, D. R., Tocheva, E. I., Wuichet, K., Li, Z., Chen, S., et al. (2009). Universal architecture of bacterial chemoreceptor arrays. Proc. Natl. Acad. Sci. U.S.A. 106, 17181–17186. doi: 10.1073/pnas.0905181106

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruaset, A. M., and Langtangen, H. P. (1997). “Diffpack: a software environment for rapid protoptying of PDE solvers,” in Proceedings of the 15th IMACS World Congress on Scientific Computation, Modeling and Applied Mathematics (Berlin), 553–558.

Google Scholar

Calì, T., Ottolini, D., and Brini, M. (2014). Calcium signaling in parkinson's disease. Cell Tissue Res. 357, 439–454. doi: 10.1007/s00441-014-1866-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Caparelli-Daquer, E. M., Zimmermann, T. J., Mooshagian, E., Parra, L. C., Rice, J. K., Datta, A., et al. (2012). A pilot study on effects of 4 x 1 high-definition tDCS on motor cortex excitability. Conf. Proc. IEEE Eng. Med. Biol. Soc. 2012, 735–738. doi: 10.1109/EMBC.2012.6346036

CrossRef Full Text | Google Scholar

Chang, K.-J., and Rasband, M. N. (2013). “Chapter five - excitable domains of myelinated nerves: axon initial segments and nodes of ranvier,” in Functional Organization of Vertebrate Plasma Membrane, volume 72 of Current Topics in Membranes, ed V. Bennett (Oxford: Academic Press), 159–192.

Google Scholar

Das, S., Holland, P., Frens, M. A., and Donchin, O. (2016). Impact of transcranial direct current stimulation (tdcs) on neuronal functions. Front. Neurosci. 10:550. doi: 10.3389/fnins.2016.00550

PubMed Abstract | CrossRef Full Text | Google Scholar

Datta, A., Zhou, X., Su, Y., Parra, L. C., and Bikson, M. (2013). Validation of finite element model of transcranial electrical stimulation using scalp potentials: implications for clinical dose. J. Neural Eng. 10:036018. doi: 10.1088/1741-2560/10/3/036018

PubMed Abstract | CrossRef Full Text | Google Scholar

Dione, I., Deteix, J., Briffard, T., Chamberland, E., and Doyon, N. (2016). Improved simulation of electrodiffusion in the node of ranvier by mesh adaptation. PLoS ONE 11:e0161318. doi: 10.1371/journal.pone.0161318

PubMed Abstract | CrossRef Full Text | Google Scholar

Dougherty, E., Turner, J., and Vogel, F. (2014). Multiscale coupling of transcranial direct current stimulation to neuron electrodynamics: modeling the influence of the transcranial electric field on neuronal depolarization. Comput. Math. Methods Med. 2014, 1–14. doi: 10.1155/2014/360179

PubMed Abstract | CrossRef Full Text | Google Scholar

Dougherty, E. T., and Turner, J. (2016). An object-oriented framework for versatile finite element based simulations of neurostimulation. J. Comput. Med. 2016, 1–15. doi: 10.1155/2016/9826596

CrossRef Full Text | Google Scholar

Faria, P., Hallett, M., and Miranda, P. (2011). A finite element analysis of the effect of electrode area and inter-electrode distance on the spatial distribution of the current density in tDCS. J. Neural Eng. 8:066017. doi: 10.1088/1741-2560/8/6/066017

PubMed Abstract | CrossRef Full Text | Google Scholar

Fregni, F., Boggio, P. S., Santos, M. C., Lima, M., Vieira, A. L., Rigonatti, S. P., et al. (2006). Noninvasive cortical stimulation with transcranial direct current stimulation in Parkinson's disease. Mov. Disord. 21, 1693–1702. doi: 10.1002/mds.21012

PubMed Abstract | CrossRef Full Text | Google Scholar

Gasca, F., Marshall, L., Binder, S., Schlaefer, A., Hofmann, U., and Schweikard, A. (2011). “Finite element simulation of transcranial current stimulation in realistic rat head model,” in Neural Engineering (NER), 2011 5th International IEEE/EMBS Conference on, Cancun 36–39.

Google Scholar

Geuzaine, C., and Remacle, J.-F. (2009). Gmsh: a 3-d finite element mesh generator with built-in pre- and post-processing facilities. Int. J. Numer. Methods Eng. 79, 1309–1331. doi: 10.1002/nme.2579

CrossRef Full Text | Google Scholar

Hairer, E., and Wanner, G. (1996). Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems. Heidelberg: Springer.

Hindmarsh, A. C. (1983). Odepack, a systematized collection of ode solvers. Sci. Comput. 1, 55–64.

Google Scholar

Hodgkin, A. L., and Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117, 500–544. doi: 10.1113/jphysiol.1952.sp004764

PubMed Abstract | CrossRef Full Text | Google Scholar

Horng, T. L., Lin, T. C., Liu, C., and Eisenberg, B. (2012). PNP equations with steric effects: a model of ion flow through channels. J. Phys. Chem. B 116, 11422–11441. doi: 10.1021/jp305273n

PubMed Abstract | CrossRef Full Text | Google Scholar

Islam, N., Aftabuddin, M., Moriwaki, A., Hattori, Y., and Hori, Y. (1995). Increase in the calcium level following anodal polarization in the rat brain. Brain Res. 684, 206–208. doi: 10.1016/0006-8993(95)00434-R

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, M. D., Miocinovic, S., McIntyre, C. C., and Vitek, J. L. (2008). Mechanisms and targets of deep brain stimulation in movement disorders. Neurotherapeutics 5, 294–308. doi: 10.1016/j.nurt.2008.01.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, E., Oliphant, T., and Peterson, P. (2001). SciPy: Open Source Scientific Tools for Python. Available online at: http://www.scipy.org (accessed June 10, 2018).

Kay, A. R., and Wong, R. K. (1987). Calcium current activation kinetics in isolated pyramidal neurones of the ca1 region of the mature guinea-pig hippocampus. J. Physiol. 392, 603–616. doi: 10.1113/jphysiol.1987.sp016799

PubMed Abstract | CrossRef Full Text | Google Scholar

Langtangen, H. P. (2003). Computational Partial Differential Equations: Numerical Methods and Diffpack Programming. Texts in Computational Science and Engineering. Berlin; Heidelberg: Springer.

Langtangen, H. P., and Tveito, A. (2003). Advanced Topics in Computational Partial Differential Equations: Numerical Methods and Diffpack Programming. Lecture Notes in Computational Science and Engineering. Berlin; Heidelberg: Springer.

Liebetanz, D., Nitsche, M. A., Tergau, F., and Paulus, W. (2002). Pharmacological approach to the mechanisms of transcranial dc–stimulation–induced after–effects of human motor cortex excitability. Brain 125, 2238–2247. doi: 10.1093/brain/awf238

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopreore, C. L., Bartol, T. M., Coggan, J. S., Keller, D. X., Sosinsky, G. E., Ellisman, M. H., et al. (2008). Computational modeling of three-dimensional electrodiffusion in biological systems: application to the node of ranvier. Biophys. J. 95, 2624–2635. doi: 10.1529/biophysj.108.132167

PubMed Abstract | CrossRef Full Text | Google Scholar

Mandonnet, E., and Pantz, O. (2011). The role of electrode direction during axonal bipolar electrical stimulation: a bidomain computational model study. Acta Neurochir. (Wien) 153, 2351–2355. doi: 10.1007/s00701-011-1151-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Marambaud, P., Dreses-Werringloer, U., and Vingtdeux, V. (2009). Calcium signaling in neurodegeneration. Mol. Neurodegenerat. 4:20. doi: 10.1186/1750-1326-4-20

PubMed Abstract | CrossRef Full Text | Google Scholar

Mardal, K. A., Sundes, J., Langtangen, H. P., and Tveito, A. (2003). “Systems of pdes and block preconditioning,” in Advanced Topics in Computational Partial Differential Equations:Numerical Methods and Diffpack Programming, Lecture Notes in Computational Science and Engineering, eds H. P. Langtangen and A. Tveito (Berlin; Heidelberg: Springer), 200–236.

Google Scholar

Maxwell, W. (2014). “Nodes of ranvier,” in Encyclopedia of the Neurological Sciences, 2nd Edn., eds M. J. Aminoff and R. B. Daroff (Oxford: Academic Press), 601–604.

Google Scholar

McIntyre, C. C., Grill, W. M., Sherman, D. L., and Thakor, N. V. (2004). Cellular effects of deep brain stimulation: model-based analysis of activation and inhibition. J. Neurophysiol. 91, 1457–1469. doi: 10.1152/jn.00989.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

McIntyre, C. C., Svasta, M., Goff, L. K., and Vitek, J. L. (2003). Uncovering the mechanism(s) of action of deep brain stimulation: activation, inhibition, or both. Clin. Neurophysiol. 115, 1239–1248. doi: 10.1016/j.clinph.2003.12.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Miniussi, C., Cappa, S. F., Cohen, L. G., Floel, A., Fregni, F., Nitsche, M. A., et al. (2008). Efficacy of repetitive transcranial magnetic stimulation/transcranial direct current stimulation in cognitive neurorehabilitation. Brain Stimulat. 1:326. doi: 10.1016/j.brs.2008.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Miranda, P. C., Lomarev, M., and Hallett, M. (2006). Modeling the current distribution during transcranial direct current stimulation. Clin. Neurophysiol. 117, 1623–1629. doi: 10.1016/j.clinph.2006.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Nitsche, M., Liebetanz, D., Antal, A., Lang, N., Tergau, F., and Paulus, W. (2003). Modulation of cortical excitability by weak direct current stimulation-technical, safety and functional aspects. Suppl. Clin. Neurophysiol. 56, 255–276. doi: 10.1016/S1567-424X(09)70230-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Nitsche, M. A., Cohen, L. G., Wassermann, E. M., Priori, A., Lang, N., Antal, A., et al. (2008). Transcranial direct current stimulation: state of the art 2008. Brain Stimul. 1, 206–223. doi: 10.1016/j.brs.2008.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Paulus, W. (2011). Transcranial electrical stimulation (tes – tdcs; trns, tacs) methods. Neuropsychol. Rehabil. 21, 602–617. doi: 10.1080/09602011.2011.557292

PubMed Abstract | CrossRef Full Text | Google Scholar

Pods, J., Schönke, J., and Bastian, P. (2013). Electrodiffusion models of neurons and extracellular space using the poisson-nernst-planck equations—numerical simulation of the intra- and extracellular potential for an axon model. Biophys. J. 105, 242–254. doi: 10.1016/j.bpj.2013.05.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Rogers, K. Team of Encyclopedia (eds.). (2018). Node of Ranvier. Encyclopedia Britannica. Chicago, IL: Encyclopedia Britannica, Inc.

Sadleir, R. (2010). A bidomain model for neural tissue. Int. J. Bioelectromagn. 12, 2–6.

Google Scholar

Sosinsky, G. E., Deerinck, T. J., Greco, R., Buitenhuys, C. H., Bartol, T. M., and Ellisman, M. (2005). Development of a model for microphysiological simulations: small nodes of ranvier from peripheral nerves of mice reconstructed by electron tomography. Neuroinformatics 3, 133–162. doi: 10.1385/NI:3:2:133

PubMed Abstract | CrossRef Full Text | Google Scholar

Stagg, C. J., and Nitsche, M. A. (2011). Physiological basis of transcranial direct current stimulation. Neuroscientist 17, 37–53. doi: 10.1177/1073858410386614

PubMed Abstract | CrossRef Full Text | Google Scholar

Sundnes, J., Lines, G. T., Cai, X., Bjorn, F. N., Mardal, K. A., and Tveito, A. (2006). Computing the Electrical Activity in the Heart. Berlin; New York: Springer.

Surmeier, D. J., Schumacker, P. T., Guzman, J. D., Ilijic, E., Yang, B., and Zampese, E. (2017). Calcium and parkinson's disease. Biochem. Biophys. Res. Commun. 483, 1013–1019. doi: 10.1016/j.bbrc.2016.08.168

PubMed Abstract | CrossRef Full Text | Google Scholar

Tuckwell, H. C. (2012). Quantitative aspects of l-type ca²⁺ currents. Prog. Neurobiol. 96, 1–31. doi: 10.1016/j.pneurobio.2011.09.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Wu, D., Chen, Y., Yuan, Y., and Zhang, M. (2013). Effects of transcranial direct current stimulation on language improvement and cortical activation in nonfluent variant primary progressive aphasia. Neurosci. Lett. 549, 29–33. doi: 10.1016/j.neulet.2013.06.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Yun, K., Song, I. U., and Chung, Y. A. (2016). Changes in cerebral glucose metabolism after 3 weeks of noninvasive electrical stimulation of mild cognitive impairment patients. Alzheimer's Res. Ther. 8:49. doi: 10.1186/s13195-016-0218-6

PubMed Abstract | CrossRef Full Text | Google Scholar