# On the Motion of Spikes: Turbulent-Like Neuronal Activity in the Human Basal Ganglia

- Science and Technology School, National University of San Martin, Buenos Aires, Argentina

Neuronal signals are usually characterized in terms of their discharge rate, a description inadequate to account for the complex temporal organization of spike trains. Complex temporal properties, which are characteristic of neuronal systems, can only be described with the appropriate, complex mathematical tools. Here, I apply high order structure functions to the analysis of neuronal signals recorded from parkinsonian patients during functional neurosurgery, recovering multifractal properties. To achieve an accurate model of such multifractality is critical for understanding the basal ganglia, since other non-linear properties, such as entropy, depend on the fractal properties of complex systems. I propose a new approach to the study of neuronal signals: to study spiking activity in terms of the velocity of spikes, defining it as the inverse function of the instantaneous frequency. I introduce a neural field model that includes a non-linear gradient field, representing neuronal excitability, and a diffusive term to consider the physical properties of the electric field. Multifractality is present in the model for a range of diffusion coefficients, and multifractal temporal properties are mirrored into space. The model reproduces the behavior of human basal ganglia neurons and shows that it is like that of turbulent fluids. The results obtained from the model predict that passive electric properties of neuronal activity, including ephaptic coupling, are far more relevant to the human brain than what is usually considered: passive electric properties determine the temporal and spatial organization of neuronal activity in the neural tissue.

## Introduction

The basal ganglia are a circuit of densely interconnected subcortical nuclei, whose disease is related to human movement disorders (Obeso et al., 2008). Current models of the basal ganglia are partly successful in the prediction of neurophysiologic alterations occurring in movement disorders, including Parkinson’s disease. However, no current model allows to predict and control deep brain stimulation (DBS), one of the major therapeutic approaches to Parkinson’s disease (Montgomery, 2016). The lack of a description of the complex properties of basal ganglia neuronal activity might be one cause of this failure. The fundament of the classic model of pathophysiology of Parkinson’s disease lies on an excessive discharge rate of the output nuclei of the basal ganglia: the globus pallidus interna (GPi), and substantia nigra reticulata (SNr) (Albin et al., 1989; DeLong, 1990). In primates, GPi and SNr neurons fire in a tonic manner, keeping the motor thalamus inhibited, and momentary reductions of their discharge rate facilitate movement. In humans, movement onset and velocity are positively correlated to the rate of discharge of the motor thalamus, and negatively to the rate of discharge of the GPi and SNr (Utter and Basso, 2008). Although these observations seem to confirm the classic model, evidence obtained during functional neurosurgery in patients with Parkinson’s disease is contradictory, since high frequency DBS increases the rate of discharge of the GPi instead of lowering it (Anderson et al., 2003; Montgomery, 2006; Andres and Darbin, 2017). Applying microdialysis techniques it has been shown that clinically effective DBS increases cyclic guanosine monophosphate (cGMP) in the putamen and GPi, supporting an increase of activity in these nuclei (Galati et al., 2006; Stefani et al., 2011a,b). This suggests that other features of neuronal activity than the rate of discharge need to be taken into account to understand the origin of parkinsonian symptoms and their amelioration by DBS. In summary, the predictions of the classic model of the basal ganglia are inconsistent with the effects of DBS and ablative surgery, raising concerns on the validity of the rate of discharge as a measure of the state of the basal ganglia (Tang et al., 2005; Montgomery and Gale, 2008).

Non-linear time series analysis has inspired an alternative view. A growing body of evidence indicates that complex properties, in particular fractality, are crucial for the understanding of basal ganglia pathophysiology (Obeso et al., 2000; Zheng et al., 2005; Darbin et al., 2006; Rasouli et al., 2006; Lim et al., 2010; Andres et al., 2011a; Alam et al., 2015). Complex and non-linear temporal properties are present in basal ganglia spike trains from rodents, primates and patients with movement disorders (Li et al., 2008; Lim et al., 2010; Hohlefeld et al., 2012, 2015; Darbin et al., 2013). The most popular non-linear measure for the characterization of parkinsonian spike trains is entropy, which diminishes due to pharmacologic or stimulation treatment, and increases due to alertness (Lafreniere-Roula et al., 2010; Lim et al., 2010; Andres et al., 2014a; Alam et al., 2015; Darbin et al., 2015). Importantly, complex properties cannot be understood isolated, but they relate to each other. For instance, entropy can be estimated wrongly if the multifractality of a system is not considered (Costa et al., 1997; Lyra and Tsallis, 1998; Latora et al., 2000). Hence the need to find methods that can be used to study complexity at this level, and equations to reproduce it.

This paper investigates some similarities between the fractal organization of spike trains from parkinsonian neurons and turbulent fluids. It presents a descriptive study of human pathology and a simulation study. In previous work, I hypothesized that neurons produce signals with correlations on different scales (Andres et al., 2014b). This hypothesis gives a reason to analyze high order structure functions. Since structure functions of order M remove polynomials of order M-1, low order linear processes are progressively eliminated by the analysis, highlighting the effect of stationary, high order, complex dynamics (Simonetti et al., 1985). This method was developed initially to describe the behavior of turbulent fields in the inertial range, and was used successfully both in experimental and simulation studies (Kolmogorov, 1941; Benzi et al., 1993; Briscolini et al., 1994). Later structure functions were applied to the study of physiologic signals (Lin and Hughson, 2001; Nanni and Andres, 2017). Although the length of the time series analyzed must increase depending on the time lag studied, the study of structure functions is robust to short time series in comparison to other tools (Schulz-DuBois and Rehberg, 1981; Van de Water and Herweijer, 1999; Monin and Yaglom, 2013; Nanni and Andres, 2017). Time series of a length of thousands of data points have been analyzed with high order structure functions in different fields, and it was shown that this is a suitable method to look for non-linearity in the exponent function, i.e., to distinguish monofractality from multifractality (Yu et al., 2003; Huang et al., 2011). The paper is organized as follows. First it is shown that multifractality can be measured in neuronal signals from the GPi analyzing structure functions of increasing order. Then a neural field model is introduced inspired by Burgers’ equation, a well-known equation in the field of fluid dynamics (Burns et al., 1998). Previous evidence indicated that this kind of equations could be useful to study mathematical properties of neuronal activity, which is analyzed here (Andres et al., 2011b). The model is based on physical properties of neurons, and it captures essential features of neuronal activity of the human brain, reproducing multifractality as observed in the neuronal activity of the human, parkinsonian basal ganglia.

## Methods and Results

### Patients and Clinical Recordings

Six patients fulfilling the clinical criteria of the United Kingdom’s Parkinson’s Disease Society Brain Bank for idiopathic Parkinson’s (UK-PDS-BB) disease, Hoehn and Yahr (1967) IV, underwent stereotactic neurosurgery and were included in this study. All patients were on chronic treatment with L-DOPA, presented similar motor affectation, with severe dyskinesias and motor fluctuations, and fulfilled the Core Assessment Program for Surgical Interventional Therapies in Parkinson’s Disease (CAPSIT) criteria for inclusion in the surgery program (Defer et al., 1999). All patients signed informed consent prior to the procedure, previously reviewed and approved by the institutional ethics committee, and were without medication at the time of the surgery and during data acquisition. Age, gender, and other details are not shown to protect patients’ privacy.

In total, 22 GPi recordings obtained from 6 patients were analyzed. Details about the recording procedure have been already published elsewhere (Andres et al., 2011a). Briefly, microrecording, stimulation, and neurosurgical procedures were performed in patients awake, under local anesthesia. Surgical targets were planned employing magnetic resonance imaging (MRI) and using a Leksell stereotactic system (Series G, Elekta, Sweden). Microrecordings of neuronal activity were obtained only during the surgery, after the GPi was identified by an expert. Platinum/iridium (Pt/Ir 80/20%) microelectrodes with nominal impedance of 0.8–1.2 megohms (mTSPBN-LX1, FHC Inc.) were used. A differential amplifier with a built-in impedance meter (FHC IS-AM-00-01 Iso-Xcell 3 Amplifier) and an isolated stimulus generator (FHC IS-PL-06 Isolated Bipolar Pulsar Stimulator; FHC, Bowdoinham, ME, United States) were connected to a preamplifier remote probe mounted onto a motorized microdrive (FHC 65-00-1 Stepper Drive and ST-M0-00 TMS Controller), located near the electrode tip to minimize pickup of electrical noise. A dedicated acquisition system (1401plus, CED) was employed to amplify, and digitize the signal, filtering with second order 300–5,000 Hz bandpass and 50 Hz line notch analog filters. The sampling rate was 50 kHz. The total amplification including probe was set at gain ×10000 (checked with a built-in calibration signal of 1 mV at the beginning of each surgery).

### Signal Processing

Signals were processed off-line. To obtain single cell recordings, spikes were extracted from raw signals and separated into classes using wavelet analysis and clustering (Quiroga et al., 2004; Andres et al., 2014a). From these single cell data, time series of interspike intervals (ISI) were constructed. Time series had a mean length of 5668 ± 773 ISI (mean ± standard error of the mean, SEM). From these time series, the temporal structure function was computed as follows. An interval I(t) is defined as the *t*th ISI, which is used to calculate the increment ΔI(τ) = I(t + τ) - I(t), being τ the time lag, between 1 and 1000. The temporal structure function S_{q}(τ) is

where ⟨⋅⟩ is the statistical average, and q is the order of the structure function, an integer between 1 and 30. Next, S_{q} was smoothed by applying a running average over a 30 points window to obtain a smoothed structure function, ${\text{S}}_{\text{q}}^{*}$ (30 points shorter than S_{q}). The behavior of interest here is the scaling behavior of ${\text{S}}_{\text{q}}^{*}$(τ). Given the general relationship

the scaling exponent is a function, ζ_{τ}(q), which defines the fractal properties of the signal under study. For random processes, ζ_{τ}(q) = 0. In the case of monofractals the exponent function grows linearly, ζ_{τ}(q) = qζ(1). If the signal is multifractal, ζ_{τ}(q) is a non-linear function (Lin and Hughson, 2001). Note that different authors use different notations. The function that is called here ζ_{τ}(q) is called τ(q) in Touchette’s work, according to whom the spectrum of singularities f(α) does not need to be concave (Touchette and Beck, 2006). A concave f(α) can be calculated as the Legendre transform of ζ_{τ}(q), whereas a not concave f(α) cannot. In any case, a concave ζ_{τ}(q) implies multifractality.

To find ζ_{τ}(q), log(${\text{S}}_{\text{q}}^{*}$(τ)) vs. log(τ) was considered, and a linear (scaling) region was looked for. The scaling region was defined as the longest range of τ fitting a linear function with non-zero slope (ζ≠0), for at least 1 ≤ q ≤ 10. Linear regressions were considered acceptable if R^{2} ≥ 0.6, and regressions with regression coefficients smaller than that were discarded. A scaling region of the temporal structure function was found for a minimum range of 1 ≤ q ≤ 10 in every neuron. In 15 neurons, the fitting was sufficiently good for q up to 30 (q_{max} = 24 ± 2, mean ± SEM). The length of the scaling interval over τ was 79 ± 25; mean ± SEM.

Results for the analysis of neuronal recordings obtained from the patients with Parkinson’s disease are shown in Figure 1. A multifractal organization of spike trains was found in 82% of neurons, evidenced by non-linearity of the temporal exponent function (ζ_{τ}(q)), obtained from structure functions of increasing order.

**FIGURE 1.** Multifractal spectra of human neurons. Temporal multifractal spectra of sample neuronal recordings, obtained from the globus pallidus interna (GPi) of patients with Parkinson’s disease with the structure function method. The non-linearity of ζ_{τ}(q) indicates multifractality. The function ζ_{τ}(q) is calculated from temporal structure functions of increasing order q, S_{q}(τ). Since structure functions are built from time series of neuronal activity, the spectra indicate multifractal organization in the temporal domain.

### Model and Simulations

To reproduce the behavior observed in the clinical data, a model is introduced that takes the form of a neural field equation, as follows:

where u(x,t) represents the velocity of spikes (a function of time t and space x), δ is a diffusion coefficient and W is a stochastic drive with amplitude *a*.

The model describes neuronal activity in terms of the velocity of spikes. While spikes travel with constant velocity along axons, the same does not hold for neural tissue or neural networks at larger scales. Indeed, consider the following simplified situation. Take a single spike S that arrives to a neuronal soma whose membrane potential V is below the triggering threshold 𝜗. If this spike induces an excitation with an amplitude e in the neuron, such that the V + e > 𝜗, then the neuron will fire in response, which is equivalent to saying that the spike S has traveled through the neuron (a node of the network). Importantly, the condition V + e > 𝜗 can be satisfied for a range of potentials V and, since summation is a time-consuming mechanism, the spike triggered in response to the excitation e will be fired faster for membrane potentials that are closer to the threshold, giving place to a higher velocity of transmission u. Not only is velocity variable, but it is non-linearly so: the process of summation is typically a non-linear function of time, introducing non-linear variability in the profile of u. Even more, due to the close relation between the variables u (velocity of spikes) and f (neuronal frequency), it is expected that any non-linearity of u will be reflected on the time series of interspike intervals. The relation between both variables is simple: since u = Δx/Δt and f = 1/Δt, then u = f ⋅Δx. Consequently, for distances Δx = 1 the module of the velocity |u| has the same value than the frequency. Regarding the time intervals between spikes (interspike intervals I, ISI), they are equal to the inverse of the instantaneous frequency of discharge, and therefore |u| = 1/I.

Computer simulations were run with custom MATLAB and R code. The model was studied numerically in 1D. The discretization scheme selected was a Crank-Nicholson scheme, which is suitable for this kind of study (Kadalbajoo and Awasthi, 2006). The space steps and time steps were Δx = Δt = 0.001, and u(t) = 0 at the boundaries of the spatial domain. The initial condition considered was a random distribution of u(x) between -1 and 1, with a power spectrum obeying the scaling law P(ω(u)) = ω^{β}, where β < 0. Following Hayot and Jayaprakash (1996), the same kind of colored noise W was used as a drive to the equation, with an amplitude a = 10^{-6}. The temporal evolution of u(x) was studied for 1000 time and space steps and for a range of diffusion coefficients 10^{-3} < δ < 0.11, varying δ in increments of 0.02. To make the data obtained from the model more directly comparable to clinical data, the inverse of the absolute value of the velocity was considered, obtaining a variable equivalent to I(x,t). This variable $\text{I(x,t)}=\frac{1}{|\text{u|}}$ was analyzed applying the same algorithm that was described in point 2.2 for the analysis of clinical ISI time series to obtain a temporal structure function S_{q}(τ) at fixed points in space. Finally, a spatial structure function S_{q}(X) was computed in a similar manner, but considering adjacent points of u for fixed times. Scales between 1 and 50 were used to compute the structure functions both in space and time. A multifractal spectrum was looked for both in the space (ζ_{X}(q)) and time domains (ζ_{τ}(q)), by analyzing the behavior of S_{q}(τ) ∼ τ^{(q)} and S_{q}(X) ∼ X^{(q)} for q varying between 0 and 2 in steps of 0.1.

The diffusion coefficient is the critical parameter in the model. Figure 2 (left column) shows the results of computer simulations for increasing diffusion coefficient δ. At the white areas the module of the velocity of spikes (equivalent to frequency) is lowest, $\frac{1}{|\text{u(x,t)|}}\ge {10}^{8}$, in opposition to the black areas, where it is highest. As the diffusion coefficient increases, white areas are enlarged, as the total area of high velocity diminishes. At the same time the whole integration domain becomes more homogeneous, with individual areas of high velocity becoming wider. This indicates that diffusion not only contributes to the dissipation of energy (lowering the global velocity of spikes transmission), but it also modifies the temporal and spatial organization of the neuronal activity across the neural field. Importantly, hallmarks of multifractality appear both in space and time for the whole range of δ investigated. The middle and right columns of Figure 2 show typical examples of the temporal and spatial multifractal spectra obtained from the model, ζ_{τ}(q) and ζ_{X}(q), respectively. Multifractality is characteristic of the model, with non-linear exponent functions in the spatial as well as temporal domains.

**FIGURE 2.** Computer simulations. Left column: Time evolution of the velocity of spikes, u(x,t), as the diffusion coefficient δ increases (from top to bottom), with time on the vertical axis and space on the horizontal axis. White areas represent the parts of the integration domain where the module of the velocity of spikes is below an arbitrary limit (10^{8}), in opposition to black areas, where it is higher than this limit. As the diffusion coefficient increases, white areas are enlarged, as the total velocity diminishes across the integration domain. Middle column: Sample temporal multifractal spectra ζ_{τ}(q) obtained from temporal structure functions of increasing order, at fixed spatial points. Non-linearity indicates temporal multifractality. Right column: Sample spatial multifractal spectra ζ_{x}(q) obtained from spatial structure functions of increasing order, at fixed times. Non-linearity indicates spatial multifractality.

## Conclusion

In the neural field studied, complex geometrical properties arise from a combination of diffusive properties, a gradient field and a stochastic drive, in a manner similar to that of fluids undergoing turbulence. While this is an analogy, and any physical implications should be considered with care, some inferences might be interesting. For instance, the model provides a mechanism for the generation of the complex patterns observed in parkinsonian neuronal activity, i.e., that high-pass filtering of turbulent-like velocity signals produces intermittent bursts in non-linear systems with stochastic drive, such as the present one (Frisch and Morf, 1981). Here the stochastic drive can be interpreted as a high dimensional environmental input to the neural network, but the gradient field represents an intrinsic property of the network, i.e., neuronal excitability at a specific point of space, determining the change of the spikes’ velocity. Regarding the diffusive term, diffusive properties are a physical consequence of neuronal processes, given the diffusive nature of electric fields. Previous modeling work on the basal ganglia showed that a diffusion coefficient is a critical parameter for the control of information transmission and information deterioration in the parkinsonian GPi (Andres et al., 2014b). However, the role of these so-called passive electric properties of neuronal activity is not usually considered in neurophysiology studies, which more commonly focus on spiking activity solely. In the neural field introduced here: (1) multifractality of neuronal spike trains is related to diffusive properties, and (2) multifractality of temporal activity is reflected on the spatial domain. The model predicts that passive (diffusive) properties of neuronal activity determine the structure of temporal and spatial neuronal activity in the basal ganglia, and must be considered in the study and treatment of movement disorders.

## Ethics Statement

The experimental protocol was revised and approved by FLENI Ethics Committee, Buenos Aires, Argentina.

## Author Contributions

DA was responsible for the full content of this article.

## Funding

The work of DA was supported by the Ministry of Science and Innovative Production, Argentina, and by National University of San Martin. The studies of this article were conducted with the support of Society in Science, The Branco-Weiss Fellowship, administered by ETH Zurich.

## Conflict of Interest Statement

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

The author acknowledges the work of the personnel at FLENI Institute, without whom the making of this paper would not have been possible.

## References

Alam, M., Sanghera, M. K., Schwabe, K., Lutjens, G., Jin, X., Song, J., et al. (2015). Globus pallidus internus neuronal activity: a comparative study of linear and non-linear features in patients with dystonia or Parkinson’s disease. *J. Neural Transm.* 123, 231–240. doi: 10.1007/s00702-015-1484-3

Albin, R. L., Young, A. B., and Penney, J. B. (1989). The functional anatomy of basal ganglia disorders. *Trends Neurosci.* 12, 366–375. doi: 10.1016/0166-2236(89)90074-X

Anderson, M. E., Postupna, N., and Ruffo, M. (2003). Effects of high-frequency stimulation in the internal globus pallidus on the activity of thalamic neurons in the awake monkey. *J. Neurophysiol.* 89, 1150–1160. doi: 10.1152/jn.00475.2002

Andres, D. S., Cerquetti, D., and Merello, M. (2011a). Finite dimensional structure of the GPI discharge in patients with Parkinson’s disease. *Int. J. Neural. Syst.* 21, 175–186. doi: 10.1142/S0129065711002778

Andres, D. S., Cerquetti, D. F., and Merello, M. (2011b). Turbulence in globus pallidum neurons in patients with Parkinson’s disease: exponential decay of the power spectrum. *J. Neurosci. Methods* 197, 14–20. doi: 10.1016/j.jneumeth.2011.01.022

Andres, D. S., Cerquetti, D., Merello, M., and Stoop, R. (2014a). Neuronal entropy depends on the level of alertness in the parkinsonian globus pallidus in vivo. *Front. Neurol.* 5:96. doi: 10.3389/fneur.2014.00096

Andres, D. S., Gomez, F., Ferrari, F. A., Cerquetti, D., Merello, M., Viana, R., et al. (2014b). Multiple-time-scale framework for understanding the progression of Parkinson’s disease. *Phys. Rev. E* 90:062709.

Andres, D. S., and Darbin, O. (2017). Complex dynamics in the basal ganglia: health and disease beyond the motor system. *J. Neuropsychiatry Clin. Neurosci.* 30, 101–114. doi: 10.1176/appi.neuropsych.17020039

Benzi, R., Ciliberto, S., Tripiccione, R., Baudet, C., Massaioli, F., and Succi, S. (1993). Extended self-similarity in turbulent flows. *Phys. Rev. E* 48, R29–R32. doi: 10.1103/PhysRevE.48.R29

Briscolini, M., Santangelo, P., Succi, S., and Benzi, R. (1994). Extended self-similarity in the numerical simulation of three-dimensional homogeneous flows. *Phys. Rev. E* 50, R1745–R1747. doi: 10.1103/PhysRevE.50.R1745

Burns, J., Balogh, A., Gilliam, D., and Shubov, V. (1998). Numerical stationary solutions for a viscous Burgers’ equation. *J. Math. Syst. Estim. Control* 8, 1–16.

Costa, U. M. S., Lyra, M. L., Plastino, A. R., and Tsallis, C. (1997). Power-law sensitivity to initial conditions within a logistic like family of maps: fractality and nonextensivity. *Phys. Rev. E* 56:245. doi: 10.1103/PhysRevE.56.245

Darbin, O., Adams, E., Martino, A., Naritoku, L., Dees, D., and Naritoku, D. (2013). Non-linear dynamics in parkinsonism. *Front. Neurol.* 4:211. doi: 10.3389/fneur.2013.00211

Darbin, O., Jin, X., Von Wrangel, C., Schwabe, K., Nambu, A., Naritoku, D. K., et al. (2015). Neuronal entropy-rate feature of entopeduncular nucleus in rat model of Parkinson’s disease. *Int. J. Neural Syst.* 26:1550038. doi: 10.1142/S0129065715500380

Darbin, O., Soares, J., and Wichmann, T. (2006). Nonlinear analysis of discharge patterns in monkey basal ganglia. *Brain Res.* 1118, 84–93. doi: 10.1016/j.brainres.2006.08.027

Defer, G. L., Widner, H., Marie, R. M., Remy, P., and Levivier, M. (1999). Core assessment program for surgical interventional therapies in Parkinson’s disease (CAPSIT-PD). *Mov. Disord.* 14, 572–584. doi: 10.1002/1531-8257(199907)14:4<572::AID-MDS1005>3.0.CO;2-C

DeLong, M. R. (1990). Primate models of movement disorders of basal ganglia origin. *Trends Neurosci.* 13, 281–285. doi: 10.1016/0166-2236(90)90110-V

Frisch, U., and Morf, R. (1981). Intermittency in nonlinear dynamics and singularities at complex times. *Phys. Rev. A* 23, 2673–2705. doi: 10.1103/PhysRevA.23.2673

Galati, S., Mazzone, P., Fedele, E., Pisani, A., Peppe, A., Pierantozzi, M., et al. (2006). Biochemical and electrophysiological changes of substantia nigra pars reticulata driven by subthalamic stimulation in patients with Parkinson’s disease. *Eur. J. Neurosci.* 23, 2923–2928. doi: 10.1111/j.1460-9568.2006.04816.x

Hayot, F., and Jayaprakash, C. (1996). Multifractality in the stochastic Burgers equation. *Phys. Rev. E* 54, 4681–4684. doi: 10.1103/PhysRevE.54.4681

Hoehn, M. M., and Yahr, M. D. (1967). Parkinsonism: onset, progression and mortality. *Neurology* 17, 427–442. doi: 10.1212/WNL.17.5.427

Hohlefeld, F., Ehlen, F., Tiedt, H., Krugel, L., Horn, A., Kühn, A., et al. (2015). Correlation between cortical and subcortical neural dynamics on multiple time scales in Parkinson’s disease. *Neuroscience* 298, 145–160. doi: 10.1016/j.neuroscience.2015.04.013

Hohlefeld, F., Huebl, J., Huchzermeyer, C., Schneider, G. H., Schönecker, T., Kühn, A., et al. (2012). Long-range temporal correlations in the subthalamic nucleus of patients with Parkinson’s disease. *Eur. J. Neurosci.* 36, 2812–2821. doi: 10.1111/j.1460-9568.2012.08198.x

Huang, Y., Schmitt, F. G., Hermand, J. P., Gagne, Y., Lu, Z., and Liu, Y. (2011). Arbitrary-order Hilbert spectral analysis for time series possessing scaling statistics: comparison study with detrended fluctuation analysis and wavelet leaders. *Phys. Rev. E* 84:016208.

Kadalbajoo, M. K., and Awasthi, A. (2006). A numerical method based on Crank-Nicolson scheme for Burgers’ equation. *Appl. Math. Comput.* 182, 1430–1442. doi: 10.1016/j.amc.2006.05.030

Kolmogorov, A. N. (1941). The local structure of turbulence in an incompressible fluid at very high Reynolds numbers. *Proc. R. Soc. Lond. A* 434, 9–13. doi: 10.1098/rspa.1991.0075

Lafreniere-Roula, M., Darbin, O., Hutchison, W. D., Wichmann, T., Lozano, A. M., and Dostrovsky, J. O. (2010). Apomorphine reduces subthalamic neuronal entropy in parkinsonian patients. *Exp. Neurol.* 225, 455–458. doi: 10.1016/j.expneurol.2010.07.016

Latora, V., Baranger, M., Rapisarda, A., and Tsallis, C. (2000). The rate of entropy increase at the edge of chaos. *Phys. Lett. A* 273, 97–103. doi: 10.1016/S0375-9601(00)00484-9

Li, W., Jia, D., Wang, J. L., Liang, Q., Jian, Z., Wang, X. L., et al. (2008). Deterministic dynamics in neuronal discharge from pallidotomy targets. *J. Int. Med. Res.* 36, 979–985. doi: 10.1177/147323000803600514

Lim, J., Sanghera, M. K., Darbin, O., Stewart, R. M., Jankovic, J., and Simpson, R. (2010). Nonlinear temporal organization of neuronal discharge in the basal ganglia of Parkinson’s disease patients. *Exp. Neurol.* 224, 542–544. doi: 10.1016/j.expneurol.2010.05.021

Lin, D. C., and Hughson, R. L. (2001). Modeling heart rate variability in healthy humans: a turbulence analogy. *Phys. Rev. Lett.* 86, 1650–1653. doi: 10.1103/PhysRevLett.86.1650

Lyra, M. L., and Tsallis, C. (1998). Nonextensivity and multifractality in low-dimensional dissipative systems. *Phys. Rev. Lett.* 80, 53–56. doi: 10.1103/PhysRevLett.80.53

Monin, A. S., and Yaglom, A. M. (2013). *Statistical Fluid Mechanics: Mechanics of Turbulence, Vol. 2*. North Chelmsford, MA: Courier Corporation.

Montgomery, E. B. Jr. (2016). Modeling and theories of pathophysiology and physiology of the basal ganglia–thalamic–cortical system: critical analysis. *Front. Hum. Neurosci.* 10:469. doi: 10.3389/fnhum.2016.00469

Montgomery, E. B. Jr., and Gale, J. T. (2008). Mechanisms of action of deep brain stimulation(DBS). *Neurosci. Biobehav. Rev.* 32, 388–407. doi: 10.1016/j.neubiorev.2007.06.003

Montgomery, E. B. (2006). Effects of GPi stimulation on human thalamic neuronal activity. *Clin. Neurophysiol.* 117, 2691–2702. doi: 10.1016/j.clinph.2006.08.011

Nanni, F., and Andres, D. S. (2017). Structure function revisited: a simple tool for complex analysis of neuronal activity. *Front. Hum. Neurosci.* 11:409. doi: 10.3389/fnhum.2017.00409

Obeso, J. A., Rodriguez-Oroz, M. C., Benitez-Temino, B., Blesa, F. J., Guridi, J., Marin, C., et al. (2008). Functional organization of the basal ganglia: therapeutic implications for Parkinson’s disease. *Mov. Disord.* 23(Suppl. 3), S548–S559. doi: 10.1002/mds.22062

Obeso, J. A., Rodriguez-Oroz, M. C., Rodriguez, M., Lanciego, J. L., Artieda, J., Gonzalo, N., et al. (2000). Pathophysiology of the basal ganglia in Parkinson’s disease. *Trends Neurosci.* 23(10 Suppl.), S8–S19.

Quiroga, R. Q., Nadasdy, Z., and Ben-Shaul, Y. (2004). Unsupervised spike detection and sorting with wavelets and super paramagnetic clustering. *Neural Comput.* 16, 1661–1687. doi: 10.1162/089976604774201631

Rasouli, G., Rasouli, M., Lenz, F. A., Verhagen, L., Borrett, D. S., and Kwan, H. C. (2006). Fractal characteristics of human Parkinsonian neuronal spike trains. *Neuroscience* 139, 1153–1158. doi: 10.1016/j.neuroscience.2006.01.012

Schulz-DuBois, E., and Rehberg, I. (1981). Structure function in lieu of correlation function. *Appl. Phys.* 24, 323–329. doi: 10.1007/BF00899730

Simonetti, J. H., Cordes, J. M., and Heeschen, D. S. (1985). Flicker of extragalactic radio sources at two frequencies. *Astrophys. J.* 296, 46–59. doi: 10.1086/163418

Stefani, A., Fedele, E., Pierantozzi, M., Galati, S., Marzetti, F., Peppe, A., et al. (2011a). Reduced GABA content in the motor thalamus during effective deep brain stimulation of the subthalamic nucleus. *Front. Syst. Neurosci.* 5:17. doi: 10.3389/fnsys.2011.00017

Stefani, A., Fedele, E., Vitek, J., Pierantozzi, M., Galati, S., Marzetti, F., et al. (2011b). The clinical efficacy of L-DOPA and STN-DBS share a common marker: reduced GABA content in the motor thalamus. *Cell Death Dis.* 2:e154. doi: 10.1038/cddis.2011.35

Tang, J. K., Moro, E., Lozano, A. M., Lang, A. E., Hutchison, W. D., Mahant, N., et al. (2005). Firing rates of pallidal neurons are similar in Huntington’s and Parkinson’s disease patients. *Exp. Brain Res.* 166, 230–236. doi: 10.1007/s00221-005-2359-x

Touchette, H., and Beck, C. (2006). Nonconcave entropies in multifractals and the thermodynamic formalism. *J. Stat. Phys.* 125, 455–471. doi: 10.1007/s10955-006-9174-z

Utter, A. A., and Basso, M. A. (2008). The basal ganglia: an overview of circuits and function. *Neurosci. Biobehav. Rev.* 32, 333–342. doi: 10.1016/j.neubiorev.2006.11.003

Van de Water, W., and Herweijer, J. (1999). High-order structure functions of turbulence. *J. Fluid Mech.* 387, 3–37. doi: 10.1017/S0022112099004814

Yu, C. X., Gilmore, M., Peebles, W. A., and Rhodes, T. L. (2003). Structure function analysis of long-range correlations in plasma turbulence. *Phys. Plasmas* 10, 2772–2779. doi: 10.1063/1.1583711

Keywords: structure function, Parkinson’s disease, neuronal activity, turbulence modeling, neuronal modeling, basal ganglia, complexity, non-linear dynamics

Citation: Andres D (2018) On the Motion of Spikes: Turbulent-Like Neuronal Activity in the Human Basal Ganglia. *Front. Hum. Neurosci.* 12:429. doi: 10.3389/fnhum.2018.00429

Received: 27 April 2018; Accepted: 02 October 2018;

Published: 24 October 2018.

Edited by:

Ioan Opris, University of Miami, United StatesReviewed by:

Sauro Succi, Consiglio Nazionale delle Ricerche (CNR), ItalyYa-tang Li, California Institute of Technology, United States

Copyright © 2018 Andres. 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: Daniela Andres, dandres@unsam.edu.ar; daniela.s.andres@gmail.com