Dynamics of neuronal firing modulated by high-frequency electrical pulse stimulations at axons in rat hippocampus

Objective. The development of electrical pulse stimulations in brain, including deep brain stimulation, is promising for treating various brain diseases. However, the mechanisms of brain stimulations are not yet fully understood. Previous studies have shown that the commonly used high-frequency stimulation (HFS) can increase the firing of neurons and modulate the pattern of neuronal firing. Because the generation of neuronal firing in brain is a nonlinear process, investigating the characteristics of nonlinear dynamics induced by HFS could be helpful to reveal more mechanisms of brain stimulations. The aim of present study is to investigate the fractal properties in the neuronal firing generated by HFS. Approach. HFS pulse sequences with a constant frequency 100 Hz were applied in the afferent fiber tracts of rat hippocampal CA1 region. Unit spikes of both the pyramidal cells and the interneurons in the downstream area of stimulations were recorded. Two fractal indexes—the Fano factor and Hurst exponent were calculated to evaluate the changes of long-range temporal correlations (LRTCs), a typical characteristic of fractal process, in spike sequences of neuronal firing. Main results. Neuronal firing at both baseline and during HFS exhibited LRTCs over multiple time scales. In addition, the LRTCs significantly increased during HFS, which was confirmed by simulation data of both randomly shuffled sequences and surrogate sequences. Conclusion. The purely periodic stimulation of HFS pulses, a non-fractal process without LRTCs, can increase rather than decrease the LRTCs in neuronal firing. Significance. The finding provides new nonlinear mechanisms of brain stimulation and suggests that LRTCs could be a new biomarker to evaluate the nonlinear effects of HFS.


Introduction
Neuronal firing of action potentials, also known as spikes in an extracellular recording, exhibits distinct properties of nonlinear dynamics such as fractal, bifurcation, and chaos [1][2][3].Fractals represent a nested pattern, where the structures of an object at smaller scales resemble those at larger scales [4].In the time domain, fractal patterns can be characterized by long-range temporal correlations (LRTCs) among events that extend over multiple time scales, such as the fractal patterns in firing sequences of spikes, EEG and local field potentials (LFP) [5][6][7][8].For instance, in a spike sequence with LRTCs, the fractal patterns can be observed in the fluctuations of the inter-spikeinterval (ISI) or the spike rate.Smaller fluctuations at shorter time scales are nested within larger fluctuations at longer time scales, resulting in fluctuations proportional to time scales [6,9].
LRTCs in spontaneous neuronal firing have been observed in various brain regions, including the basal ganglia, hippocampus and visual cortex [1,10].
These LRTCs indicate history-effects or memory in the neuronal firing, implying that two spikes separated far away in a firing sequence can be dependent on each other [7].LRTCs are important for neuronal activity related to information processing.Fractal patterns with LRTCs in the neuronal firing have been shown to reflect a balance between stability and excitability and to represent an optimal state for information processing in the brain [8,11].With a similar number of action potentials in a given period, fractal firing can encode more information than periodic or random firing.An increase of LRTCs in neuronal firing may indicate an enhanced capacity for information processing in the brain [5,11].A decrease of LRTCs means a decrease in the efficiency of information processing in neural networks.Previous experimental and clinical studies have shown decreases of LRTCs in the activity of abnormal neural networks with brain diseases such as epilepsy, Alzheimer's disease, schizophrenia, major depression, and so on [8,12,13].Therefore, recovery of LRTCs may serve as a biomarker that reflects a therapeutic effect in improving symptoms of the diseases [8,14].
Indeed, previous studies have shown correlations between increases of LRTCs and improvements of brain situations generated by medication or external stimulations.For instance, an increase of LRTCs can be a new biomarker of the medication therapeutic effects in patients with Parkinson's disease or patients with infantile spasms [8,14].An acupuncture stimulation applied on the bilateral hind limbs, which is effective to treat pain, has been shown to increase the LRTCs in firing sequences of individual neurons in rat spinal dorsal horn.It suggests a correlation between the LRTCs changes and the curative effects [15].Visual flash stimulations have been shown to improve memory and perceptual performance in human beings.The improvement is accompanied by an increase in the LRTCs of EEG during resting states with eyes closed [16].In addition, applying stimulations of electrical pulses directly to the nervous system can also change the LRTCs of neuronal firing [17][18][19].For instance, applying sustained electrical pulses with a high-frequency of 70 Hz to spinal cords can increase the LRTCs of EEG in patients under minimally conscious state, which is accompanied by a restoration of cortical information integration [20].In contrast, applying low-frequency (0.33 Hz) pulse stimulations on the left median nerve can decrease the LRTCs of EEG in the sensorimotor and occipito-parietal region of human brain, which is accompanied by a hindrance for neuronal networks to retain memory of neuronal activity [21].The results of these studies suggest that the change in LRTCs can serve as a potential biomarker for nervous system therapies.
In addition to the above mentioned stimulations, electrical pulse stimulation in the brain has shown effectiveness in treating various brain disorders with different targets of stimulation, such as subthalamic nucleus and globus pallidus interna stimulation for Parkinson's disease and essential tremor, hippocampal stimulation for epilepsy and Alzheimer's disease, thalamus or cingulum for pain, as well as the nucleus accumbens for major depression and obsessive compulsive disorder, and so on [22,23].In general, all of these brain stimulations can be covered by deep brain stimulation (DBS).DBS typically utilizes high-frequency stimulations (HFS) with a constant pulse frequency in a range of 90-200 Hz [22,23].Studies have shown that HFS can increase the neuronal firing downstream of stimulation sites through axonal projections and modulate the neuronal firing patterns [24][25][26][27][28].It has been reported that HFS can also increase the LRTCs in LFP oscillations that represent integrated activity of a population of neurons [29].However, it is unclear whether the LRTCs in firing of individual neurons could be altered by the HFS.
The LRTCs represent nonlinear and nonstationary events, indicating a high degree of variability in neuronal firing, while HFS with purely periodical pulses of a constant inter-pulse-interval (IPI) represents regular stimulations.Additionally, previous studies have shown that HFS can modulate the neurons to fire at a specific phase in the IPI to generate time-locked neuronal firing [25,28].Furthermore, the therapeutic effects of brain HFS mainly occur during the period of sustained stimulation, not after HFS.The effects of HFS are reversible after the end of HFS.Along this line, the HFS could increase the regularity and reduce the LRTCs in neuronal firing.On the other hand, however, clinical studies have shown that effective medication treatments for patients with Parkinson's disease, epilepsy and minimally conscious state are associated with an increase of LRTCs in EEG or in LFP [8,14,20].Therefore, an interesting question arises: could the LRTCs of neuronal firing also increase during the periodic pulse stimulations of HFS?
To answer the question, we analyzed the LRTCs in the spike sequences of neuronal firing directly modulated by HFS in the hippocampal CA1 region in anesthetized rats in-vivo.The HFS was applied in the afferent fiber tracts (i.e.Schaffer collaterals) of CA1 region.Unit spikes of both the pyramidal cells and the interneurons in the downstream area of stimulations were recorded.To investigate the changes of LRTCs of the neuronal firing during HFS, the Fano factor and Hurst exponent [1,10] of the spike sequences were calculated.Results of the study at a cell level can reveal new nonlinear dynamics of real firing of individual neurons in brain in-vivo and provide information to reveal more mechanisms of HFS effects.

Animals and surgery
The animal experiment was approved by the Laboratory Animal Welfare and Ethics Committee of Zhejiang University (ZJU20210108).The data were recorded from 15 male adult Sprague-Dawley rats (325 ± 27 g).Details of the surgery have been reported previously [27,30].Briefly, the rats were anesthetized with urethane (1.25 g kg −1 , i.p.), and fixed in a stereotaxic frame (Stoelting Co., USA).A 16channel silicon electrode probe (#Poly2, NeuroNexus Technologies, USA) was positioned in left hippocampal CA1 region (AP −3.5 mm; ML 2.7 mm; DV ∼2.5 mm) for extracellular recording of neuronal signals (figures 1(a) and (b)).A concentric bipolar stainless-steel electrode (#CBCSG75, FHC Inc., USA) was placed in the afferent axons (i.e. the Schaffer collaterals) of the CA1 region (AP −2.2 mm; ML 2.2 mm; DV ∼2.8 mm) to orthodromically activate the CA1 neurons downstream (figure 1(a)).According to the unit spike signals and the waveforms of pulse-evoked potentials that appeared serially in the recording array, the correct sites of the two electrodes can be judged [31].

Stimulation and recording
Biphasic current pulses with a width of 100 µs per phase were utilized for stimulation.The pulses were generated by an isolated pulse stimulator (Model 3800, A-M Systems Inc., USA).The pulse intensity was 0.3-0.4mA that was able to induce a population spikes (PS) with an amplitude of approximately 3/4 maximal value according to an inputoutput curve of single pulse stimulation with a gradually increased intensity.The PS is an integrated waveform of action potentials synchronously fired by a population of neurons.The duration of HFS pulse trains was 2 min with a pulse frequency of 100 Hz.
The extracellular electrical signals collected by the recording electrode array were amplified 100-fold by an extracellular amplifier (Model 3600, A-M Systems Inc., USA) with a band-pass filter in a range of 0.3-5000 Hz.Then, the signals were sampled at 20 kHz using a PowerLab data acquisition system (Model PL3516, AD Instruments Inc., Australia) and were stored for offline analysis.

Data analysis for unit spikes
Signals of four neighboring channels located in the pyramidal layer of CA1 region were used to analyze unit spikes (figure 1(b)).First, the HFS stimulation artifacts were removed by an algorithm implemented by a custom MATLAB program: The artifacts were detected automatically by setting a slope threshold based on the fact that the slope of artifacts was much greater than that of neuronal signals.Subsequently, the segment of ∼0.9 ms signal including each artifact was replaced by a straight-line segment connecting the two end points of artifact segment [32].Then, the artifact-free signals were filtered (>500 Hz) to obtain multiple unit activity (MUA) by removing the LFP.The unit spikes in the MUA signals were detected by a negative amplitude threshold at 5-fold standard deviations of the MUA signals.Finally, the single unit activity (SUA) of pyramidal cells and interneurons were obtained by a spike sorting process [27]: feature vectors of the spike waveforms (including the first principal components and amplitudes) were calculated from the MUA signals of the four neighboring channels.These feature vectors were then utilized for spike sorting using the open-source software SpikeSort 3D developed by Neuralynx Inc. (www.Neuralynx.com).The sorting results were verified by examining the histograms of ISIs of the sorted spike sequences.Only the sequences with ISI histograms exhibiting clear refractory periods were considered as single unit spikes.Based on the distinctive waveform features of spikes from interneurons and pyramidal cells, unit spikes from the two types of individual neurons were distinguished [27].Briefly, unit spikes with a mean rising phase shorter than 0.4 ms were classified as interneurons, while those longer than 0.7 ms were classified as pyramidal cells [27].
To investigate the effect of HFS on the neuronal firing in the downstream region, the mean firing rates of SUA and the ISI of spike sequences were calculated during the periods of late 90 s of 2 min HFS and in 90 s baseline recording before HFS as a control.The initial 30 s periods of 2 min HFS were not used to analyze unit spikes due to the interference of the stimulationevoked PS [27,30].The firing rhythms of both types of neurons were examined using the following spectrum of auto-correlograms for binary sequences representing the spike firing of individual neurons [26,33,34].Firstly, the auto-correlogram of each binary sequence was calculated and normalized by its mean value to eliminate the impact of firing rate.And, the mean component of auto-correlogram was removed.Then, the power spectrum density (PSD) of the auto-correlogram of individual neurons was calculated using Welch's periodogram with a Hanning window and a window overlap of 50%.The frequency resolution of PSD was 0.98 Hz.This PSD for unit spikes was also termed as unit-PSD.
To evaluate the changes of LRTCs (one of fractal properties) in spike sequences of individual neurons, two indexes of LRTCs, namely the Fano factor and Hurst exponent, were calculated [1,10].The Fano factor, abbreviated as F(T), is defined as the ratio of the variance of spike numbers to the mean of spike numbers in non-overlapping time windows of size T (in seconds): For the recordings of late 90 s HFS and 90 s baseline, T is increased from a minimum of 0.01 s to a maximum of 15 s.N i (T) is the spike number in the ith window of size T. The Fano factor curve F(T) is plotted with double-logarithmic scales.
For a spike sequence with LRTCs, i.e, a fractal process, the F(T) value increases as a power-law function of T. This type of F(T) curve represents a greater variance of spike numbers in a longer window size T [6,35], indicating a nested pattern of clusters (one of fractal patterns) in spike sequences.That is, a cluster in a longer window consists of sub-clusters in shorter windows.In addition, the distributions of the clusters and the sub-clusters in a fractal process are nonuniform.The F(T) increases with a longer T because a longer window is apt to contain more clusters with various patterns [36].In time domain, this characteristic of fractal process makes the values of F(T) situate around a linear regression line with a slope α in the double-logarithmic plot of F(T), resulting a power-law function F(T) ∝ T α [35,36].Therefore, the slope α of F(T) provides a quantitative estimation of fractal scaling [1].In the present study, the α in the range with T > 10 0 s (i.e.T > 1 s) was calculated to estimate the fractal scaling of spike sequences.If a spike sequence is a Poisson process without LRTCs but with fluctuations in its spike instantaneous rates, its F(T) would be approximately 1.0 for all T because the variance of a Poisson process is equal to its mean.If a spike sequence is periodic with no changes and no LRTCs, i.e. a nonfractal process, its F(T) would approach 0 when T increases to enough long to contain enough number of cycles and to make the variance close to 0. The α of the two types of sequences is approximately 0.
The second index of LRTCs, Hurst exponent (H), is defined on the ISIs of spike sequences by the following rescaled range method [10,37].First, the ISI sequence of spikes, with a total ISI number of N, is divided into M segments, denoted as I m , each with a length of d.Here, d × M = N, and m = 1,..., M. Next, for each segment I m , calculate the mean (E m ), the standard deviation (S m ) and the cumulative difference (X k,m ) of ISIs.The X k,m is the running sum of the differences between individual ISI and the mean E m : For each segment I m , the rescaled range R m is defined as the difference between the maximum and minimum values of the X k,m normalized by S m : Thus, the mean rescaled range R d of all segments I m with a length d is: The length d is increased from a minimum of 4 (10 0.6 ) to a maximum of N/3.Finally, the Hurst exponent (H) is the slope of the linear regression line of R d against d in a double-logarithmic plot.The H value ranges from 0 to 1.0 and determines whether the spike sequence is fractal.As long as H ̸ = 0.5, the spike sequence exhibits long-range correlated [38,39].If 0 < H < 0.5, the LRTCs of spike sequence is negative: a decrease of ISI is generally followed by an increase, and an increase is followed by a decrease.Conversely, if 0.5 < H < 1, the LRTCs of spike sequence is positive: a decrease of ISI is followed by further decreases, and an increase is followed by more increases.If H = 0.5, successive changes in the ISIs are independent and do not exhibit LRTCs.Therefore, H can indicate the directions of LRTCs, while the α of F(T) cannot.All the data are presented as the mean ± SD, with 'n' representing the number of neurons.The statistical significance of the differences between the groups during HFS and at the baseline recording was tested using a paired t-test.

Simulation data
The LRTCs of neuronal firing during HFS were also confirmed by two types of simulation spike sequences: randomly shuffled sequences and surrogate sequences.A randomly shuffled sequence was created by rearranging the ISIs of the original spike sequence obtained from experimental recording (figure 2(a)).To do this, a random number was assigned to each of the original ISIs, and then the ISIs were rearranged in ascending order based on their assigned numbers (figure 2(b)).Consequently, the mean ISI, ISI variance, and ISI histogram of the randomly shuffled sequence were identical to those of the original spike sequence.However, the temporal characteristics of the original spike sequence, such as the order of ISIs and the interdependence between ISIs, were completely disrupted [6,36].
A surrogate sequence was created by replacing the ISIs in the original spike sequence with random numbers following a uniform distribution.The mean ISI of the surrogate sequence remained equal to that of the original spike sequence.The order of the random numbers followed the numerical value order of the original ISIs (figure 2(c)).As a result, the surrogate sequence preserved similar temporal patterns of the original spike sequence, but with different ISI values and an altered ISI histogram.
If the LRTCs of the original spike sequence is eliminated by its randomly shuffled sequence, but is reproduced by its surrogate sequence, it means that the LRTCs is generated by the temporal patterns (i.e. the ISI order) in the original spike sequence [6,36].

Increase of neuronal firing in the CA1 region during HFS
To investigate the fractal properties of neuronal firing during HFS, 2 min 100 Hz HFS sequences were applied at the Schaffer collateral of the rat hippocampal CA1 region.Similar to previous reports [30], large PS appeared at the initial period of the stimulations (figure 3(a)).Once the PS events disappeared, the stimulations continued to increase the firing of postsynaptic neurons in the CA1 region.In 15 experiment preparations, we recorded spike sequences from a total of 49 neurons in the CA1 pyramidal layer, including 31 pyramidal cells and 18 interneurons.During the late 90s periods of 2 min HFS, the mean firing rates of both pyramidal cells and interneurons were significantly higher than the baseline firing rate.However, the firing rates were still considerably smaller than the stimulation pulse frequency of 100 Hz (see the right plots in figures 3(b) and (c)).Furthermore, contrary to the stable firing rates

Change of LRTCs in the neuronal firing induced by HFS
The distributions of ISIs of spike sequences in CA1 neuronal firing were altered by the HFS (figure 4).During the baseline recording before HFS, most ISIs of the pyramidal cells fell within the range of 0-10 ms, resulting in a distinct peak around 5 ms in the average ISI histogram due to the bursty nature of these cells (left plot in figure 4(a) and blue curve in figure 4(b)).At baseline recording, interneurons had most of their ISIs within the range of 0-30 ms, with a small peak around 6 ms in the average ISI histograms (left plot in figure 4(c) and blue curve in figure 4(d)).However, during HFS, both types of neurons exhibited ISIs that clustered around multiples of the 10 ms IPI of the 100 Hz HFS, resulting in repeated peaks in the ISI histogram.The auto-correlograms and unit-PSD further showed that HFS altered neuronal firing rhythms to the stimulation frequency (figures 4(e) and (f)).During the period of HFS, the original rhythm in the firing of both pyramidal cells and interneurons were replaced by a new rhythm at the stimulation frequency, indicating by the peak at 100 Hz in the spectrums of unit-PSD (the plot at the lower right corner in figures 4(e) and (f)).These results indicated stimulation modulations on the firing patterns of the neurons.Next, we investigated whether the modulations could change the LRTCs of neuronal firing.
Although the stimulation of HFS was strictly periodic and the neuronal firing was modulated at a 100 Hz rhythm, the firing was non-stationary during HFS.The curves of firing rate per second for both types of neurons exhibited fractal structures with similar patterns in different temporal scales (figures 5(a) and (b)).We next calculated the curves of Fano-factor to quantify the fractal structures in the firing changes.
During the late 90s period of 2 min HFS, the F(T) of pyramidal cells' firing increased with the increase of window size T and reached a maximum value of approximately 10 (figure 6(a)).Importantly, a powerlaw relationship started from T ≈ 10 0 s (∼1.0 s) with a large positive slope α ≈ 0.7 in the doublelogarithmic plot.During baseline recording before HFS, the F(T) of pyramidal cells also increased with the increase of T and exhibited a power-law relationship starting from T ≈1 s, but with a smaller slope α = 0.17.Both the increase speed and the maximum value of F(T) were greater during HFS than during baseline.The F(T) slope (i.e.scaling exponent) during HFS (α = 0.72 ± 0.17) was significantly greater than α = 0.24 ± 0.15 at baseline (figure 6(b), P < 0.01, paired t-test, n = 31 pyramidal cells).
Similarly, in the window size of 10 0 s (1.0 s) to 10 1.18 s (15 s), the F(T) of interneurons also exhibited a power-law relationship in double-logarithmic plots both at baseline and during HFS (figure 6(c)).The mean F(T) slope during HFS (α = 0.83 ± 0.15, close to the theoretical maximum 1.0 for the slope) was significantly greater than the baseline value of α = 0.39 ± 0.17 (figure 6(d), P < 0.01, paired t-test, n = 18 interneurons).The power-law relationship in the F(T) and the large scaling exponents (α) during

Confirming the HFS-induced LRTCs by simulation sequences
One type of the simulation sequences (the randomly shuffled sequence) had the identical ISI values and distribution as the experimental spike sequence (the red curve in figures 8(a) and (b)).However, the temporal patterns in ISIs was completely disrupted due to the randomly rearrangement of original ISIs (figure 2(b)).For the pyramidal cells, both the autocorrelograms and unit-PSD of randomly shuffled sequence were similar to experimental spike sequence with a rhythm at the stimulation frequency 100 Hz (top and middle rows in figure 8(c)).However, in contrast to the rapid increase after T > 1 s in the F(T) of experimental spike sequences (the red curve in figure 8(e)), the F(T) of the randomly shuffled sequences hovered around a level line without substantial increase (the gray curve in figure 8(e)).For the firing of interneurons, the randomly shuffled sequence also showed a stimulation rhythm (100 Hz) in its auto-correlogram and unit-PSD, although the 100 Hz rhythm was decreased compare to the experimental spike sequence (top and middle rows in figure 8(d)).However, the F(T) of randomly shuffled sequences (the gray curve in figure 8(f)) fluctuated dramatically after T > 1 s with values much smaller than the F(T) of experiment spike sequences (the red curve in figure 8(f)).The F(T) slope α of randomly shuffled sequences of both pyramidal cell and interneuron approached zero, indicating no LRTCs in this type of simulation sequences [6,35].
The second type of simulation sequences (the surrogate sequences) preserved the temporal patterns in the original experiment spike sequences, but with a different ISI distribution (green curves in figures 8(a) and (b)).The surrogate sequences eliminated the rhythm at the stimulation frequency (bottom rows in figures 8(c) and (d)).However, the F(T) curves of the surrogate sequences for both pyramidal cells and interneurons showed a power-law relationship in the double-logarithmic plots (green curves in figures 8(e) and (f)).The mean F(T) slope α of the surrogate sequences for pyramidal cells (0.66 ± 0.12) and interneurons (0.81 ± 0.17) were respectively similar to the values of the original experiment spike sequences (right plots in figures 8(e) and (f), P = 0.25, paired t-test for n = 31 pyramidal cells; P = 0.55, paired t-test for n = 18 interneurons).
The eliminated power-law relationship in F(T) of randomly shuffled sequences as well as the preserved power-law relationship in F(T) of surrogate sequences further confirmed the LRTCs in the experiment spike sequences induced by HFS.Additionally, the increase of LRTCs by HFS was generated by the temporal structure (i.e. the ISI order) in neuronal firing rather than the ISI values and the distribution of ISIs, or the rhythm in spike sequences.

Discussion
The major finding of this study is that the activations of periodic axonal HFS without LRTCs can increase rather than decrease the LRTCs in the firing of pyramidal cells and interneurons in the postsynaptic area downstream of the stimulation site in the rat hippocampus.The increase of LRTCs was demonstrated by: (1) the power-law relationship and the slope increase in both F(T) and Hurst curves of the firing during HFS (figures 6 and 7), (2) the elimination of power-law relationship in F(T) curves

Possible mechanisms underlying the fractal neuronal firing induced by HFS
A single electrical pulse with a suprathreshold intensity can activate the axonal fibers to fire action potentials simultaneously, which then propagate to the terminals of axons and activate the downstream neurons through synaptic transmissions [27,30].However, sustained HFS of pulses with a constant frequency over 50 Hz can generate intermittent depolarization block in the axonal membranes [30, 40,41].Under this situation, the axons cannot follow each HFS pulse to fire due to an extension of the refractory period of axons caused by the depolarization block [27,40].In addition, HFS can also induce synaptic failures and neurotransmitter depletions [32,42].These effects prevent the postsynaptic neurons to follow each pulse of HFS to firing an action potential [25,27,28].Therefore, the firing rates of both pyramidal cells and interneurons during HFS are significantly lower than the 100 Hz pulses frequency of HFS but greater than those of baseline (figures 3(b) and (c)).
Intriguingly, the HFS-induced firing of downstream neurons was fractal rather than random or regular.This firing pattern can be caused by the intrinsic properties of neurons and the mechanism of HFS-induced axonal block and synaptic failure.Previous studies have shown that the LRTCs of neuronal firing biophysically originate from the intrinsic properties of ionic channels (e.g.Na + and K + channels) on neuronal membranes [10,36].The activation-induced opening and closing of ionic channels have fractal properties with LRTCs [43,44].Under the situation of axonal HFS, the intensive activations from pulses can increase the concentration of potassium ions [K + ] o in the peri-axonal spaces thereby elevating the membrane potential to a certain depolarization level to result in intermittent block [42,45].However, at the elevated potential level, the behaviors of Na + and K + channels still follow the intrinsic properties of nonlinear dynamics and initiate action potentials at an even more nonlinear level [46].Therefore, the HFS can increase the LRTCs in the induced firing (figures 5-7).In addition, fractal properties have been also observed in the secretion of neurotransmitters from the presynaptic membrane [47].The fractal properties in synaptic transmissions can further contribute to the generation of fractal firing in the post-synaptic neurons during HFS.Furthermore, the increase of firing rates (figures 3(b) and (c)) may also contribute to the increase of LRTCs during HFS.
In addition, interactions in the neuronal networks can also contribute to the HFS-induced fractal firing according to previous reports [6].For example, fractal firing can propagate from the sensory areas of cortex to the basal ganglia and induce fractal firing in substantia nigra cells [10,36,48].In the hippocampal CA1 region of the present study, pyramidal cells and interneurons are heavily interconnected in local inhibitory circuits [49,50].There are also interconnections among pyramidal cells and among interneurons [50].The feedback of inhibitory and/or excitatory interactions can result in the dependences of subsequent neuronal firing on the preceding firing thereby increasing the LRTCs in the firing sequence [51].
Although the spikes in firing sequences were timelocked with the HFS pulses thereby resulting in the peaks in n-folds of IPI (n × 10 ms) in the ISI histograms (figures 4(b) and (d)) and resulting in the rhythm at the stimulation frequency in the autocorrelograms and unit-PSD (figures 4(e) and (f)).This did not necessarily mean regular firing.The curves of firing rate showed fluctuations in different time scales (figure 5).And, the fluctuations were greater during HFS than at baseline (figures 3(b) and (c)).The results of simulation sequences further confirmed the non-uniform fractal property in the HFSinduced firing (figure 8).The fractal property of neuronal firing is generated by the ISI order (temporal structure) in the spike sequence rather than the distribution of ISIs or the rhythm in the spike sequence.The fractal property indicates nested clusters in the spike sequences, i.e. clusters in longer time scales consisting of subclusters in shorter time scales [1,48].
In summary, the increase of fractal properties in neuronal firing during HFS can be related with multiple factors, including the nonlinear dynamics of ionic channels in neuronal membranes with an elevated membrane potential, the fractal properties in synaptic transmissions, the increased firing rates, as well as the feedback in interactions of neuronal networks.

Implications of the study
First, our study provides new evidence of nonlinear mechanisms in the neuronal responses to electrical pulse stimulations of HFS in brain.These mechanisms are far from well understood.We found that the periodic stimulation of HFS pulses can increase rather than decrease the LRTCs in the firing of individual neurons directly under stimulations.The LRTCs in cell level can be the origin of LRTCs in the neuronal signals related with a population of neurons in networks or in local brain regions, such as LFP and EEG [15,16].The increases of LRTCs in LFP and EEG in patient brains under DBS, spinal cord stimulation as well as closed-loop neurofeedback stimulation have been reported to correlate with therapy effectiveness [8,14,16,20,29].For example, an effective therapy of thalamic DBS can increase the LRTCs in the high beta-band (21-30 Hz) of EEG recorded in patients with essential tremor [29].Therefore, an increase of LRTCs in the neuronal firing may be an important mechanism of brain stimulations.
Second, the increase of LRTCs in the HFSinduced firing suggests an improvement of information capacity and procession of neurons.HFSinduced axonal block may prevent transmissions of some signals, including pathological signals [28,30,52].However, the increase of LRTCs indicates that HFS can facilitate the processes of certain neuronal information which may be due to the increase of neuronal excitability by HFS.The LRTCs of neuronal firing mean a 'memory' or a history in the neuronal firing.That is, the firing timing of a specific spike is determined not only by the timing of preceding spikes, but also by the spikes that have occurred in 'distant' past [7,53].The historical effect results in fluctuations in neuronal firing over multiple time scales.Such fluctuations may enable neurons to maximize their capacity of information and may enhance the error tolerance of nervous system [10,54].Therefore, the increase of LRTCs by HFS may indicate an advantage in neuronal information processing.It provides new evidence to support the point that the effects of HFS in brain are different from destruction of brain tissue by an ablation.
Third, the study suggests that LRTCs can be a potential biomarker to assess the nonlinear effects of HFS that cannot be assessed using linear indexes, such as firing rates and ISIs [10].Since many effective treatments for brain diseases are associated with the restoration or increase of LRTCs in neuronal activity [8,14,20,29], LRTCs may be used to assess the outcome of brain therapies, including various electrical stimulations in brain.However, further investigations are needed to gather enough evidence to support the implementation of this application.
Taken together, our present study reveals a fractal mechanism at the level of individual neurons, provides a clue to support that HFS can facilitate neuronal information processing and suggests a potential to utilize the indexes of LRTCs to assess neuromodulation therapies.

Limitations of the study
First, we investigated the fractal neuronal firing induced by axonal stimulations in the hippocampal region because the region is one of the potential targets of DBS for treating brain disorders such as refractory epilepsy and Alzheimer's disease, and also because the lamellar structure of hippocampal CA1 region facilitates accurate locations of stimulating and recording sites [55,56].However, the axonal HFS-induced fractal neuronal firing found here needs to be duplicated in other brain regions to verify its universality.
Second, the HFS-induced LRTCs may be artificial, not identical to biological ones.Further studies are needed to investigate the differences between the artificial LRTCs and biological LRTCs and to confirm the benefits of HFS-induced LRTCs to the neural information processing in brain, as well as to the health of human beings.
Third, the efficiency and efficacy of brain stimulation therapy are greatly dependent on the pulse frequency of HFS [22,28].Although the 100 Hz frequency used in the present study is within the range of commonly used frequencies of DBS pulses, further studies are required to explore the fractal neuronal firing patterns induced by HFS with other pulse frequencies.In addition, further investigations using animal models of neurological disorders are needed to verify the relationship between the increase of LRTCs in the neuronal firing and the therapeutic efficacy of HFS.

Conclusions
The present study firstly showed that the HFS of afferent axonal fiber with a fixed IPI can enhance the fractal properties of the firing of individual neurons in the rat hippocampal CA1 region.The finding reveals a new potential mechanism of DBS to improve instead of to destroy the information processing of brain.It also suggests that a fractal analysis of neuronal firing may be useful to evaluate the efficacy of neuronal modulation therapies.

Figure 1 .
Figure 1.Electrode placements and recording segments of spontaneous neuronal potentials in the rat hippocampal CA1 region.(a) Schematic diagram of the implant positions of recording electrode (RE) in the pyramidal layer and the stimulation electrode (SE) in the Schaffer collaterals of CA1 region.(b) An example of raw recordings collected by four neighboring contacts (denoted by the red dots) on the RE array locating in the pyramidal layer.The signals of unit spikes are denoted by red arrow heads in the expanded insets.Unit spikes from the four channels were used for spike sorting.

Figure 2 .
Figure 2. Illustration of the generation of simulation sequences.(a) An example segment of an original spike sequence.The black vertical bars denote the 6 spikes, and the colored horizontal bars denote the 5 ISIs in the segment.Each ISI is assigned a random number, written in red.(b) The corresponding randomly shuffled sequence.(c) The corresponding surrogate sequence.

Figure 3 .
Figure 3. Increase of neuronal firing during 100 Hz HFS of Schaffer collaterals in the hippocampal CA1 region.(a) A typical example of neuronal response to 2 min 100 Hz HFS.The red bar denotes the stimulation period.The high-pass filtered signal (>500 Hz) shows the increase of MUA during the late 90 s HFS period without obvious PS events.(b)-(c) Left: the dynamic curves of firing rate of a pyramidal cell (Pyr) and an interneurons (Int), right: comparison of the firing rates of pyramidal cells and interneurons between the pre-stimulation baseline and during HFS.* * P < 0.01, paired t-test.

Figure 4 .
Figure 4. HFS-induced changes of firing patterns of pyramidal cells and interneurons.(a) Typical examples of raster plots of ISIs of a pyramidal cell at baseline and during HFS.The resolutions of ISI axis and time axis are 1 ms and 1 s, respectively.(b) Comparison of the average ISI curves of the total 31 pyramidal cells between baseline and during HFS.The blue and red curves represent the data from baseline recording and during HFS, respectively.(c)-(d) Corresponding plots as (a) and (b) for interneurons.(e)-(f) The auto-correlogram and unit-PSD of the firing of pyramidal cell and interneuron shown in (a) and (c).

Figure 5 .
Figure 5. Fractal structures in the firing rates of pyramidal cells and interneurons.(a)-(b) Examples of the curves of firing rate per second for a pyramidal cell (a) and an interneuron (b) during the late 90s HFS (top), along with graded enlargements of the curves at temporal scales of 30 s (middle) and 10 s (bottom), to illustrate the nested patterns across the three temporal scales.

Figure 6 .
Figure 6.Analysis of Fano factor in the firing of pyramidal cells and interneurons.(a) Examples of Fano factor F(T) curves of the firing of a pyramidal cell during baseline and HFS.(b) Comparisons of F(T) slope α between the baseline and during HFS.* * P < 0.01, paired t-test.(c)-(d) F(T) curves and comparisons of α for interneurons.* * P < 0.01, paired t-test.The red and yellow dashed lines on the curves in (a) and (c) represent the fitting lines.

Figure 7 .
Figure 7. Analysis of Hurst exponent in the firing of pyramidal cells and interneurons.(a) For the firing of pyramidal cells, example of the double-logarithmic plots of R d against the ISI number d (left) and the comparisons of Hurst exponents (H, right) between the baseline and during HFS.(b) Corresponding plots as (a) for interneurons.** P < 0.01, paired t-test.

Figure 8 .
Figure 8. Comparisons of the LRTCs among the experimental spike sequences and the simulation sequences.(a) For pyramidal cells, example of ISI curves of experimental spike sequences (red curve), randomly shuffled sequences (overlapped red curve), and surrogate sequences (green curve).(b) Corresponding plots as (a) for interneurons.(c)-(d) The auto-correlogram and unit-PSD for experimental spike sequences (top), randomly shuffled sequences (middle), and surrogate sequences (bottom) of pyramidal cells and interneurons showed in (a) and (b).(e) For pyramidal cells, left: example of F(T) curves of experiment spike sequence (red), randomly shuffled sequence (gray) and surrogate sequence (green).Right: comparison of the F(T) slope α between the experiment spike sequences (red) and the surrogate sequences (green).(f) Corresponding plots as (e) for interneurons.