Impact Factor 3.707 | CiteScore 5.1
More on impact ›

Original Research ARTICLE

Front. Neurosci., 19 May 2020 |

The Appearance Order of Varying Intervals Introduces Extra Modulation Effects on Neuronal Firing Through Non-linear Dynamics of Sodium Channels During High-Frequency Stimulations

  • 1Key Laboratory of Biomedical Engineering for Ministry of Education, College of Biomedical Engineering and Instrument Science, Zhejiang University, Hangzhou, China
  • 2Department of Biomedical Engineering, The College of New Jersey, Ewing, NJ, United States

Electrical pulse stimulation in the brain has shown success in treating several brain disorders with constant pulse frequency or constant inter-pulse interval (IPI). Varying IPI may offer a variety of novel stimulation paradigms and may extend the clinical applications. However, a lack of understanding of neuronal responses to varying IPI limits its informed applications. In this study, to investigate the effects of varying IPI, we performed both rat experiments and computational modeling by applying high-frequency stimulation (HFS) to efferent axon fibers of hippocampal pyramidal cells. Antidromically evoked population spikes (PSs) were used to evaluate the neuronal responses to pulse stimulations with different IPI patterns including constant IPI, gradually varying IPI, and randomly varying IPI. All the varying IPI sequences were uniformly distributed in the same interval range of 10 to 5 ms (i.e., 100 to 200 Hz). The experimental results showed that the mean correlation coefficient of PS amplitudes to the lengths of preceding IPI during HFS with random IPI (0.72 ± 0.04, n = 7 rats) was significantly smaller than the corresponding correlation coefficient during HFS with gradual IPI (0.92 ± 0.03, n = 7 rats, P < 0.001, t-test). The PS amplitudes induced by the random IPI covered a wider range, over twice as much as that induced by the gradual IPI, indicating additional effects induced by merely changing the appearance order of IPI. The computational modeling reproduced these experimental results and provided insights into these modulatory effects through the mechanism of non-linear dynamics of sodium channels and potassium accumulation in the narrow peri-axonal space. The simulation results showed that the HFS-induced increase of extracellular potassium ([K+]o) elevated the membrane potential of axons, delayed the recovery course of sodium channels that were repeatedly activated and inactivated during HFS, and resulted in intermittent neuronal firing. Because of non-linear membrane dynamics, random IPI recruited more neurons to fire together following specific sub-sequences of pulses than gradual IPI, thereby widening the range of PS amplitudes. In conclusion, the study demonstrated novel HFS effects of neuronal modulation induced by merely changing the appearance order of the same group of IPI of pulses, which may inform the development of new stimulation patterns to meet different demands for treating various brain diseases.


Deep brain stimulation (DBS) generally refers to utilizing sequences of electrical pulses to treat neurological and psychiatric disorders in the brain (Lozano and Lipsman, 2013). In theory, electrical stimulations on the neuronal membrane could generate various modulation effects on neuronal activity to treat various brain diseases. However, current DBS therapy has only been popularized in treating a few diseases such as movement disorders and refractory epilepsy (Bergey, 2013; Wichmann and DeLong, 2016; Cury et al., 2017). Besides the lack of understanding of the exact pathophysiological mechanisms of many brain diseases, extending DBS treatment to other diseases is limited by a lack of versatility in stimulation paradigms, as well as a lack of understanding of various stimulation effects on neurons and their mechanisms.

Major parameters involved in the design of pulse stimulation paradigms include pulse amplitude, pulse width, and pulse repetition rate (frequency). In particular, the efficiency and efficacy of DBS therapy are strongly dependent on the pulse frequency [or inter-pulse interval (IPI)]. Common DBS therapy utilizes pulse sequences with a constant frequency of around 130 Hz (Cagnan et al., 2017; Dandekar et al., 2018), termed as high-frequency stimulation (HFS). For example, a frequency in 90–185 Hz was proven effective in suppressing tremor of Parkinson’s disease in clinic (Rizzone et al., 2001; Bittar et al., 2005; Dayal et al., 2017). Frequencies in a similar range have also been applied to control epileptic seizures (Cukiert et al., 2011; Fisher and Velasco, 2014; Laxpati et al., 2014).

Besides constant IPI, irregular temporal patterns of stimulation (i.e., sequences with varying IPI) have been exploited to improve the DBS therapies as well as to probe the underlying mechanisms of DBS (Grill, 2018). For example, stimulations with random IPI have emerged as a strategy to suppress the epileptic activity of neural networks (Cota et al., 2009; Mesquita et al., 2011; Medeiros et al., 2012; Santos-Valencia et al., 2019) and to affect central nervous system arousal (Quinkert et al., 2010; Quinkert and Pfaff, 2012; Tabansky et al., 2014). These studies demonstrated a significant effect of temporally random IPI on behavioral outcomes of seizure suppression or CNS arousal evaluated on experimental animals. Also, stimulations with varying IPI of different patterns have been used to suppress tremor (Birdno et al., 2007; Kuncel et al., 2012; Swan et al., 2016). However, the effects of stimulations and mechanisms of action with varying IPI are uncertain yet, limiting the application and development of irregular HFS in clinic.

Previous studies have shown that HFS with constant IPI can only intermittently generate action potentials in neurons because the HFS-induced depolarization block on the neuronal membrane prevents the neurons from firing in response to each stimulation pulse (Jensen and Durand, 2009; Feng et al., 2013; Rosenbaum et al., 2014). Intermittent activation of individual neurons may generate asynchronous firing in neuronal populations (Popovych and Tass, 2014; Feng et al., 2017; Wang et al., 2018). Pulses with a higher frequency (i.e., a shorter IPI) cause a deeper depolarization block, allowing less neurons to be recruited to fire by every pulse, as indicated by smaller evoked population spikes (PSs) during steady-state period of HFS (Feng et al., 2013, 2014). However, a recent study has shown that even with higher instantaneous frequencies for all pulses, stimulations of varying IPI with a random order can occasionally recruit more, rather than less, neurons to fire simultaneously than stimulations with constant IPI (Feng et al., 2019). This suggests that the recruitment of neuronal firing could be affected not only by pulse frequencies but also by the temporal orders of varying IPI. However, the underlying mechanisms are not clear. The purpose of this study is to unveil the cellular mechanisms of neuronal responses to HFS with different temporal orders of varying IPI, which may promote the development of novel stimulation paradigms of DBS for treating more diseases.

Membrane dynamics of ionic channels, such as Na+ and K+ channels, change non-linearly with time and with transmembrane potential (Hodgkin and Huxley, 1952). We hypothesize that a random arrangement of varying IPI could introduce additional effects of stimulation on neuronal population due to non-linear responses of the neuronal membrane to stimulation. To investigate the relationship between dynamic responses of neurons and instantaneous IPI lengths, we performed both rat experiments in vivo and computational modeling on stimulations of efferent axon fibers of hippocampal pyramidal cells with HFS of gradually varying and randomly varying IPI. Antidromically evoked unit spikes and PSs were used to evaluate the neuronal responses directly to pulse stimulations without involving synaptic transmission.

Hippocampal structure contains high density of neurons with clear lamellar organization of neuronal elements, allowing manipulation of distinct stimulation and recording from local cells and their axons as well as evaluation of responses of neuronal populations. Additionally, the hippocampus plays an important role in brain diseases such as epilepsy and Alzheimer’s disease (Fisher and Velasco, 2014; Deeb et al., 2019). Therefore, the responses of hippocampal neurons to electric stimulation have gained increasing attention in DBS research and in clinic. In addition, we investigated the neuronal responses initiated from axons because axons are more prone to excitation by extracellular pulses than other structural elements of neurons (Nowak and Bullier, 1998). The results of this study provide new clues to design effective novel stimulation paradigms with temporal arrangement of IPI for extending the application of DBS.

Materials and Methods

Animal Experiments

The animal experiment was approved by the Institutional Animal Care and Ethics Committee, Zhejiang University. Seven adult Sprague-Dawley rats (250–350 g) were anesthetized with urethane (1.25 g/kg, i.p.) and placed in a stereotaxic apparatus (Stoelting Co., United States). Details of surgical procedures and electrode placements have been reported previously (Feng et al., 2013, 2014). Briefly, a 16-channel array of recording electrodes (#Poly2, NeuroNexus Technologies, United States) was inserted into the hippocampal CA1 region [anterior–posterior, −3.5; medio-lateral, 2.7; dorso-ventral, 2.5]. A stimulation electrode (#CBCSG75, a concentric bipolar electrode, FHC Inc., United States) was inserted in the CA1 alveus [anterior–posterior, −4.8; medio-lateral, 2.7; dorso-ventral, 2.3] for antidromically activating the neurons upstream in the recording site. Based on the clear lamellar organizations of neuronal structures in the hippocampus, the recordings of both unit activity and the specific waveforms of evoked potentials along the 16 channels were used to correct the final positions of the two electrodes. The channel closest to the pyramidal layer that was able to record an evoked-PS with a maximum amplitude was used to analyze the changes of evoked potential during HFS.

Electrical signals collected by the recording electrode were amplified by a 16-channel extracellular amplifier (Model 3600, A-M Systems Inc., United States) with a frequency range of 0.3–5000 Hz. The amplified signals were then sampled by an ML880 Powerlab 16/30 data acquisition system (ADInstruments Inc., Australia) at a sampling rate of 20 kHz/channel. PSs were extracted by a high-pass digital filter with a cutoff of 50 Hz.

Stimulations were sequences of biphasic current pulses with a strength of 0.3 mA and a width of 0.1 ms per phase. To compare the steady-state neuronal responses with those of two temporal arrangements of varying IPI, we utilized two patterns of pulse sequences: a same prelude 20-s period of 100-Hz HFS with constant IPI of 10 ms allowing the neuronal responses to reach steady state, followed by a 10-s period of 100–200 Hz varying IPI either quasi-linearly decreasing from 10 to 5 ms or randomly varying in the range of 5–10 ms. The mean frequency of varying IPI was 133 Hz, similar to the most used DBS frequency 130 Hz. The deceasing IPI was “quasi-linear” because it was designed to have a uniform distribution of IPI counts but not with an exact linear relationship between IPI lengths and time course.

The 10-s sequences of gradually varying IPI and randomly varying IPI were composed of the same set of 1,334 IPIs in the range of 5–10 ms with a sampling resolution of 0.05 ms, that is, 20 kHz, same as the sampling rate of the signal recordings in the experiments. Under this resolution, the 5- to 10-ms range was divided into 101 different lengths of IPI. With a uniform distribution, each IPI length was repeated 13 times in the sequences with an addition of one IPI in every five IPI lengths to get a total of 1,334 IPIs (101 × 13 + 21) with a duration of 10.005 s. The 1,334 IPIs were arranged from the longest to the shortest in the sequence of gradual IPI or were randomized by the MATLAB function randperm() to get the sequence of random IPI. Four different sequences of random IPI were used in the rat experiments, and one of the sequences was used in the computational simulation.

Simulation Models

The computational model includes a population of pyramidal neurons in hippocampus with the mechanisms of potassium ion accumulation and clearance in the narrow peri-axonal space, distributed application of stimulation currents induced by a bipolar electrode, and recording of extracellular potential. The simulation was performed by using the NEURON simulation package (Hines and Carnevale, 1997).

Model of Pyramidal Neurons


The pyramidal neuron consisted of a cell body, an axonal initial segment (AIS), dendrites, and a long myelinated axon (Figure 1A). The cell body was an asymmetric spindle with a length of 20 μm (Migliore et al., 1999), consisting of 10 cylinder segments with a maximum diameter of 20 μm in the center that linearly decreased to 5 μm to the dendrite side and to 1.7 μm to the AIS side. The AIS was 50 μm long with a diameter of 1.22 μm (Hu et al., 2009). The dendrites were simplified as five separate cylinder segments, each with 20-μm length and 5-μm diameter. The myelinated axon consisted of 20 Ranvier nodes (abbreviated as Node) and 20 internode segments with detailed structures similar to those of a previous axon model (Bellinger et al., 2008; Guo et al., 2018). The outside diameter of axon was ∼1 μm wrapped by 15 layers of myelin membranes. The internode segments were made up of paranodal junction (PNJ), juxtaparanode (JXP), and internode (IND). Details of sizes and morphologic parameters of each neuronal element are shown in Table 1.


Figure 1. Schematic diagram of the simulation model. (A) The neuron consisting of a cell body, simplified dendrites, an axonal initiation segment (AIS) and a long myelinated axon. The axon consisted of 20 Ranvier nodes (Node) and 20 internode sections including paranodal junction (PNJ), juxtaparanode (JXP), and internode (IND). Voltage-gated channels of K+ (Kv7) and Na+ (Nav1.6) are on the membrane of Node, whereas fast K+ channels (Kv1.1) are on the JXP membrane. The green lines denote the outwards flow of K+ through Kv1.1 channels, accumulating in the peri-axonal space, back through NaK pump, and diffusing to the outside of Node. (B) Simulation of the electric field generated by the negative phase of a stimulation pulse in COMSOL. A layer of 256 neurons are arranged in parallel with small axial shifts randomly distributed within 20 μm. The stimulation electrode is located 10 μm over the surface of neuronal layer above the central node (Node10) of all neurons, delivering biphasic current pulses (0.1 ms, 0.3 mA). A mimic recording electrode is set at the center of cell body end.


Table 1. Morphologic parameters of the neuron model.

Electrical properties

The neuronal membrane consisted of both passive and active electrical properties. The parameters of passive electrical properties include membrane capacitance, leakage conductance, leakage reverse potential, and intracellular resistivity. The capacitance and leakage conductance were distributed on the entire membrane of the neuron including the myelin. The intracellular resistivity was 300 Ω cm in cell body and 70 Ω cm in axon (McIntyre et al., 2004). The details of the passive electrical parameters are shown in Table 2.


Table 2. Electrical parameters of the cell body, dendrite, and axon.

The parameters of active electrical properties including sodium, potassium, and calcium channels are distributed in different membrane sections (Table 2 and Figure 1A). The ion channels distributed on the neuronal membrane include the followings: voltage-gated Na+, non-inactivating K+ (M current) channel, voltage-gated K+ channel, Ca2+-dependent K+ channel, and high-voltage activated Ca2+ channel for cell body and dendrites; voltage-gated K+ channel (Kv) and voltage-gated Na+ channel (Nav1.2 and Nav1.6) for the AIS (Hu et al., 2009); fast K+ channel (i.e., Kv1.1) for the membrane of JXP in the internode sections; and voltage-gated Na+ channel (Nav1.6) and slow K+ channel (Kv7) for the Node. Specifically, the Nav1.2 subunit was present in proximal AIS, whereas Nav1.6, although present in proximal AIS, was dominant in distal AIS and the axonal Node. The Kv density increased linearly with distance in the AIS to a maximum value of 1,000 pS/μm2. A previous model (Hu et al., 2009) may be referred to for the detailed distribution of Na+ and K+ channels in the AIS. On the Node membrane, the density of Nav1.6 channels was about 6,000 pS/μm2 (Lorincz and Nusser, 2010) and the K+ channel density was 150 pS/μm2 (Röper and Schwarz, 1989; Battefeld et al., 2014).

Accumulating and Clearing Mechanisms of Potassium Ions in the Model

Potassium accumulated in the narrow peri-axonal space during stimulation because of the outflow of K+ from the intracellular space through K+ channels in JXP (Bellinger et al., 2008; Guo et al., 2018). The accumulated K+ was cleared by NaK pumps distributed uniformly on the membrane of internode sections and by axial diffusion of K+ from the peri-axonal space to the extracellular space outside Nodes (indicated by the green arrow in Figure 1A). The axial diffusion of K+ followed the Fick law:

J = D × A × d [ K + ] o d x (1)

where J is the diffusion flux; D is the diffusion coefficient (1.85 μm2/ms); A is the cross-sectional area of the peri-axonal space of 0.025, 0.076, and 0.076 μm2 outside PNJ, JXP, and IND, respectively; [K+]o is the K+ concentration in peri-axonal space; and x is the diffusion distance. The concentration of potassium ion outside the Node (extracellular concentration) was 3 mM.

The clearance of K+ by NaK pump, which pumped out three Na+ and pumped in two K+ per ATP hydrolyzed, was described as follows:

I NaK = I NaK max ( [ K + ] o [ K + ] o + KmK )
× ( [ Na + ] i 1.5 [ Na + ] i 1.5 + KmNa 1.5 ) ( V + 150 V + 200 ) (2)

where V is the membrane potential; [Na+]i is intracellular Na+ concentration; INaKmax is the maximum transport current per unit area, set as 2.46 μA/cm2; and KmK and KmNa are the equilibrium binding constants of K+ and Na+, set as 5.3 and 27.9 mM, respectively.

Models of Stimulation and Recording

Pulse sequences were applied by a stimulation electrode located 10 μm above the top surface of the axial center of an axon layer or approximately over the center Nodes (Node10) of axons (Figure 1A). The layer included 256 axons distributed uniformly in an array of 8 × 32 in a cross section of 70 μm × 310 μm with an inter-axon distance of 10 μm in both spatial orientations (Figure 1B). The axial positions of the parallel axons shifted randomly within a range of −10 to +10 μm, mimicking the lamellar organization of alveus axons. The cell bodies of the neurons were distributed uniformly in an array of 16 × 16 in a cross section of 300 μm × 300 μm and in random one-to-one correspondence with the 256 axons.

The electrical field generated by biphasic pulses was simulated with finite element software COMSOL Multiphysics 5.3 (COMSOL Inc., Sweden), thereby determining the stimulation potentials imposing immediately outside every Ranvier node (Node) of the axonal fibers. Both the electrode parameters (size and material) and stimulation parameters (current biphasic pulses) used in COMSOL simulation were consistent with experimental settings. The stimulation sequences applied in the model were similar to those used in rat experiments, that is, following a prelude of 100-Hz HFS with constant IPI by 100–200 Hz varying IPI of either gradual IPI or random IPI. A recording site was located in the center of the cell body layer to record extracellular potential (Ep) including the extracellular action potentials either from a specific single neuron (i.e., unit spike, abbreviated as spike) or from a population of neurons (i.e., PS). The equation for Ep was as follows (Varona et al., 2000; Guo et al., 2018):

Ep = 1 4 π σ I r (3)

where σ is the extracellular conductivity (0.286 S/m), I is the transmembrane current, and r is the distance between the recording site and each neuronal element generating I, including the compartments of cell body and AISs. Although the transmembrane currents from dendrites may contribute to the shape of the PS waveforms (Varona et al., 2000), we did not take into account the dendrite currents because the exact PS waveform was not critical for the study. We only utilized the PS amplitude, especially the relative changes of PS amplitudes during an HFS sequence, to estimate the relative changes in the number of firing neurons. The major contributions from cell body and AIS were adequate.

Other signals such as intracellular membrane potentials (Vm), concentration of K+ in peri-axonal spaces ([K+]o), and dynamics of Na+ channels (e.g., inactivation variable of Na+ channel, h_Na) were extracted directly from simulation results in NEURON.

Data Analyses for Both Rat Experiments and Modeling

The amplitude of PS was used to evaluate the number of neurons simultaneously generating action potentials (Kloosterman et al., 2001). The PS amplitude was normalized by the mean PS amplitude induced during steady-state of HFS with constant IPI of 10 ms. The mean value of normalized PS amplitudes induced by each pulse during varying IPI was calculated by dividing the sum of PS amplitudes by the pulse number. The range (maximum–minimum) of normalized PS amplitudes was used to describe the variation of PS amplitudes. Because the history of IPIs plays a role in how many neurons are recruited in the PS, the two IPIs immediately preceding a stimulus (“1-back IPI” and “2-back IPI,” termed as IPI1 and IPI2, respectively) were examined. Correlation analysis was implemented to analyze the relationships between the PS amplitude and the preceding IPIs (including IPI1 and ΔIPI = IPI1−IPI2). Student t-test was used to show the statistical significance of the difference between data groups. Additionally, in simulation, the threshold of h_Na at the center node (Node10) for a specific neuron was defined as the minimal value of h_Na at the onset of a pulse that was able to induce an action potential successfully propagating to the cell body.


Different Neuronal Responses to High-Frequency Stimulation With Gradually and Randomly Varying Inter-Pulse Interval in Rat Hippocampus in vivo

Stimulations applied in the alveus of hippocampal CA1 region can antidromically excite a population of CA1 pyramidal neurons without involving synaptic transmissions (Figure 2A1). A single pulse evoked a large PS (Figure 2A2) with a short latency of 1.27 ± 0.09 ms and an amplitude 9.1 ± 1.9 mV (n = 7 rats), just as that evoked by the very first pulse in any pulse sequences. However, after seconds of prolonged 100-Hz HFS with constant IPI, during the steady state, each pulse only evoked a small PS with an amplitude ∼20% of the original amplitude. Afterward, when the IPI gradually decreased from 10 ms (100 Hz) to 5 ms (200 Hz) in a 10-s period, the PS amplitude decreased further to ∼6% of the original amplitude (Figure 2B1). The PS amplitudes induced by gradual IPI were compared with those with constant IPI in a 10-s period immediately before the gradual IPI (Figure 2B2, the switch moment of the two IPI patterns is set as time “0,” with the time duration of −10 to 0 s for constant IPI and 0–10 s for gradual IPI). The normalized PS amplitude decreased from 1 to 0.3 during the 10-s period of gradual IPI. During this period, the IPI decreased quasi-linearly and distributed uniformly in the range of 5–10 ms (Figure 2B3). The PS amplitude significantly correlated with its 1-back IPI (IPI1) with a correlation coefficient of 0.94 (Figure 2B4). When the above stimulation was followed by a pulse sequence with a reverse change of IPI from 5 to 10 ms, a reverse change of PS increase was observed, as expected. For simplicity, the data are not repeatedly presented here.


Figure 2. Neuronal responses to stimulations with different inter-pulse interval (IPI) sequences in rat hippocampal CA1 region in vivo. (A) A1: Schematic diagram of the locations of a recording electrode in pyramidal layer and a stimulation electrode at alveus to antidromically activate the pyramidal cells. A2: A population spike (PS) evoked by a single pulse (0.3 mA). (B) PS evoked by pulses with gradually decreasing IPI. B1: A typical PS recording during 30-s HFS consisted of 20 s of constant IPI (10 ms, 100 Hz) followed by 10 s of gradually decreasing IPI (10 to 5 ms, i.e., 100 to 200 Hz). B2: Changes of normalized PS amplitudes during 10 s of constant IPI and 10 s of gradual IPI. The switch moment from constant IPI to gradual IPI is set as time “0.” B3: IPI curve along stimulation time (left) and distribution probability of gradual IPI (right). B4: Correlation between normalized PS amplitude and preceding IPI during 10-s period with gradual IPI. (C) PS evoked by pulses with random IPI. C1–C4: Corresponding to B1–B4. Bottom of C4: Correlation between the normalized PS amplitude and the difference of IPI1–IPI2.

When the varying IPI of a gradual order was replaced by a random order, the normalized amplitudes of evoked PS varied randomly in a much wider range from 0 (no PS) to ∼2.4 (Figures 2C1,C2), as compared with the range of 0.3–1 induced by gradual IPI. Although the random IPI was in the same range of 5–10 ms and with the same uniform distribution (Figure 2C3), the correlation coefficient (0.73) between the PS amplitudes and the preceding IPI1 of random IPI was smaller than that of gradual IPI. The less correlation was caused by a lot of small PSs appearing following the longer IPIs (Figure 2C4). However, the correlation coefficient between the PS amplitudes and ΔIPI (=IPI1−IPI2) increased to 0.88 (Figure 2C4), indicating that a longer IPI1 accompanied by a short IPI2 (2-back IPI) would result in a larger PS with a higher probability.

Statistical data (n = 7 rats) showed that the mean normalized PS amplitude during 10-s period of random IPI was similar to that of gradual IPI (Table 3), whereas the range of normalized PS amplitude with random IPI was significantly greater than the value with gradual IPI. However, because of the same pulse number and similar mean of normalized PS amplitude, the sums of normalized PS amplitude during gradual IPI and random IPI were approximately equal, indicating similar amounts of total neuronal firing. Additionally, the mean correlation coefficient between the PS amplitude and the preceding IPI1 (or IPI1−IPI2) was significantly smaller with random IPI than with gradual IPI. The results of large PS indicated that pulses with random order of varying IPI could occasionally recruit a larger fraction of neurons to fire together.


Table 3. Normalized PS amplitudes distribution during gradually and randomly varying IPI.

Previous studies with constant IPI have shown that prolonged HFS can extend the refractory periods of the neuronal membranes owing to accumulation of extracellular potassium, leading to the intermittent firing of individual neurons and a decrease of PS amplitude with higher frequency (i.e., shorter IPI) (Feng et al., 2014; Guo et al., 2018). We hypothesize that the same mechanism of potassium accumulation in the narrow peri-axonal space together with non-linear dynamics of ionic channels may generate the irregular neuronal responses with simultaneous recruitment of a larger fraction of neurons by pulses with random order of varying IPI. Because current techniques of in vivo experiments have limitations in simultaneous intracellular recordings of individual axons, we next tested the hypothesis by a modeling study.

Simulated Neuronal Responses to High-Frequency Stimulation With Varying Inter-Pulse Interval of Gradual and Random Orders

The simulated evoked PS (Figure 3) closely matched the experimental recordings (Figure 2). At the onset of a pulse sequence, the evoked PS was also as large as a PS evoked by a single pulse (Figures 3A,B1). Following seconds of prolonged 100-Hz HFS, during the steady-state period, the PS amplitude decreased to ∼20% of the original amplitude. Afterward, when IPI gradually decreased from 10 to 5 ms in a 10-s period, the PS amplitude further decreased to ∼4% of the original amplitude (Figure 3B1). The normalized PS amplitude decreased from ∼1 at the start (corresponding to 10-ms IPI) to ∼0.2 at the end (corresponding to 5-ms IPI) during the 10-s period with gradual IPI (Figure 3B2). In addition, the PS amplitude was significantly correlated with its preceding IPI with a correlation coefficient of 0.92 (Figure 3B3 and Table 3).


Figure 3. Simulation of neuronal responses to stimulations with different inter-pulse interval (IPI) sequences. (A) Schematic diagram of the locations of stimulation electrode and recording electrode together with a population spike (PS) waveform evoked by a single pulse. (B) Simulated PS evoked by pulses with gradual IPI. B1: Simulated PS signal during 30-s HFS including 20 s of constant IPI (10 ms, 100 Hz) followed by 10 s of gradually decreasing IPI (10 to 5 ms, i.e., 100 to 200 Hz). B2: Changes of normalized PS amplitudes during 10 s of constant IPI and 10 s of gradual IPI. The switch moment from constant IPI to gradual IPI is set as time “0.” B3: Correlation between normalized PS amplitude and its preceding IPI during 10-s period with gradual IPI. (C) Simulated PS evoked by pulses with random IPI. C1–C3: Corresponding to B1–B3. C4: Correlation between the normalized PS amplitude and the difference of IPI1–IPI2.

Similar to the experimental results, when the varying IPI in a gradually decreasing order was replaced with a random order, the normalized PS amplitude varied in a much wider range from 0 to 2.95 (Figures 3C1,C2), more than double of the value with gradual IPI (2.95 vs. 1.12, Table 3). However, the mean normalized PS amplitudes during the periods with two different orders of varying IPI were relatively close (0.62 vs. 0.73). Additionally, the correlation coefficient between the PS amplitude and the preceding IPI1 was smaller with random order IPI than with gradual IPI (0.73 vs. 0.92); but the correlation coefficient between the PS amplitude and the ΔIPI increased to 0.8 (Figures 3C3,C4 and Table 3).

The PS only reflects the summation of firing from a population of neurons. To investigate the underlying mechanisms, we next examined unit spike firing of individual neurons.

Simulated Membrane Dynamics of Individual Neurons During High Frequency Stimulation With Different Inter-Pulse Interval Patterns

Neurons at different distances from the stimulation electrode had different periodical firing patterns because they received different stimulation intensities (Guo et al., 2018). We took one neuron as an example, with its middle axonal node (Node10) located at a distance 104.4 μm (30 μm below and 100 μm left) from the stimulation center (Figure 4A). During the steady state of 100 Hz of constant HFS, the inter-spike interval (ISI) of the neuronal firing kept constant at two times the IPI as the neuron fired a spike at every other pulses (Figure 4B). During the period of gradually decreasing IPI, the firing ISI exhibited two alternating stages: gradually decreasing stage and jumping stage (Figure 4C). At a gradually decreasing stage, the ISI decreased gradually following the decrease of IPI as the neuron fired steadily at every n-th pulses. The stage was interrupted by sporadic and short jumping stages when the ISI increased by one or more IPI. Therefore, during most of the stimulation period with gradual IPI, the firing of individual neurons was similar to the firing under constant IPI. However, during the entire 10-s stimulation period with random IPI, the firing ISI changed randomly in a wider range (14.5–64.5 ms) (Figure 4D).


Figure 4. Simulation of the response firing of a typical neuron to different pulse sequences. (A) Schematic diagram of the neuron’s location with its axon center (Node10) 104.4 μm (30 μm below and 100 μm left) from the bottom center of the stimulation electrode. (B) During the steady-state period of constant inter-pulse interval (IPI) of 10 ms, the neuron fires with a constant inter-spike interval (ISI) of 20 ms. (C) During the stimulation period with gradual IPI, the ISI of neuronal firing gradually decreases along with the decrease of IPI in most time except several jumping stages (denoted by shadows) to elevate the ISI in stairs by adding one or more IPI into ISI. (D) During the stimulation period with random IPI, the ISI of neuronal firing changes irregularly all the time.

To analyze the mechanisms underlying the behaviors of individual neurons during HFS with gradual IPI and random IPI, we extracted the following variables from the simulation results: the accumulated [K+]o in peri-axonal space near Node10, the membrane potential (Vm) on Node10, the inactivation variable (h_Na) of Na+ channel on Node10, and the extracellular spike signal of the soma (Figure 5A).


Figure 5. Simulated dynamics of neuronal membrane during high-frequency stimulation (HFS) with different patterns of inter-pulse interval (IPI). (A) Schematic diagram of simulated signals collected from a neuron, including K+ concentration in peri-axonal space near Node10 ([K+]o), membrane potentials of Node10 (Vm), inactivation variable of Na+ channels on Node10 (h_Na), and spike recording near soma. (B) Dynamics of four signals during stimulations with gradually varying IPI (B1) and randomly varying IPI (B2). Circles on the curves of Vm, [K+]o and h_Na denote the values at the onset of stimulation pulses. Thin horizontal lines on the three curves denote thresholds for neuronal firing. Three typical IPI sequences (shadowed) generate different lengths of inter-spike interval (ISI) of 16.0, 27.1, and 34.8 ms. (C) Generation of different ISI. C1: The three examples of IPI sequences from B2 show that the levels of inactivation variable h_Na at the onset of pulses determine whether AP can be induced at Node10 and whether the AP can successfully propagate to soma. The red dots and cyan numbers on the h_Na curves denote h_Na values at the onset of pulses inducing AP at Node10. C2: A pulse (P2) failing to induce AP can decrease h_Na and delay the recovery course of h_Na, thereby causing no spike at soma following the next pulse (P3). A spike would otherwise be evoked by P3 in the situation without P2. C3: Even an AP is induced on Node10 immediately under the stimulation electrode, due to an inadequate recovery of h_Na, the AP may not be strong enough to propagate along the axon to generate a spike at soma. Hollow triangles in B2 and C1 denote that the APs of Node10 fail to generate a somata spike but can cause a substantial decrease of h_Na.

Following the prelude 20-s period of 100-Hz HFS with constant IPI of 10 ms, the [K+]o had increased from original 3 to above 12 mM. The accumulated K+ elevated Vm from original −70 to around −60 mV, thereby inactivating Na+ channel as h_Na declining from baseline value 0.95 to below 0.5. Nodal membrane depolarized upon each incoming pulse and repolarized during the interval before the incoming of next pulse.

During most of the period with gradual IPI, pulses of HFS changed the membrane dynamics of the neuron periodically. Figure 5B1 shows an episode of gradually decreasing IPI of ∼10 ms. At every other pulse, an action potential (AP) was generated and then resulted in a pronounced elevation of [K+]o as well as a pronounced decline of h_Na, consequently. Accordingly, the Vm, [K+]o and h_Na changed periodically. At the onset moment of each pulse, Vm on Node10 was either −59.0 or −60.6 mV, [K+]o 12.9 or 12.4 mM, and h_Na 0.330 or 0.454, respectively corresponding to a subthreshold state or a suprathreshold state of firing. The alternation between the two states gave rise to an intermittent and regular firing of spikes with an ISI of ∼20 ms, two times the IPI.

In contrast, HFS with random IPI (5–10 ms) induced irregular firing of spikes in the same neuron owing to the non-linear dynamics of the membrane excitation. Because the inactivation of Na+ channel (h_Na) has a much longer time constant than the activation of Na+ channel (Hodgkin and Huxley, 1952), the excitability of Na+ channels during HFS was mainly determined by the recovery of h_Na. A small difference in h_Na could change the membrane excitability and flip its state from subthreshold to suprathreshold or vice versa, thereby drastically changing the firing pattern of the neuron (Figure 5B2). For instance, for the three firing patterns of type 1, type 2, and type 3 illustrated in Figures 5B2,C1, one, three, and four pulses failed to activate spikes in soma and resulted in an ISI of 16.0, 27.1, and 34.8 ms, respectively. Failed pulses still caused the membrane depolarization that delayed the recovery course of h_Na. In this case, the h_Na threshold for generating a spike at soma was 0.386 (Figure 5B2). For the type 2 firing pattern shown in Figure 5C1, with a preceding pulse P2 (a failed pulse), a recovery of h_Na to 0.370 at the onset moment of pulse P3 only allowed the pulse to generate an AP at Node10, but the AP failed to activate the soma to fire a spike. Otherwise, in a case without the preceding failed pulse P2, the h_Na would have recovered to 0.421 to exceed the threshold, thereby driving the soma to fire a spike (Figure 5C2).

It should be noted that a fraction of AP induced at axonal nodes were not strong enough to travel to the soma and generate an AP (i.e., a spike) at soma (see Figures 5B2,C1). Nevertheless, the failed AP would decrease h_Na substantially, delaying the recovery of h_Na and possibly preventing the next pulse from generating a spike (denoted by hollow triangles in Figures 5B2,C1). The repeated delays of h_Na recovery resulted in a prolonged interval between two adjacent spikes (Figure 5C3).

These simulation results indicated that stimulations with regular IPI (including most of the period with gradually varying IPI) can generate AP with regular ISI in individual neurons. However, stimulations with a random order of varying IPI can postpone the firing of individual neurons and generate AP with various ISI. The different patterns of ISI could result in different numbers of neurons to fire spikes simultaneously. Therefore, we next investigated the firing of neuron populations by modeling.

Comparing Firing of Neuronal Population Induced by High-Frequency Stimulation With Regular and Random Inter-Pulse Interval

During episodes of stimulation with gradually varying IPI, individual neurons fired spikes in a fixed periodical pattern. For example, five neurons with different distances (at Node10) to the stimulation center of 28, 70, 104, 117, and 130 μm had different periodical patterns of spike firing (Figure 6A1). Neurons located nearer to the stimulation electrode received a greater stimulation intensity and had a lower h_Na threshold for firing. Because of the overstimulation received, the nearest neuron (Neuron-1) was not able to fire spike, with its h_Na at Node10 only fluctuating below the firing threshold at the onset of each incoming pulse. The three neurons at other distances, Neuron-2, Neuron-3, and Neuron-4 fired spikes at a same induction ratio of 50% (i.e., every other pulse induced a spike), but Neuron-2 and Neuron-3 followed the same pulses whereas Neuron-4 followed the others. The farthest neuron (Neuron-5) had an induction ratio of 33%.


Figure 6. Simulated responses of multiple neurons to form population spikes (PSs) during stimulations with different patterns of inter-pulse interval (IPI). (A) Somata spikes and inactivation variable h_Na curves of five neurons with different distances (28–130 μm) from the stimulation electrode. Thin horizontal lines on the h_Na curves denote spike thresholds of individual neurons (0.280–0.538). A1: During an episode of stimulation with gradually decreasing IPI, except that Neuron-1 does not fire spikes, the other four neurons generate spikes regularly with h_Na regularly switching between suprathreshold and subthreshold levels (denoted by green and black dots, respectively) upon a pulse incoming. A2: During stimulation with random IPI, all of the five neurons generate spikes irregularly with h_Na surpassing threshold occasionally upon a pulse incoming. (B) During an episode of stimulation with gradually decreasing IPI, regular firing of individual neurons with different induction ratios (ratio of spike number to pulse number) produces uniform PS. (C) During stimulation with random IPI, a large PS formed by synchronous neuronal firing is induced by a typical IPI sequence. C1: Different spike sequences of individual neurons integrate into PS sequences with various PS amplitudes including a large PS (yellow zone) with surrounding silent periods of no PS (gray zones). C2: Inactivation variable h_Na curves show the recovery course of h_Na before synchronous firing of neurons following an IPI sequence including two shorter IPI (6.75 and 6.35 ms) and one longer IPI (9.6 ms). The thickening and dark-blue segments of h_Na curves show that the series of preceding pulses fail to generate somata spikes in individual neurons before the generation of the large PS.

However, during the stimulation with random IPI, the five neurons all fired spikes irregularly with their h_Na at Node10 exceeding respective threshold occasionally (Figure 6A2). Interestingly, Neuron-1 changed from no firing with gradual IPI into intermittent firing with random IPI.

The differences in spike firing of individual neurons resulted in different responses of a neuronal population to different patterns of varying IPI. During an episode of stimulation with gradually varying IPI (Figure 6B), the periodical spikes of individual neuron firing were integrated to form PSs with similar amplitudes following every pulse. However, during stimulation with random IPI, owing to the irregular firing of individual neurons, the amplitudes of integrated PS varied randomly in a large range (Figure 6C1). Repeatedly decreasing h_Na by several successive pulses with shorter IPI could postpone the time for h_Na to recover to its threshold (Figure 6C2, thick blue lines in h_Na curves), finally giving rise to simultaneous firing of a large number of neurons to form a large PS following a longer IPI (the yellow zone in Figure 6C). The large PS was accompanied by very small PS induced by neighboring pulses (the two gray zones in Figure 6C). The mean PS amplitudes induced by gradual IPI and by random IPI were similar (Table 3 simulation), indicating that randomization of IPI could only redistribute the firing timing of neurons but not significantly change their firing amount.

These simulation results indicated that the stimulation pulses with random IPI could recruit neurons to fire at specific moments because of the occasional occurrence of specific sub-sequences of pulses, thereby generating extremely large and extremely small PS events.


The novel finding of the present study with experiments and computation is that the HFS pulses with varying IPI in random order had additional recruitment effects on the firing of neuronal population. Possible mechanisms underlying the experimental finding are unveiled by our simulation study and analyzed below.

Possible Cellular Mechanisms Underlying the Firing Patterns of Neurons Induced by Different Patterns of Varying Inter-Pulse Interval

Previous experiment studies both in vitro and in vivo have shown that HFS with a constant frequency over 50 Hz can induce partial depolarization block of axons in the brain regions such as the hippocampus and subthalamus (Jensen and Durand, 2009; Zheng et al., 2011; Rosenbaum et al., 2014). Our previous study has shown that with a substantial decrease of antidromic PS induced by HFS applied to efferent fibers of CA1 neurons as shown in the present study, the neuronal somata preserved the capacity to respond to an orthodromic stimulus with large PS (Feng et al., 2013), thereby confirming the HFS-caused failures in the axons, not in the soma. In vitro studies have shown that repeated pulse stimulation can elevate [K+]o, thereby causing a decreased speed of axonal conduction until a complete block. Artificially increasing [K+]o can also lead to axonal block (Förstl et al., 1982; Poolos et al., 1987). Computational models on HFS with constant IPI suggest the underlying mechanisms of the block: the outflow of K+ during stimulations may be trapped in the narrow peri-axonal space and increasing [K+]o, thereby prolonging depolarization of membrane and inactivation of Na+ channels and leading to intermittent depolarization block (Bellinger et al., 2008; Guo et al., 2018). HFS at the commonly used frequency range (∼100–200 Hz) for DBS can only induce partial but not complete axonal block. The axons can still intermittently generate action potentials to form PS events (Jensen and Durand, 2009; Feng et al., 2013, 2014). During an HFS sequence with a higher frequency (i.e., shorter IPI), each pulse can recruit fewer neurons to fire, thereby only generating a smaller PS (Guo et al., 2018).

The above mechanisms may explain the IPI-dependent PSs evoked by pulses with gradually varying IPI in the present study, as shown by the highly linear correlation between PS amplitudes and lengths of immediate preceding IPI (1-back IPI) in both the experiment and simulation results (Figures 2B, 3B and Table 3). Based on these mechanisms, it seems reasonable to speculate that the IPI-dependent PSs would have also appeared with IPI arrangements in any other orders as long as the IPI set was the same. However, our study shows that even with the same set of IPI varying only in a small range of 5–10 ms (i.e., 200–100 Hz), substantially different responses of population neurons appeared during HFS with random IPI (Figures 2C, 3C and Table 3). Here, we reveal the mechanisms of these different neuronal responses by the recovery course of Na+ channel with non-linear dynamics.

Owing to the much longer time constant of Na+ channel inactivation than that of Na+ activation (Hodgkin and Huxley, 1952), under the situation of intermittent depolarization block, the excitability of Na+ channels in response to HFS pulses would be mainly determined by the non-linear recovery of Na+ inactivation (h_Na). An IPI in the range of 5–10 ms was not long enough for a complete recovery of h_Na after a firing spike under the situation of HFS. It took several IPIs for h_Na to recover to a suprathreshold level. In addition, each failed pulse would also lower the h_Na and postpone the recovery of h_Na (Figures 5, 6). The firing history of a neuron, the lengths of several preceding IPIs, and the pulse impulses following the previous firing affect the non-linear recovery course of h_Na, hence the next firing of the neuron.

During HFS periods with regular IPI (e.g., gradually varying IPI or constant IPI), although the number of IPI cycles needed for the h_Na recovery varied for neurons at various distances from the stimulation point, the firing of individual neurons was regular at most time because of the approximately fixed IPI length within a relative short episode. The integration of regular firing from a large population of neurons remained steady over time (Figures 6A1,B). However, when the varying IPI appeared randomly, the number of IPI cycles as well as the sum of IPI lengths for h_Na to recover to suprathreshold was uncertain, and the firing of individual neurons became irregular, thereby occasionally generating extremely large and extremely small PSs as a result of the firing integration of a neuronal population (Figures 6A2,C).

Taken together, through the mechanisms of intermittent [K+]o elevation and Na+ channel inactivation induced by HFS, as well as the non-linear recovery course of Na+ channel, HFS with different patterns of varying IPI in a small range can modulate the firing of neuronal populations differently. To our knowledge, this is the first simulation study addressing the effect of irregular pulse patterns on the neuronal firing by incorporating the mechanism of submyelin K+ accumulation and non-linear gating dynamics of Na+ channel activation.

Implications of the Neuronal Firing Induced by Different Patterns of Varying Inter-Pulse Interval

Pulse frequency is the most commonly used parameter to adjust the effects of DBS applications. For example, different temporal patterns of stimulation with varying frequency (IPI) have been investigated to modulate central nervous system arousal in intact mice (Quinkert et al., 2010; Quinkert and Pfaff, 2012) and in mice with traumatic brain injury (Tabansky et al., 2014) and to suppress epileptogenesis (Santos-Valencia et al., 2019). Irregular temporal patterns of stimulation clearly offer the opportunity to improve the efficacy and efficiency of DBS therapies (Grill, 2018). However, the vast design space for arbitrary sequences of varying IPI poses a clear challenge in designing irregular temporal patterns of stimulation. It is not clear how patterns of stimulation should be selected for therapeutic advantage in treating various diseases. In this study, we find that even rearranging the order of the same set of varying IPI can change the effect of stimulation pulses substantially. Notably, the IPIs need to vary only in a small range such as 5–10 ms (i.e., 200–100 Hz) with a mean pulse frequency of ∼130 Hz, a commonly used frequency in DBS therapy. In addition, the changes in the orders of varying IPI did not significantly change the total amount of neuronal firing during HFS (Table 3) but only redistributed the firing timing of individual neurons.

Furthermore, our investigation of the neuronal responses to specific IPI episodes in a random IPI sequence provides clues to designing novel temporal patterns of stimulation. The salient changes in the ability of a pulse sequence to recruit a population of neurons by the specific IPI episodes (several shorter IPIs followed by a longer IPI as in Figure 6C) may exert stronger effects on neurons than regular and mild impulses from stimulation with constant IPI. This suggests that adjusting the amount of such IPI episodes may quantitatively regulate the stimulation effect, which acts as a “dose” of stimulation. The recruitment of a large population of neurons may be beneficial for treating diseases caused by excessive suppression of neuronal activity, such as disorders of consciousness after traumatic brain injury. In another perspective, the enhanced effects by random IPI suggest a potential to lower the electrical energy consumption, because only a decreased strength of stimulation may be needed to obtain a similar level of effect as constant IPI. Therefore, the present study with random IPI may inform the design of novel temporal patterns of stimulation that provide therapeutic advantage, less energy consumption, and more safety over regular stimulation for extending DBS applications in treating more brain diseases.

Limitations of the Study

In the study, the effects of varying IPI were generated from the direct action of stimulation on the neuronal membrane, not involving specific synaptic transmission or neuronal network. Because the mechanisms of action involving the non-linear dynamics of ion channels are universal to the membrane of various neurons in brain, the effects of varying IPI may also apply to neurons in other brain regions directly activated by stimulations. However, synaptic and network involvements may complicate these effects and impose additional modulation to projection neurons in orthodromic propagations, which needs further investigations. Furthermore, even for the antidromic situation, in addition to the depolarization block of axons investigated here, sustained HFS might induce changes in other features along the axon to soma. Also, the stochastic network activity received by soma and the complexity of axonal structures might affect soma’s responses to the axonal HFS. Therefore, more accurate modeling with considerations of these factors is needed to further confirm the results of the present study.

In addition, the rat experiments were performed under anesthesia. The anesthetic urethane could reduce the synaptic responses, elevate the threshold of action potential firing, and decrease the firing rate of neurons (Mercer et al., 1978; Shirasaka and Wasterlain, 1995). Although the antidromic activation paradigm used in the present study did not involve synaptic transmissions, the excitability of soma may be changed by the anesthetic, thereby possibly affecting the responses of neurons to the stimulation. Therefore, further studies are needed to duplicate the results in awake animals.


The present study suggests that with the same amount of pulses together with the same distribution of IPI, only changing the order of varying IPI can introduce substantial changes in neuronal responses to HFS. Non-linear recovery course of membrane excitability might be the underlying cellular mechanisms resulting in the differences in neuronal responses. The finding provides a new strategy for developing novel stimulation paradigms to modulate the firing patterns of neuron populations, which may improve the DBS efficacy and extend the application of DBS therapy to more brain diseases.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Ethics Statement

The animal study was reviewed and approved by the Institutional Animal Care and Ethics Committee, Zhejiang University, Hangzhou, China.

Author Contributions

ZF and LZ designed the study. HH and ZW performed the animal experiments. LZ, HH, and YY analyzed the experimental data. LZ created the computational model, and XW contributed to the model modification. LZ performed the simulations. ZF, LZ, and XW interpreted the results and wrote the manuscript. All authors approved the final version for submission.


This work was supported by the National Natural Science Foundation of China (No. 30970753).

Conflict of Interest

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.


Battefeld, A., Tran, B. T., Gavrilis, J., Cooper, E. C., and Kole, M. H. (2014). Heteromeric Kv7.2/7.3 channels differentially regulate action potential initiation and conduction in neocortical myelinated axons. J. Neurosci. 34, 3719–3732. doi: 10.1523/JNEUROSCI.4206-13.2014

CrossRef Full Text | Google Scholar

Bellinger, S. C., Miyazawa, G., and Steinmetz, P. N. (2008). Submyelin potassium accumulation may functionally block subsets of local axons during deep brain stimulation: a modeling study. J. Neural Eng. 5, 263–274. doi: 10.1088/1741-2560/5/3/001

CrossRef Full Text | Google Scholar

Bergey, G. K. (2013). Neurostimulation in the treatment of epilepsy. Exp. Neurol. 244, 87–95. doi: 10.1016/j.expneurol.2013.04.004

CrossRef Full Text | Google Scholar

Birdno, M. J., Cooper, S. E., Rezai, A. R., and Grill, W. M. (2007). Pulse-to-pulse changes in the frequency of deep brain stimulation affect tremor and modeled neuronal activity. J. Neurophysiol. 98, 1675–1684. doi: 10.1152/jn.00547.2007

CrossRef Full Text | Google Scholar

Bittar, R. G., Burn, S. C., Bain, P. G., Owen, S. L., Joint, C., Shlugman, D., et al. (2005). Deep brain stimulation for movement disorders and pain. J. Clin. Neurosci. 12, 457–463. doi: 10.1016/j.jocn.2004.09.001

CrossRef Full Text | Google Scholar

Cagnan, H., Pedrosa, D., Little, S., Pogosyan, A., Cheeran, B., Aziz, T., et al. (2017). Stimulating at the right time: phase-specific deep brain stimulation. Brain 140, 132–145. doi: 10.1093/brain/aww286

CrossRef Full Text | Google Scholar

Cota, V. R., Medeiros, D. D. C., Doretto, M. C., and Moraes, M. F. D. (2009). Distinct patterns of electrical stimulation of the basolateral amygdala influence pentylenetetrazole seizure outcome. Epilepsy Behav. 14, 26–31. doi: 10.1016/j.yebeh.2008.09.006

CrossRef Full Text | Google Scholar

Cukiert, A., Cukiert, C. M., Argentoni-Baldochi, M., Baise, C., Forster, C. R., Mello, V. A., et al. (2011). Intraoperative neurophysiological responses in epileptic patients submitted to hippocampal and thalamic deep brain stimulation. Seizure 20, 748–753. doi: 10.1016/j.seizure.2011.07.003

CrossRef Full Text | Google Scholar

Cury, R. G., Fraix, V., Castrioto, A., Pérez Fernández, M. A., Krack, P., Chabardes, S., et al. (2017). Thalamic deep brain stimulation for tremor in Parkinson disease, essential tremor, and dystonia. Neurology 89, 1416–1423. doi: 10.1212/WNL.0000000000004295

CrossRef Full Text | Google Scholar

Dandekar, M. P., Fenoy, A. J., Carvalho, A. F., Soares, J. C., and Quevedo, J. (2018). Deep brain stimulation for treatment-resistant depression: an integrative review of preclinical and clinical findings and translational implications. Mol. Psychiatry 23, 1094–1112. doi: 10.1038/mp.2018.2

CrossRef Full Text | Google Scholar

Dayal, V., Limousin, P., and Foltynie, T. (2017). Subthalamic nucleus deep brain stimulation in parkinson’s disease: the effect of varying stimulation parameters. J. Parkinsons Dis. 7, 235–245. doi: 10.3233/JPD-171077

CrossRef Full Text | Google Scholar

Deeb, W., Salvato, B., Almeida, L., Foote, K. D., Amaral, R., Germann, J., et al. (2019). Fornix-region deep brain stimulation-induced memory flashbacks in Alzheimer’s Disease. N. Engl. J. Med. 381, 783–785. doi: 10.1056/NEJMc1905240

CrossRef Full Text | Google Scholar

Feng, Z., Ma, W., Wang, Z., Qiu, C., and Hu, H. (2019). Small changes in inter-pulse-intervals can cause synchronized neuronal firing during high-frequency stimulations in rat hippocampus. Front. Neurosci. 13:36. doi: 10.3389/fnins.2019.00036

CrossRef Full Text | Google Scholar

Feng, Z., Wang, Z., Guo, Z., Zhou, W., Cai, Z., and Durand, D. M. (2017). High frequency stimulation of afferent fibers generates asynchronous firing in the downstream neurons in hippocampus through partial block of axonal conduction. Brain Res. 1661, 67–78. doi: 10.1016/j.brainres.2017.02.008

CrossRef Full Text | Google Scholar

Feng, Z., Yu, Y., Guo, Z., Cao, J., and Durand, D. M. (2014). High frequency stimulation extends the refractory period and generates axonal block in the rat hippocampus. Brain Stimul. 7, 680–689. doi: 10.1016/j.brs.2014.03.011

CrossRef Full Text | Google Scholar

Feng, Z., Zheng, X., Yu, Y., and Durand, D. M. (2013). Functional disconnection of axonal fibers generated by high frequency stimulation in the hippocampal CA1 region in-vivo. Brain Res. 1509, 32–42. doi: 10.1016/j.brainres.2013.02.048

CrossRef Full Text | Google Scholar

Fisher, R. S., and Velasco, A. L. (2014). Electrical brain stimulation for epilepsy. Nat. Rev. Neurol. 10, 261–270. doi: 10.1038/nrneurol.2014.59

CrossRef Full Text | Google Scholar

Förstl, J., Galvan, M., and Ten, B. G. (1982). Extracellular K+ concentration during electrical stimulation of rat isolated sympathetic ganglia, vagus and optic nerves. Neuroscience 7, 3221–3229. doi: 10.1016/0306-4522(82)90244-5

CrossRef Full Text | Google Scholar

Grill, W. M. (2018). Temporal pattern of electrical stimulation is a new dimension of therapeutic innovation. Curr. Opin. Biomed. Eng. 8, 1–6. doi: 10.1016/j.cobme.2018.08.007

CrossRef Full Text | Google Scholar

Guo, Z., Feng, Z., Wang, Y., and Wei, X. (2018). Simulation study of intermittent axonal block and desynchronization effect induced by high-frequency stimulation of electrical pulses. Front. Neurosci. 12:858. doi: 10.3389/fnins.2018.00858

CrossRef Full Text | Google Scholar

Hines, M. L., and Carnevale, N. T. (1997). The NEURON simulation environment. Neural Comput. 9, 1179–1209. doi: 10.1162/neco.1997.9.6.1179

CrossRef Full Text | 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.1007/bf02459568

CrossRef Full Text | Google Scholar

Hu, W., Tian, C., Li, T., Yang, M., Hou, H., and Shu, Y. (2009). Distinct contributions of Na(v)1.6 and Na(v)1.2 in action potential initiation and backpropagation. Nat. Neurosci. 12, 996–1002. doi: 10.1038/nn.2359

CrossRef Full Text | Google Scholar

Jensen, A. L., and Durand, D. M. (2009). High frequency stimulation can block axonal conduction. Exp. Neurol. 220, 57–70. doi: 10.1016/j.brainres.2017.02.008

CrossRef Full Text | Google Scholar

Kloosterman, F., Peloquin, P., and Leung, L. S. (2001). Apical and basal orthodromic population spikes in hippocampal CA1 in vivo show different origins and patterns of propagation. J. Neurophysiol. 86, 2435–2444. doi: 10.1152/jn.2001.86.5.2435

CrossRef Full Text | Google Scholar

Kuncel, A. M., Birdno, M. J., Swan, B. D., and Grill, W. M. (2012). Tremor reduction and modeled neural activity during cycling thalamic deep brain stimulation. Clin. Neurophysiol. 123, 1044–1052. doi: 10.1016/j.clinph.2011.07.052

CrossRef Full Text | Google Scholar

Laxpati, N. G., Kasoff, W. S., and Gross, R. E. (2014). Deep brain stimulation for the treatment of epilepsy: circuits. Targets Trials Neurotherap. 11, 508–526. doi: 10.1007/s13311-014-0279-9

CrossRef Full Text | Google Scholar

Lorincz, A., and Nusser, Z. (2010). Molecular identity of dendritic voltage-gated sodium channels. Science 328, 906–909. doi: 10.1126/science.1187958

CrossRef Full Text | Google Scholar

Lozano, A. M., and Lipsman, N. (2013). Probing and regulating dysfunctional circuits using deep brain stimulation. Neuron 77, 406–424. doi: 10.1016/j.neuron.2013.01.020

CrossRef Full Text | 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

CrossRef Full Text | Google Scholar

Medeiros, D. C., Cota, V. R., Vilela, M. R. S. P., Mourão, F. A. G., Massensini, A. R., and Moraes, M. F. D. (2012). Anatomically dependent anticonvulsant properties of temporally-coded electrical stimulation. Epilepsy Behav. 23, 294–297. doi: 10.1016/j.yebeh.2012.01.004

CrossRef Full Text | Google Scholar

Mercer, L. F. Jr., Remley, N. R., and Gilman, D. P. (1978). Effects of urethane on hippocampal unit activity in the rat. Brain Res. Bull. 3, 567–570. doi: 10.1016/0361-9230(78)90089-8

CrossRef Full Text | Google Scholar

Mesquita, M. B. S., Medeiros, D. C., Cota, V. R., Richardson, M. P., Williams, S., and Moraes, M. F. D. (2011). Distinct temporal patterns of electrical stimulation influence neural recruitment during PTZ infusion: an fMRI study. Prog. Biophys. Mol. Biol. 105, 109–118. doi: 10.1016/j.pbiomolbio.2010.10.005

CrossRef Full Text | Google Scholar

Migliore, M., Hoffman, D. A., Magee, J. C., and Johnston, D. (1999). Role of an A-Type K + conductance in the back-propagation of action potentials in the dendrites of hippocampal pyramidal neurons. J. Comput. Neurosci. 7, 5–15. doi: 10.1023/a:1008906225285

CrossRef Full Text | Google Scholar

Nowak, L. G., and Bullier, J. (1998). Axons, but not cell bodies, are activated by electrical stimulation in cortical gray matter I. Evidence from chronaxie measurements. Exp. Brain Res. 118, 477–488. doi: 10.1007/s002210050304

CrossRef Full Text | Google Scholar

Poolos, N. P., Mauk, M. D., and Kocsis, J. D. (1987). Activity-evoked increases in extracellular potassium modulate presynaptic excitability in the CA1 region of the hippocampus. J. Neurophysiol. 58, 404–416. doi: 10.1152/jn.1987.58.2.404

CrossRef Full Text | Google Scholar

Popovych, O. V., and Tass, P. A. (2014). Control of abnormal synchronization in neurological disorders. Front. Neurol. 5:268. doi: 10.3389/fneur.2014.00268

CrossRef Full Text | Google Scholar

Quinkert, A. W., and Pfaff, D. W. (2012). Temporal patterns of deep brain stimulation generated with a true random number generator and the logistic equation: effects on CNS arousal in mice. Behav. Brain Res. 229, 349–358. doi: 10.1016/j.bbr.2012.01.025

CrossRef Full Text | Google Scholar

Quinkert, A. W., Schiff, N. D., and Pfaff, D. W. (2010). Temporal patterning of pulses during deep brain stimulation affects central nervous system arousal. Behav. Brain Res. 214, 377–385. doi: 10.1016/j.bbr.2010.06.009

CrossRef Full Text | Google Scholar

Rizzone, M., Lanotte, M., Bergamasco, B., Tavella, A., Torre, E., Faccani, G., et al. (2001). Deep brain stimulation of the subthalamic nucleus in Parkinson’s disease: effects of variation in stimulation parameters. J. Neurol. Neurosurg. Psychiatry 71, 215–219. doi: 10.1136/jnnp.71.2.215

CrossRef Full Text | Google Scholar

Röper, J., and Schwarz, J. R. (1989). Heterogeneous distribution of fast and slow potassium channels in myelinated rat nerve fibres. J. Physiol. 416, 93–110. doi: 10.1113/jphysiol.1989.sp017751

CrossRef Full Text | Google Scholar

Rosenbaum, R., Zimnik, A., Zheng, F., Turner, R. S., Alzheimer, C., Doiron, B., et al. (2014). Axonal and synaptic failure suppress the transfer of firing rate oscillations, synchrony and information during high frequency deep brain stimulation. Neurobiol. Dis. 62, 86–99. doi: 10.1016/j.nbd.2013.09.006

CrossRef Full Text | Google Scholar

Santos-Valencia, F., Almazan-Alvarado, S., Rubio-Luviano, A., Valdes-Cruz, A., Magdaleno-Madrigal, V. M., and Martinez-Vargas, D. (2019). Temporally irregular electrical stimulation to the epileptogenic focus delays epileptogenesis in rats. Brain Stimul. doi: 10.1016/j.brs.2019.07.016

CrossRef Full Text | Google Scholar

Shirasaka, Y., and Wasterlain, C. G. (1995). The effect of urethane anesthesia on evoked potentials in dentate gyrus. Eur. J. Pharmacol. 282, 11–17. doi: 10.1016/0014-2999(95)00244-f

CrossRef Full Text | Google Scholar

Swan, B. D., Brocker, D. T., Hilliard, J. D., Tatter, S. B., Gross, R. E., Turner, D. A., et al. (2016). Short pauses in thalamic deep brain stimulation promote tremor and neuronal bursting. Clin. Neurophysiol. 127, 1551–1559. doi: 10.1016/j.clinph.2015.07.034

CrossRef Full Text | Google Scholar

Tabansky, I., Quinkert, A. W., Rahman, N., Muller, S. Z., Lofgren, J., Rudling, J., et al. (2014). Temporally-patterned deep brain stimulation in a mouse model of multiple traumatic brain injury. Behav. Brain Res. 273, 123–132. doi: 10.1016/j.bbr.2014.07.026

CrossRef Full Text | Google Scholar

Varona, P., Ibarz, J. M., López-Aguado, L., and Herreras, O. (2000). Macroscopic and subcellular factors shaping population spikes. J. Neurophysiol. 83, 2192–2208. doi: 10.1152/jn.2000.83.4.2192

CrossRef Full Text | Google Scholar

Wang, Z., Feng, Z., and Wei, X. (2018). Axonal stimulations with a higher frequency generate more randomness in neuronal firing rather than increase firing rates in rat hippocampus. Front. Neurosci. 12:783. doi: 10.3389/fnins.2018.00783

CrossRef Full Text | Google Scholar

Wichmann, T., and DeLong, M. R. (2016). Deep brain stimulation for movement disorders of basal ganglia origin: restoring function or functionality? Neurotherapeutics 13, 264–283. doi: 10.1007/s13311-016-0426-6

CrossRef Full Text | Google Scholar

Zheng, F., Lammert, K., Nixdorf-Bergweiler, B. E., Steigerwald, F., Volkmann, J., and Alzheimer, C. (2011). Axonal failure during high frequency stimulation of rat subthalamic nucleus. J. Physiol. 589(Pt 11), 2781–2793. doi: 10.1113/jphysiol.2011.205807

CrossRef Full Text | Google Scholar

Keywords: high-frequency stimulation, inter-pulse interval, in vivo experiment, computational modeling, hippocampal pyramidal cell, membrane dynamics, non-linearity, population spike

Citation: Zheng L, Feng Z, Hu H, Wang Z, Yuan Y and Wei X (2020) The Appearance Order of Varying Intervals Introduces Extra Modulation Effects on Neuronal Firing Through Non-linear Dynamics of Sodium Channels During High-Frequency Stimulations. Front. Neurosci. 14:397. doi: 10.3389/fnins.2020.00397

Received: 29 August 2019; Accepted: 30 March 2020;
Published: 19 May 2020.

Edited by:

Caleb Kemere, Rice University, United States

Reviewed by:

Martin Han, University of Connecticut, United States
Daniel Llewellyn Rathbun, Henry Ford Health System, United States

Copyright © 2020 Zheng, Feng, Hu, Wang, Yuan and Wei. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Zhouyan Feng,;