Long-range temporal correlation in Auditory Brainstem Responses to Spoken Syllable/da/

The speech auditory brainstem response (sABR) is an objective clinical tool to diagnose particular impairments along the auditory brainstem pathways. We explore the scaling behavior of the brainstem in response to synthetic /da/ stimuli using a proposed pipeline including Multifractal Detrended Moving Average Analysis (MFDMA) modified by Singular Value Decomposition. The scaling exponent confirms that all normal sABR are classified into the non-stationary process. The average Hurst exponent is H = 0:77 ± 0:12 at 68% confidence interval indicating long-range correlation which shows the first universality behavior of sABR. Our findings exhibit that fluctuations in the sABR series are dictated by a mechanism associated with long-term memory of the dynamic of the auditory system in the brainstem level. The q-dependency of h(q) demonstrates that underlying data sets have multifractal nature revealing the second universality behavior of the normal sABR samples. Comparing Hurst exponent of original sABR with the results of the corresponding shuffled and surrogate series, we conclude that its multifractality is almost due to the long-range temporal correlations which are devoted to the third universality. Finally, the presence of long-range correlation which is related to the slow timescales in the subcortical level and integration of information in the brainstem network is confirmed.

The scalp-recorded Auditory Brainstem Response (ABR) is the most common auditory evoked potential that reflects the dynamics of the large populations of neurons along the auditory brainstem to simple acoustic sounds (e.g., tones, click) [1][2][3] . This evoked potential can be used for oto-neurological diagnosis, particularly possible lesions along the auditory brainstem pathways 4-6 . However, the ABR to click sounds cannot anticipate encoding of complex sounds because of the nonlinear dynamic behavior of the auditory system. Therefore, in the recent studies, more complex stimuli such as speech or music have been used to evaluate the behavior of the brainstem to the more complex stimuli 7,8 . The speech ABR (sABR) comprises transient and sustained responses. Transient and non-periodic characteristics of the stimulus are seen in the transient responses, while sustained time-locked responses including periodic characteristics are shown in the sustained responses [8][9][10] . The sABR signal is a non-invasive and objective tool that is suited for evaluating individuals with developmental and learning problems and also for studying the role of the brainstem in encoding complex sounds [9][10][11][12][13] .
The sABR signal is mostly contaminated by several types of external artifacts such as muscular movement, electroencephalogram, non-biological noises, and trends. A more general definition of trend is a part of a series representing a pattern or dominant behavior. As an illustration, monotonous and periodic features can be considered as well-known trends 14 . Over the past two decades, various methods in time and frequency domains have been used to analyze sABR using peak latency and amplitude, cross-correlation, Fourier analysis, and cross-phasogram 15,16 . However, these methods conceal any grasp of the exact dynamical features of the sABR signals. Furthermore, most quantitative methods to evaluate the temporal dynamics of sABR need the sABR series to be stationary in which the mean and variance of the signal do not change with time 17 . While, the recent studies have described non-linear and non-stationary dynamics of sABR signals, limiting the utility of these methods for examining generic features of sABR signals 18,19 . In addition, due to the complex nature of sABR signals and the presence of several types of noises, underlying signals are mimicked by noises and trends. Examining such data based on linear analysis is not reliable encouraging taking into account non-linear methods which are effective ways of explaining these complex relationships. Therefore, it is crucial to implement robust methods for analyzing sABR signals in order to remove the destructive effects of various trends and noises. The integration of information over large time scales is an essential ingredient of neural networks. The dynamics of this network are usually characterized by slow power-law decay or long-range temporal correlations. The long-range temporal correlations intuitively guarantee some memory about the past neural activities across different cortical and sub-cortical areas. Mathematical description of long-range correlation corresponds to the divergence of correlation function integration when the size of data becomes infinitely long 20 . It has been exhibited that various time series extracted from biological systems, including electrocardiography (ECG) 21,22 , electroencephalography (EEG) [23][24][25][26][27][28][29][30][31][32][33][34][35] , signal neuron discharge 36 and human gait 37 behave as scale-invariant processes. The power-law scaling studies have shown the universality of the behaviors in the complex biological systems [38][39][40][41][42] . Theuniversality indicates that there are properties for a large portion of systems which are independent of the dynamics of systems. In addition, it represents the fact that, a few essential factors are necessary to determine the scaling exponents of a complex system. Such scaling exponents characterize the behavior of a typical system. Consequently, various systems which may seem to be independent of each other can be classified in a category and they behave in a considerably similar manner 43 . For a typical system experiencing a phase transition, the special value of the parameter at which the system changes its phase is the system's critical point. For systems that exhibit universality, the closer the parameter is to its critical value, the less sensitively the order parameter depends on the details of the system. However, several other studies have emphasized extensive fluctuations in scaling exponents across subjects 44 . Similarly, it is revealed that long-range temporal correlations can change dramatically with small changes in the connectivity of the original networks 45 . So, it motivates us to evaluate whether there is a universality of the behaviors in sABR signals.
The sABR series is contaminated by non-stationary sources including trends and non-biological noises. To achieve reliable results, these spurious effects should be removed from the intrinsic fluctuations. Since there is no universal definition for the trends, various methods have been proposed to eliminate the trends from the underlying time series 14 . In order to examine long-range temporal correlation in non-stationary time series, some mathematical and computational methods have been proposed recently. Detrended fluctuation analysis (DFA) is the most well-known nonlinear methods for studying non-linear, non-stationary time series [46][47][48] . DFA is a scaling analysis that offers a quantitative parameter to characterize the correlation properties in non-stationary data by detrending the data on various time scales to remove spurious discovery of temporal correlations arising as an artifact of non-stationarity. Therefore, this approach allows us to determine an exact scaling exponent of the time series.
The generalized form of DFA which is known as Multifractal Detrended Fluctuations Analysis (MFDFA) is one of the best-known methods to capture multifractality in series 24,49 , and used in various fields, ranging from cosmic microwave background radiations 50 , sunspot fluctuations 51,52 , plasma fluctuations 53 , astronomy 54 , economic time series [55][56][57] , music 58,59 , traffic jamming 60 to image processing 61,62 and biological time series 21,[29][30][31]63 . However this approach is not suitable to completely remove sinusoidal and power-law trends [64][65][66] . In the presence of a sinusoidal trend superimposed on the data, the fluctuation functions derived by MFDFA or MFDMA contain at least one cross-over. Therefore, one can not assign a unique scaling exponent to clarify the fractality nature of underlying process 64,[67][68][69] . On the other hand, MFDFA contains discontinuity in its internal algorithm for capturing local trends leading to discrepancy in computed scaling exponents. Therefore, MFDMA has been proposed to quantify the statistical properties of mono (multi) fractal time series [70][71][72] . In addition, several robust methods have been also proposed to remove trends such as Fourier Detrended Fluctuations Analysis (FDFA) 73 , Adaptive Detrending method (AD) 74 , Singular Value Decomposition (SVD) 67,75 , and Empirical Mode Decomposition (EMD) 76 . In this study, we have used MFDMA and SVD-MFDMA approaches to remove crossover in our results. The crossover is a changing point in a typical scaling function. More precisely, in the presence of a crossover, one can assign different scaling exponents for two regimes. In our paper, we denote the point where a changing is recognized in the fluctuation function behavior by s × 77 . The SVD method can remove trends corresponding to the sinusoidal trends in the results 75,78 .
To the best of our knowledge, the multifractal analysis has not been used on the sABR series in the literature. In this paper, for the first time, MFDMA 72 method is used to capture the intrinsic multiscaling dynamics and assessment of universality behavior of sABR signals. We will also utilize SVD detrending algorithm to remove or at least decrease the influence of trends and noises as much as possible 75,79 . Then, the multifractal behavior of sABR signals is checked and in the case of multifractality, we determine the sources of multifractality based on the comparison the MFDMA results to those obtained via the MFDMA for shuffled and surrogate series. Furthermore, we evaluate the pattern of sABR signals and measure long-range temporal correlations in sABR series which has not been obtained as yet using other methods.
The rest of this paper is organized as follows. In section 2, the theoretical overview of MFDMA and SVD will be explained in details. In section 3, experimental results of the multifractal methods on sABR series are discussed. The multifractality of sABR series is assessed in this section. Finally, section 4 is devoted to the conclusion.

Results
In this section, we utilize MFDMA accompanying SVD method to evaluate the multifractal nature and the complexity of sABR data leading to clarify the statistical behavior of observed sABR series.
The upper panel of Fig. 1 depicts a typical normal observed sABR signal and associated trend determined by the SVD method. The corresponding residue between the observed time series and trend is plotted in the lower panel of Fig. 1.
We use DMA on the observed sABR data sets. Our results demonstrate that there is a crossover time scale approximately equal to the 128 Hz in the fluctuation function versus scale. This crossover corresponds to . ×s 7 8 msec and it is associated with the fundamental frequency of spoken syllable /da/. The scaling exponent for s < s × is h(q = 2) = 1.75 ± 0.15 demonstrating non-stationary nature of time series in this regime (Fig. 2). The slope of fluctuation function versus scale for s > s × is h(q = 2) = 1.31 ± 0.19. As illustrated in Fig. 2, SVD as the pre-processing algorithm can almost remove crossover which is almost devoted to the fundamental frequency. It is worth noting that mentioned behavior has been obtained for all q values as well as for other sABR samples used in this paper. In Fig. 3, we illustrate values of Hurst exponent determined by SVD-MFDMA for q = 2 for all normal sABR series including associated error-bar at 68% level of confidence. Our results show that Hurst exponent determined by SVD-MFDMA for q = 2 for all normal sABR series is H > 0.5 leading to long-range correlation behavior (equation (7)). The slope of 2  for small s is dominated by intrinsic fluctuations and consequently, it can be considered as a robust value for Hurst exponent [64][65][66] .
The generalized Hurst exponents averaged on 40 normal sABR signals applied by the MFDMA method for s < s × is shown in Fig. 4  obtain the nature of intrinsic fluctuations. Subsequently, to find reliable scaling exponents, we should implement a pre-processing method to reduce the effect of trends in sABR series.
The singularity spectrum f(α) of observed sABR series are shown in Fig. 5 (equation (9)). The strength of multifractality nature of sABR is examined by width of singularity spectrum, Δα = α max − α min . The values of Δα for all normal sABR series including associated error-bar at 68% level of confidence value are shown in Fig. 6. Now, we compare the fluctuation function of original sABR with the results of the corresponding shuffled and surrogate series to determine the nature of multifractality. In Fig. 7, the q dependence of the h(q) are shown for original, shuffled and surrogate series averaged on all normal sABR data sets. The sABR series values have been randomly shuffled to destroy the long-range correlations in the data. The shuffled data behaves as a monofractal signal and the h(q) does not change in general with q (Fig. 7). The sABR series have long-range correlations in different scales as it is obvious from the variance of h(q) values corresponding to different qs. The presence of the long-range power-law correlations exhibit the fractal dynamics of the under investigation system 80 .
The values of the Hurst exponent h(q = 2) for original, surrogate and shuffled sABR signal with MFDMA method averaging on 40 normal individuals are indicated in Table 1.

Discussion
In this study, we have examined the large-scale dynamics of neural oscillations in the normal auditory system in the brainstem level. The presence of the long-range power-law correlations exhibit the fractal dynamics of the under-investigation system 80 . Long-range temporal correlation has been obtained in a wide range of complex biological systems, including DNA sequences 81 , heart rate 82 , medullary sympathetic neurons in neurophysiology 83,84 , human brain oscillations 85 and long memory in human coordination 86 . All the studies lend considerable credence to the argument that an intrinsic part of the mechanism of neural information processing is the scale-free temporal correlation 44 . Various studies have revealed long-range temporal correlations in the neural oscillations of the normal human brain reflecting a memory of the underlying dynamics 42,80,87 . In the present study, we have used the non-invasively recorded sABR to examine whether such scaling behavior occurs in the subcortical auditory structure. We proposed a pipeline containing a combination of multifractal analysis (MFDMA) with singular value decomposition (SVD) to determine the multifractal characterization of the sABR series in the presence of trends. Our results based on the MFDMA, which does not have the discontinuity like the MFDFA method, reveal a  as a function of s averaged on all sABR series for s < s × is 1.75 ± 0.15, while this slope for s > s × is 1.31 ± 0.19 demonstrating that MFDMA cannot remove sinusoidal trends which are almost devoted to the fundamental frequency (Fig. 2). To determine the reliable generalized Hurst exponents of normal sABR signals, we applied SVD as a pre-processing algorithm to remove sinusoidal trends (Fig. 2). After applying SVD on our data sets, the MFDMA showed self-similarity features in normal sABR series. Also, the self-similarity parameters attained with the MFDMA method for s < s × and by SVD-MFDMA for all available scales were invariant across subjects revealing a universality class for normal sABR samples subjected to the long-range temporal correlated signal. The presence of the q−dependency of h(q) and Δα ≠ 0 indicate the presence of robust multifractality in the sABR series corresponding to the second universality property. The width of the singularity spectrum on average equates to Δα = 0.36 ± 0.22 at 1σ confidence interval quantifying the amount of multifractality of sABR series. The results demonstrate the complex structure of the sABR series revealing the presence of long-range correlation, which has been related to the slow timescales in subcortical level and integration of information in brainstem network 35 .
We have also found the source of multifractality in the sABR is almost related to the long-range temporal correlations confirmed by comparing generalized Hurst exponent of original sABR series with the results of the corresponding shuffled which is devoted to the third universality. This means that we have a generic property in all data sets from the source of multifractality point of view. By comparing the scaling exponents of original data with that of computing for surrogate series, we obtained that the contribution of probability distribution function broadness in complexity behavior is not significant. Our findings exhibited that fluctuations in the sABR series are dictated by a mechanism associated with long-term memory of the dynamic of the auditory system in the brainstem level revealing third universality class 80 . This evaluation provides deeper insight into the time series originating from the auditory system in the brainstem level. Further studies are required to examine how the long-range correlation is affected by different pathological conditions and evaluate the hypothesis that the MFDMA analysis of sABR series would reveal the difference in the long-range temporal correlation between pathological and normal conditions quantitatively. This pipeline can be used as a pre-processor to remove different types of noises and trends and also evaluation of other features of underlying sABR series. Evaluation of this approach as a diagnostic method is in progress.

participants and procedures. Forty volunteers from Iran University of Medical Sciences (18 women and
22 men), aged 20-28 years (mean ± SD = 22.77 ± 2.05) contributed to this experiment. Based on their self-report, they were right-handed and monolingual Persian speakers with no history of auditory, learning or neurologic problems. The study has been explained to subjects and then the informed consent has been collected from them. All procedures have been approved by the deputy of research review board, Iran University of medical sciences. It is worth mentioning that we followed the relevant approved guidelines in this research. The stimuli  consisted of a 40 ms speech syllable /da/ provided with Biologic Navigator instrument (Natus Medical Inc., San Carlos, CA, USA). The speech signals are presented to the subjects at a sampling frequency of 44.1 kHz through the computer's internal sound card. The fundamental frequency (F 0 ) of this syllable is 128 Hz. Complex nature combining both transient and sustained features, the universality of the syllable in most phonetic inventories and its research potential in hearing and learning disorders motivate researchers to utilize Consonant-vowel stimulus /da/ in auditory brainstem study 88 . The contact electrodes were positioned at the vertex (Cz) as noninverting, earlobes (inverting) and forehead (Fpz) as ground. The stimulus was presented monoaurally at 80 dB nHL via Biologic insert earphone (580-SINSER), with a repetition rate of 7.1/s. The sABR responses were recorded using 1024 digital sampling points over an 85.33 ms time window in the right ear. Finally, the average of total 6000 artifact free responses was collected from every volunteer. Figure 8 indicates a typical synthesis stimulus of /da/ and sABR signals used for further analysis. It is worth noting that, here we collect s-ABR series in the brainstem level (superior olivary complex and lateral lemniscus) to investigate the brainstem temporal encoding of speech.

MFDMA.
The multifractal detrended moving average analysis is a modified version of multifractal detrended fluctuation analysis (MFDFA). MFDMA method has been introduced to solve the presence of a discontinuity for fitting a polynomial at the boundary of each partition in MFDFA 72,89 . This algorithm is described as follows: • Consider sABR(i) as a time series where i = 1, …, N, then, construct the sequence of cumulative sums as follow: ( 1) (1 ) where s is the window size, and θ is the position parameter varying in the range [0, 1]. Hence, the moving average θ = 0 is called the backward moving average, θ = 0.5 corresponds to the centered moving average, and θ = 1 refers to the forward moving average 71,72 . Throughout this paper we set θ = 0.0 due to its robustness reported in various papers 72,90 • Construct the detrended data by subtracting computed moving average function from the cumulative series X as: where s − s 1 ≤ i ≤ N − s 1 • Now, divide ε(i) into N s = int[N/s] non-overlapping windows with the size of s and then we have fluctuation function as follows:  for q = 0 according to L'Höspital's rule, we have: To compute the reliable generalized Hurst exponent (equation 7), generally the likelihood statistics has been implemented 95,96 . For a Gaussian distribution, 68% confidence interval corresponds to the integration over variable on the interval represented by [−σ, +σ]. In this paper we computed the standard deviation of scaling exponents and reported by ±σ in order to clarify the statistical uncertainty in derived exponents (to make more sense see also 68 ). There is a non-linear dependency between τ(q) and q in multifractal signals. The Hurst exponent h(q) is related to the scaling exponent τ(q) by the following formula: To characterize multifractality more quantitatively, the so-called singularity spectrum is defined by Legendre transformation as 97,98 : and where α is the hölder exponent and f(α) denotes the dimension of the subset series that is characterized by α. In other word, the singularity spectrum reveals a value indicating the scaling behavior of the signal. The domain of Hölder spectrum, α ∈ [α min , α max ], becomes 68,90,99 : q max The width of the spectrum, Δα, gives a measure of the multifractality (complexity) of the spectrum. It can be explained as follows: max m in The higher values of the multifractality nature represent the complexity of underlying time series. The origin of multifractality in a sABR series can be confirmed by evaluating the corresponding shuffled and the surrogate time series. As explained in more details by J.W. Kantelhardt et al. 77 , generally, there is two different types of multifractality in a time series: (i) Multifractality due to a broadness probability density function of the time series. In Scientific RepoRts | (2019) 9:1751 | https://doi.org/10.1038/s41598-018-38215-w this case, the multifractality of the time series cannot be removed by random shuffling. (ii) Multifractality due to variations of long-range correlations in small and large fluctuations. Here, the probability density function of the time series can have a distribution with finite moments. The shuffling procedure destroys all long-range correlations. In another word, the auto-correlation of a shuffled series behaves as δ ′ ∼ − ′ ⟨ ⟩ x t x t t t ( ) ( ) ( ) Dirac . Hence, if the multifractality of the original time series belongs to the long-range correlation, the shuffled data will show the non-fractal scaling. However, the multifractality due to the broadness of the probability density function is not affected by the shuffling procedure. To determine the multifractality due to the broadness of the probability density function, the surrogate method is used. If the multifractality of the original time series belongs to a broad probability density function, h(q) obtained by the surrogate method will be independent of q. If sABR series has both kinds of multifractality, then the shuffled and surrogate series will have weaker multifractality than the original time series. Singular Value Decomposition (SVD). SVD method can be determined in the following steps 67,69,78 : (1): Construct a matrix which its elements are sABR series with the following order: where Σ d×(N−(d−1)τ) is a diagonal matrix and its elements are the singular values. We set 2p + 1 the dominant eigenvalues in the matrix Σ to zero to eliminate long periods. Then, set the dominant eigenvalues in the matrix Σ to zero. The resulting matrix will be Σ * . where 1 ≤ j ≤ N − (d − 1)τ and 1 ≤ i ≤ d. Finally, the cleaned sABR series will be used as input for MFDMA.

Data Availability
The datasets analyzed during the current study are not publicly available due to legal reasons but are available from the corresponding author on a reasonable request.