**Figure 1** shows the averaged (across all sensors into a single value) developmental trajectory of the spectral power (**Figure 1**, upper panels) and MSE (**Figure 1**, lower panels) across five infants aged from 5 to 36 months old. All infants demonstrated an increase in the MSE value with age. The Jonckheere-Terpstra test was used to test for an age-related trend in MSE values, and statistically significant age-related trends were identified for coarse time scales (31–40 ms, scales: 16–20) (TJT = 459.0, standard error = 34.8, z = 5.6, p < 0.001). **Figure 2** shows the averaged developmental trajectory in each time scale bin (**Figure 2A**) and the topography of MSE values across different ages (**Figure 2B**). A more detailed examination of our results revealed that the remarkable increase in MSE identified for longer time scales (31–40 ms, scale: 16–20) was predominantly observed at ages up to 15 months and was found across brain regions (**Figure 2**). After 15 months of age, this increase tended to slow. However, the power spectral analysis also showed an increase in power in the theta and alpha bands. This increase was more prominent after 15 months of age, while the MSE change was more prominent during the earlier infancy periods. Regarding the shorter time scales (2–10 ms, scale: 1–5), the developmental trajectory of MSE varied across subjects. For instance, some infants showed constant MSE values across development, whereas one infant showed a gradual decrease (**Figure 1**, bottom panels).

In the surrogate analysis, we found region- and scale-specific entropy alterations in the surrogate data, which may suggest an inherent non-linearity in the MEG data (data not shown). Specifically, in the surrogate data, the SampEn increased near the frontal and temporo-occipital regions. Interestingly, this regionspecific SampEn alteration was more prominent for smaller temporal scales (less than 20 ms) and was frequently identified during early infancy (5–10 months of age).

### DISCUSSION

The neurodevelopmental trajectory of infancy has received much attention because infancy is a critical period of brain development in which cognitive and behavioral abilities are enhanced (Cao et al., 2017) and neurodevelopmental disorders, such as autism spectrum disorder (ASD), are predicted to develop. This is the first longitudinal investigation of how brain signal complexity, which represents neural system dynamics, changes during infancy. The analysis revealed scale-dependent developmental trajectories of MEG signal complexity. Specifically, we found an increase in signal complexity for longer time scales, whereas the changes in complexity varied across infants for shorter time scales.

Many studies have investigated age-related signal complexity changes from late childhood into adulthood. Polizzotto et al. (2015) examined the MSE of resting-state EEG results in healthy subjects aged 8–22 years old. They reported an age-related increase in entropy in lower scales and a decrease in entropy for higher scales. McIntosh et al. (2008) calculated MSE changes in EEG during a face recognition visual memory task in children (8–15 years old) and young adults (20–33 years old). They found an age-related increase in EEG complexity that was significantly correlated with the accuracy of task performance. This observation was replicated by the same group using MEG (Misic et al., 2010), confirming the characteristic shape of the MSE curve and its prominent task-dependent increase during development. We have also demonstrated an age-related increase in MEG signal complexity. However, enhanced complexity was identified in children with ASD, particularly in earlier childhood (Takahashi et al., 2016). Compared to the changes that occur during the period from childhood to adolescence, brain signal complexity during infancy has been addressed by few studies. Lippe et al. (2009) investigated EEG signal complexity in response to visual and auditory stimulation in children ranging from

FIGURE 1 | Array plot showing the developmental trajectory of spectral power (top panels) and MSE (bottom panels) for each infant and their average (X-axis, age, months; Y-axis, frequency and time scales).

months of age).

1 month to 5 years of age. They found a task-dependent increase in EEG complexity with aging. However, these studies were based on a cross-sectional study design. Despite the small number of subjects, a unique aspect of this study is that we longitudinally investigated the development of MEG complexity during infancy.

In these contexts, our study provides the longitudinal underpinnings for the concept of significant shifts in brain signal complexity with aging (Garrett et al., 2013). Notably, for higher scales, we captured a robust developmental MSE profile across infants and across MEG sessions despite conditional inconsistencies (i.e., selected videos, emotions, or physical conditions), which may indicate the potential usefulness of MSE as a reliable and clinically useful trait biomarker of the infant brain. For instance, we have demonstrated a linear age-related increase in complexity at higher scales across 40–110 monthold children (Takahashi et al., 2016). Additionally, enhanced brain signal variability was observed in children with ASD, which was conformed for younger children. On the other hand, Bosl et al. (2011) examined resting-state EEG complexity by MSE in typically developing infants and infants with a high risk of ASD

across the ages of 6–24 months, and they found consistently lower EEG complexity at higher scales in the high-risk group, particularly at 9–12 months of age. These inconsistent results may be attributed to the different age ranges of these two previous studies on children with ASD.

Considering biological background, the observed rapid increase in MEG complexity in the high time scale (i.e., lower frequency range) at approximately 5–15 months old might demonstrate the development of long-range network-related cognitive processing. Given that long-range communication between multiple brain areas is driven by slow waves (i.e., theta and beta waves) (Wang, 2010), MEG complexity in a high time scale (i.e., lower frequency range) may be useful and a non-invasive biomarker of brain maturation in infants. Pujol et al. (2006) assessed myelination from birth to 3 years of age in children's brains using three-dimensional MRI imaging. Intriguingly, this volumetric study demonstrated that a period of rapid myelination started after the 5th month and reached the mature appearance by the 18th month, and the study revealed the relationship with vocabulary acquisition in children. This period of rapid myelination is almost the same as the period in which we observed a rapid change in the present study.

However, contrary to the developmental trajectory for higher scales, the developmental trajectory of MSE for lower scales is diverse across infants, and the reason for this difference remains unclear. Lippe et al. (2009) reported a rapid increase in complexity at lower scales, especially during the early stage of infancy (1–2 months old vs. 2–8 months old) that is followed by a gradual increase. This may suggest the possibility that complexity at lower scales (corresponding to ≤16 ms) saturates by 8 months of age. This may partially explain our finding of a constant complexity value across age after 5 months of age at lower scales. Theoretically, SampEn at finer (i.e., lower) time scales is based on wider frequencies, whereas coarser (i.e., higher) time scales are based on narrower frequencies (i.e., high frequency is filtered out). Signal variabilities in different frequencies must be reflected by differences in time scale. Therefore, a frequency-specific role in the differentiation of cognitive processing (Fries, 2015) and differences in maturational speed (Uhlhaas et al., 2009) may underlie these contradictory findings between the results from high and low time scales.

Surrogate analysis showed a region- and scale-specific increase in surrogate data compared to that in MSE from original data, which may suggest an inherent non-linearity in the MEG data. Furthermore, the developmental trajectory of the spectral power and MSE differed. Specifically, an increase in the power spectral seemed to be prominent after 12 months of age, whereas an

#### REFERENCES

Bolton, P. F., Golding, J., Emond, A., and Steer, C. D. (2012). Autism spectrum disorder and autistic traits in the avon longitudinal study of parents and children: precursors and early signs. J. Am. Acad. Child Adolesc. Psychiatry 51, 249 e25–260 e25. doi: 10.1016/j.jaac.2011.12.009

increase in MSE emerged from early infancy until 15 months of age. Therefore, we assume that the enhancement in MSE with development may be associated with non-linear processes and may be independent of spectral power. In addition, as the outputs of neuronal networks are produced by interactions due to both local dense interconnectivity and sparse long-range excitatory projections (Friston et al., 1995), the resulting dynamics could be expected to operate at multiple scales.

Some potential limitations of the present study must be considered. First, despite frequent MEG recording, we were only able to follow five infants, which precluded statistical evaluation. Second, we did not correct for cognitive behavioral or psychological assessments, which might strengthen our claims. Third, the confounding influence of head motion cannot be excluded from potentially influencing the MSE results. Finally, as a technical consideration, the recent advent of cortical source localization techniques was not applied due to difficulties in performing MRI on infants. Although several limitations must be considered, our findings for the examination of MEG signal variability using MSE may add another dimension to the previously identified neural dynamics of development and may provide useful biomarkers for typically and abnormally developing brains.

#### AUTHOR CONTRIBUTIONS

YY, YM, and MK designed the study. CH, YY, and HK recruited the participants. CH, YY, TI, and DS performed the experiments. CH, SN, and TT analyzed the results and wrote the manuscript. All authors participated in revising the manuscript and approved the final draft of the manuscript.

### FUNDING

This work was supported by JSPS KAKENHI Grant Number JP15K19717 to CH and JP16K10206 to TT and was partially supported by the Center of Innovation Program from the Japan Science and Technology Agency, JST and JST CREST Grant Number JPMJCR17A4, Japan.

### ACKNOWLEDGMENTS

We wish to thank the study participants and their families. We also thank Sachiko Kitagawa for technical assistance and Wasin Hawaree from Fukui University of Technology for data analysis.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer C-KP and the handling Editor declared their shared affiliation.

Copyright © 2018 Hasegawa, Takahashi, Yoshimura, Nobukawa, Ikeda, Saito, Kumazaki, Minabe and Kikuchi. 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.

# Do Complexity Measures of Frontal EEG Distinguish Loss of Consciousness in Geriatric Patients Under Anesthesia?

Sarah L. Eagleman1 \* † , Don A. Vaughn2,3 \* † , David R. Drover1 , Caitlin M. Drover4 , Mark S. Cohen2,5, Nicholas T. Ouellette6 and M. Bruce MacIver1

1 Department of Anesthesiology, Perioperative and Pain Medicine, Stanford University, Palo Alto, CA, United States, 2 UCLA Semel Institute for Neuroscience and Human Behavior, Los Angeles, CA, United States, 3 Department of Psychology, University of Santa Clara, Santa Clara, CA, United States, 4 University of Washington, Seattle, WA, United States, 5 UCLA Departments of Psychiatry, Neurology, Radiology, Psychology, Biomedical Physics and Bioengineering, California Nanosystems Institute, Los Angeles, CA, United States, 6 Department of Civil and Environmental Engineering, Stanford University, Stanford, CA, United States

#### Edited by:

Kay Jann, University of Southern California, United States

#### Reviewed by:

Keiichiro Nishida, Kansai Medical University, Japan Thomas Koenig, Universität Bern, Switzerland

#### \*Correspondence:

Sarah L. Eagleman<EMAIL_ADDRESS>Don A. Vaughn<EMAIL_ADDRESS>†These authors have contributed equally to this work

#### Specialty section:

This article was submitted to Brain Imaging Methods, a section of the journal Frontiers in Neuroscience

Received: 10 June 2018 Accepted: 29 August 2018 Published: 20 September 2018

#### Citation:

Eagleman SL, Vaughn DA, Drover DR, Drover CM, Cohen MS, Ouellette NT and MacIver MB (2018) Do Complexity Measures of Frontal EEG Distinguish Loss of Consciousness in Geriatric Patients Under Anesthesia? Front. Neurosci. 12:645. doi: 10.3389/fnins.2018.00645 While geriatric patients have a high likelihood of requiring anesthesia, they carry an increased risk for adverse cognitive outcomes from its use. Previous work suggests this could be mitigated by better intraoperative monitoring using indexes defined by several processed electroencephalogram (EEG) measures. Unfortunately, inconsistencies between patients and anesthetic agents in current analysis techniques have limited the adoption of EEG as standard of care. In attempts to identify new analyses that discriminate clinically-relevant anesthesia timepoints, we tested 1/f frequency scaling as well as measures of complexity from nonlinear dynamics. Specifically, we tested whether analyses that characterize time-delayed embeddings, correlation dimension (CD), phase-space geometric analysis, and multiscale entropy (MSE) capture loss-ofconsciousness changes in EEG activity. We performed these analyses on EEG activity collected from a traditionally hard-to-monitor patient population: geriatric patients on beta-adrenergic blockade who were anesthetized using a combination of fentanyl and propofol. We compared these analyses to traditional frequency-derived measures to test how well they discriminated EEG states before and after loss of response to verbal stimuli. We found spectral changes similar to those reported previously during loss of response. We also found significant changes in 1/f frequency scaling. Additionally, we found that our phase-space geometric characterization of time-delayed embeddings showed significant differences before and after loss of response, as did measures of MSE. Our results suggest that our new spectral and complexity measures are capable of capturing subtle differences in EEG activity with anesthesia administration—differences which future work may reveal to improve geriatric patient monitoring.

Keywords: anesthesia, geriatric, electrophysiology (EEG), propofol, fentanyl, nonlinear dynamics, multiscale entropy, 1/f

## INTRODUCTION

fnins-12-00645 September 18, 2018 Time: 16:50 # 2

About 50% of geriatric patients aged 65 and older will require anesthesia for a surgical procedure at some time in their remaining years (Kim et al., 2015). It has been suggested that geriatric patients have an increased risk of dementia, delirium, and neurocognitive dysfunction after exposure to anesthetic agents (Avidan and Evers, 2011; Chen et al., 2013; Strøm et al., 2014; Purdon et al., 2015a; Yang and Fuh, 2015). Given this potential risk, several investigators have suggested that geriatric patients may benefit from maintenance at lighter anesthetic levels (Lindholm et al., 2009; Kalkman et al., 2011; Strøm et al., 2014; Petsiti et al., 2015), which may reduce the risks of developing dementia (Chen et al., 2014). However, a downside of light anesthesia is the possibility that a patient may not be unconscious during the procedure. It is estimated that 1–2 out of every 1000 patients undergoing general anesthesia for surgery are insufficiently anesthetized, and yet are immobilized and unable to respond (Bischoff and Rundshagen, 2011). There would be tremendous advantage in the ability to titrate an anesthetic dose accurately to balance an improved medical outcome against the risk of intraoperative awareness.

Electroencephalogram (EEG) signals, especially from the frontal cortex, exhibit stereotypical responses to some anesthetics. For example, loss of consciousness (LOC) correlates with a transition from low amplitude, high frequency EEG waveforms to high amplitude, low frequency patterns, resembling the transition to sleep. In anesthesia, low frequency delta (∼1– 4 Hz) rhythms are replaced gradually by burst suppression patterns as patients transition to deeper surgical planes of anesthesia (Pilge et al., 2014). Since the 1990s, several monitoring devices have been developed to capitalize on these EEG frequency domain transitions (Fahy and Chau, 2018). Unfortunately, the resultant measures are inconsistent as EEG changes depend on the anesthesia and the patient (Akeju et al., 2015; Purdon et al., 2015a,b).

In particular, geriatric patients pose a unique challenge for electrophysiological monitoring because brain activity (as measured with EEG) attenuates with advancing age (Akeju et al., 2015; Purdon et al., 2015b; Lee et al., 2017). This makes it difficult to quantify the frequency changes between awake and anesthetized states accurately in individual patients. None of the EEG intraoperative monitoring methods account for the differences that exist in the geriatric population. In addition, geriatric surgical patients are often medicated with betaadrenergic blockers (beta-blockers). Beta-blockers may hide the cardiovascular signs of inadequate anesthesia (tachycardia and hypertension) (Ghosh et al., 2008), and geriatric patients treated with beta-blockers may require less anesthesia (Zaugg et al., 2003; Ghosh et al., 2008). Importantly, beta-blockers themselves can cause changes in the EEG activity, and thus are likely to cause problems with anesthetic depth monitoring (Johansen, 2001; Zaugg et al., 2003; Ghosh et al., 2008). The result is that EEG monitoring has not been adopted as standard of care. The development of efficacious, anesthetic- and patient-invariant EEG processing techniques could mitigate this and provide a better way to monitor patients and improve outcomes. Thus, in the current study, we analyzed EEG data collected from geriatric patients who received beta-blockers for at least 24 h prior to surgery, and who were subsequently anesthetized with fentanyl and propofol.

Fentanyl and propofol are used in combination routinely to obtain balanced anesthesia for induction during surgical procedures. Fentanyl is a potent opioid that decreases the intensity of response to intubation and provides general pain relief during surgical procedures. Administration of fentanyl correlates with a shift in the EEG pattern from high frequency, low amplitude to low frequency, high amplitude (Scott et al., 1985). Propofol has been shown to increase frontal EEG activity in the alpha (8–14 Hz) and slow (0.1–1 Hz) bands after LOC (Gugino et al., 2001; Feshchenko et al., 2004; Purdon et al., 2013; Akeju et al., 2014). This change in activity occurs in both young and elderly patients (Purdon et al., 2015a); however, the reduction in amplitude with age causes more subtle differences in spectral characteristics with anesthesia onset, which may not be detected by modern EEG monitoring devices (Purdon et al., 2015b).

A potentially more sensitive spectral measure of anesthetic depth might be 1/f. Electrophysiological signals demonstrate 1/flike frequency scaling (He, 2011): the power declines relative to increases in frequency composition. Importantly, this measure has demonstrated sensitivity to the electrophysiological brain state changes associated with sleep (Bédard et al., 2006). Further, 1/f scaling in EEG changes with age (Voytek et al., 2015). To our knowledge, 1/f frequency scaling has not been tested on EEG protocols that include anesthesia.

In addition to spectral measurements, previous studies have demonstrated that complexity measures from nonlinear dynamics correlate with anesthetic depth (Watt and Hameroff, 1988; Widman et al., 2000; van den Broek et al., 2006; Walling and Hicks, 2006; MacIver and Bland, 2014; Eagleman et al., 2018). Structural changes in time-delayed embeddings (attractors) of EEG signals have been reported in both rodents (MacIver and Bland, 2014) and humans (Watt and Hameroff, 1988; Widman et al., 2000; Walling and Hicks, 2006; Eagleman et al., 2018). The awake attractors appear in 3D as spheroids, and then flatten to ellipsoids with LOC (Watt and Hameroff, 1988; Walling and Hicks, 2006; MacIver and Bland, 2014; Eagleman et al., 2018). Previous work has quantified these attractor changes using the correlation dimension (CD) (Grassberger and Procaccia, 1983; Widman et al., 2000; Walling and Hicks, 2006) and a phasespace geometric fit called the ellipse radius ratio (ERR), which fits the attractor with an ellipsoid solid of revolution, and then reports the ratio of the lengths of the minimum and maximum symmetry axes (Eagleman et al., 2018). Another complexity measure, multiscale entropy (MSE), has also been reported to correlate well with existing anesthetic depth measures (Li et al., 2010; Liu et al., 2015). In the current study, we tested whether these complexity measures could distinguish between subtle changes in EEG activity occurring before and after LOC.

We performed a retrospective analysis of data collected from geriatric patients (aged 65 and older) who were on betaadrenergic blockers at least 24 h prior to surgery. Patients were anesthetized with fentanyl and propofol. We identified and analyzed 20 s clips of frontal EEG from before and after loss of response to verbal commands (LOR, considered here as loss of consciousness). We compared frequency-derived measures, as well as complexity measures, of the EEG signal in the before-andafter clips to identify measures that discriminated between the two states.

### MATERIALS AND METHODS

fnins-12-00645 September 18, 2018 Time: 16:50 # 3

#### Study Protocol

All procedures took place under an approved protocol from the Stanford School of Medicine Administrative Panel on Human Subjects in Medical Research (ClinicalTrials.gov, NCT00938782). The database consisted of 67 surgical patients all older than 65 years of age. Patients underwent anesthesia for non-cardiac procedures classified as status 1–3 by the American Society of Anesthesiologists. All patients in the original study were receiving beta-adrenergic treatment for a minimum of 24 h preoperatively and received their medication prior to surgery. We reviewed the database retrospectively and analyzed data under a separate Stanford-approved IRB. Out of the 67 patients, 28 patients were used for analysis (details of selection below). These 28 patients had an average age of 76 (± 6) years and included 18 males and 10 females. They had an average body mass index (BMI) of 27.3 (± 4.7). The healthy BMI range for patients aged 65 and older has been suggested as 23 to 30 (Porter Starr and Bales, 2015). Our analysis included 8 patients with a BMI greater than 30 and 5 patients with a BMI less than 23.

The surgical procedures and anesthetic administration have been described in detail previously (Drover et al., 2011). Briefly, patients were induced with fentanyl (1–3 mcg/kg), propofol (1–2 mg/kg), and muscle relaxant (if required) using either rocuronium (0–1 mg/kg) or vecuronium (0.1 mg/kg). After intubation, sevoflurane in oxygen with 50–60% nitrous oxide was initiated and used to maintain anesthesia. Medical staff recorded significant clinical events including induction, anesthetic administration and dosage, LOR to verbal commands, and time.

The average time of fentanyl use prior to LOR was 2.1 ± 1.6 min. With propofol, the average before LOR was 1.0 ± 0.6 min. A total of 6 patients received muscle relaxant before LOR. One patient received vecuronium 0.7 min before LOR, and 5 patients received rocuronium an average of 5.2 min before LOR. A total of 14 patients received muscle relaxant within a 4 min window following LOR. Three patients received vecuronium and 11 received rocuronium, both an average of 0.8 min following LOR. In addition, 18 patients were started on sevoflurane in oxygen with 50–60% nitrous oxide an average of 1.7 min following LOR.

#### EEG Recording and Preprocessing

Electroencephalogram recordings were acquired using a SedLine Legacy EEG monitor (Masimo, Irvine, CA, United States). The manufacturer's standard adhesive electrode was attached to the patient prior to starting the anesthetic, as per the manufacturer's instructions. EEGs were recorded at approximately F7 grounded to Fpz and referenced to ∼1 cm above Fpz (**Figure 1A**). Data was digitized at 250 Hz. Records of the surgical events including time stamps of start of induction, LOR to verbal stimuli, and administration and dosages were de-identified and then used for analysis. The EEG recordings were de-trended, and notch filtered using a second-order Butterworth Infinite Impulse Response (IIR) filter to remove 60 and 120 Hz noise prior to analysis. Additionally, a second-order Butterworth IIR filter was used to remove a 78.125 Hz impedance measurement pulse generated by the EEG monitor system.

We identified each subject's 20 s pre-LOR and post-LOR EEG clip via a two-step algorithm which was a mix of exclusion and inclusion criteria (**Figure 1A**). First, three of the authors (SLE, DRD, and MBM) visually inspected the EEG traces, spectrums, processed spectrograms to identify EEG periods containing burst suppression or artifacts (exclusion criteria). Second, given the fundamentally imprecise metric of LOR, from these remaining clips we used the clips most temporally-distant from the LOR timepoint (within a 2 min window) to obtain the clips most representative of pre-LOR and post-LOR (inclusion criteria). From the original 67 patients, 28 were selected who had LOR timestamps and artifact and noise free EEG clips for at least 20 continuous seconds before and after LOR.

#### Spectral Analyses

To visualize the spectral changes that occurred in our clips before and after LOR, we performed multitaper spectral analysis on the 20 s clips during pre and post LOR period using the MATLAB Chronux toolbox (**Figure 1B**; Mitra and Bokil, 2008).1 Specifically, we used a time-bandwidth product of 5 with 9 tapers, limited the frequency ranges calculated to 0 to 50 Hz, and computed the theoretical error range at the 95% confidence interval. Power values were expressed in decibels.

To give an example of the temporal profile of the spectral changes that occurred during the windows surrounding the LOR transitions, we computed a normalized spectrogram (**Figure 1C**). We calculated the Fourier transform using Hann windows with half window overlaps. We then cutoff the frequencies above 50 Hz, converted the magnitude to decibels (dB), and scaled the spectrogram output by its maximum magnitude.

We calculated the spectral edge frequency (i.e., the frequency bounding 95% of the power from above) and total power using multitaper spectral analysis (Mitra and Bokil, 2008) without limiting the frequency range to below 50 Hz. We calculated the percentage of total power for individual frequency bands per condition. The percentage of total power was used because of prior reports of significant changes in total power with exposure to anesthetics. The ranges we used for the frequency bands were as follows: delta: 0.1 to 4 Hz; theta: 4 to 8 Hz; alpha: 8 to 14 Hz; beta: 14 to 30 Hz; and gamma: above 30 Hz (Gugino et al., 2001; Purdon et al., 2013). To test whether we observed similar changes in before and after clips compared to other studies, we calculated the slow frequency component (0.1–1 Hz) separately.

To determine whether 1/f characteristics change before and after LOR, we fit each patient's spectral power to c/f α , where α and

1http://chronux.org/

FIGURE 1 | Characterizing loss of response (LOR) using electroencephalogram (EEG) analysis. (A) EEGs were recorded from F7 during induction with propofol and fentanyl anesthesia in geriatric patients receiving beta-adrenergic blockade. We selected continuous artifact-free 20 s clips from before (red) and after (blue) patients lost response to verbal commands. (B) We analyzed these clips using multitaper spectral analysis. In this patient (PT09) an increase in lower frequencies and a decrease in higher frequencies from before to after LOR can be observed. (C) To view the dynamics of the spectral changes through time we plotted a normalized spectrogram starting 2 min before LOR to 2 min after LOR. We can see the increase in lower frequencies and decrease in higher frequencies coordinates well with the LOR timestamp. (D) We tested several complexity measures of time-delayed embeddings (attractors). An awake attractor (red) looks less ellipsoidal than the anesthetized attractor (blue). (E) One method we used to quantify change in attractors before and after LOR is by fitting the attractor with an ellipsoid solid of revolution. We then calculate the ratio of the minimum and maximum radii.

c were free parameters representing the quickness of frequency decay and an arbitrary constant, respectively. We determined α and c for each patient by minimizing the L2 norm of the residuals between predicted and actual values.

#### Multiscale Entropy

fnins-12-00645 September 18, 2018 Time: 16:50 # 5

As complex signals often have meaningful relationships at multiple timescales, we used MSE to characterize the relationship and complexity of the EEG time series. MSE utilizes an algorithm that calculates a traditional entropy metric at several timescales (Costa et al., 2002). In this case, we used sample entropy as our multiscale metric (Pincus, 1991) as it has been shown to extract meaningful complexity differences in a variety of physiological signals (Costa et al., 2003; Bian et al., 2009; Liu et al., 2013). We calculated MSE using custom MATLAB (The MathWorks, Inc.) scripts.

## Characterization of Dynamical Attractors

#### Correlation Dimension

We constructed three-dimensional time-delayed embeddings (attractors) of the EEG signal before and after LOR using an 8 ms delay (**Figure 1D**) as previously described (Watt and Hameroff, 1988; Walling and Hicks, 2006; MacIver and Bland, 2014; Eagleman et al., 2018). We chose this delay as we observed shape changes in attractors when plotted at this timescale. We tested whether the CD captured this shape change with LOR. CD is used to determine the non-integer (fractal) dimensionality of irregular objects (e.g., a point cloud in our case) as described previously (Grassberger and Procaccia, 1983; Widman et al., 2000; Walling and Hicks, 2006). We also tested whether significant differences in CD could be observed when we increased the dimensionality of the embeddings from 3D to 5D.

We tested whether the embedding delay of the attractor impacted the estimate of the CD value. To estimate the optimal delay for creating the attractors we calculated the first zerocrossing of the autocorrelations of the pre-LOR and post-LOR signals. This is the first time lag that the EEG signal differs maximally with itself. We used this value to set our maximum range of embedding delays and tested multiple delays between the shortest time window (shifting the EEG by 1 point) to the largest (set by the autocorrelation zero crossing). We tested delays of 4, 8, 12, 52, 100, 500, 1000, 1500, 2000, and 2500 ms. We tested these delays for both 3D and 5D time-delayed embeddings. We tested whether a trend existed between embedding delay and CD using the Spearman correlation.

#### Ellipse Radius Ratio

We constructed three-dimensional time-delayed embeddings, as described above (**Figure 1D**). We quantified this shape change by fitting the three-dimensional attractor to an ellipsoidal solid of revolution (**Figure 1E**; Khachiyan, 1980; Eagleman et al., 2018). The lengths of the symmetry axes of the ellipsoid were calculated and the ratio of the minimum and maximum axes (which we term the ellipsoid radius ratio, ERR) was used to quantify the shape change. A radius ratio of 1 implies a sphere, while smaller ratios imply more strongly ellipsoidal shapes.

Similar to the CD, we tested also whether the ERR was changed by the embedding delay time. We created time delayed embeddings using the same delays as were used to calculate the CD. We tested whether a trend existed between embedding delay and ERR using the Spearman correlation statistic.

#### Correlations Between EEG Measures and Effect Size of EEG Measures

To test whether our EEG measures correlated with patient age or body-mass index (BMI), we calculated the Spearman correlation statistic between the change in ERR (LORpost–LORpre) at the shortest delay (4 ms) with corresponding patient ages and BMIs. We did the same for the MSE results.

To test whether our EEG measures correlated with each other or with spectral changes, we calculated the Spearman correlation statistic between the changes in these measures before and after LOR. Specifically, we tested the correlation between changes in ERR and MSE as well as between ERR or MSE and percentage of power change (LORpost–LORpre) in the individual frequency bands that we measured (delta, theta, alpha, beta, and gamma). For the change in ERR, we chose the values calculated at the shortest delay (4 ms) as these showed the most significant changes before and after LOR.

We calculated a paired-data Cohen's D on our EEG measures before and after LOR for the novel spectral (1/f) and complexity (MSE and ERR phase space analysis) measures that showed significant differences before and after LOR. We also calculated Cohen's D for the percentage of power in each of the frequency bands (delta, theta, alpha, beta, and gamma) for comparison.

#### Statistics

We corrected significance values for multiple univariate statistical comparisons within a particular analysis type, by using the Holm-Bonferroni method—a sequentially-rejective procedure (Holm, 1979). Specifically, we corrected p values for pre vs post metrics within each of the following analyses: (1) the power percentage across all 5 frequency bands (delta, theta, alpha, beta, gamma), (2) the CD across all 10 embedding delays, (3) the ERR across all 10 embedding delays, and (4) the correlation of MSE, ERR, and frequency band power. We report our results as medians (25, 75 percentiles), and significance values (p) are calculated from Wilcoxon Signed Rank Tests.

#### RESULTS

#### Spectral Analyses

We performed multitaper spectral analysis to quantify the changes that occur before and after LOR (**Figure 2A**). To compare our measures to previous reports, we computed a common spectral measure that correlates with anesthetic depth: the spectral edge frequency. We also tested whether total power differed from before and after LOR. We did not observe a significant change in spectral edge frequency before and after LOR (**Figure 2B**, pre-LOR 14.1 Hz [9.6, 19.4], post-LOR 13.0 Hz

[11.1, 17.6], p = 0.64). We did not observe a significant difference in total power before and after LOR (pre-LOR 35.2 dB [33.0, 36.2], post-LOR 34.6 dB [32.8, 36.3], p = 0.21).

To determine the changes in frequency bands, we calculated the changes in delta (0.1–4 Hz), theta (4–8 Hz), alpha (8– 14 Hz), beta (14–30 Hz), and gamma (above 30 Hz) ranges. Significant differences were found for alpha (**Figure 2C**) and gamma (**Figure 2D**), and all spectral results are summarized in **Table 1**. Similar to previous reports, we observed a significant increase in alpha from before to after LOR. In addition, we observed a significant decrease in gamma power from before to after LOR. In a separate analysis, we separated out the slow (0.1–1 Hz) frequency component from the EEG signal to see if the percentage of slow activity changed before to after LOR. We did not observe a significant difference in this frequency band (pre-LOR 48.6% [37.6, 62.9], post-LOR 54.3% [34.0, 64.4], p = 0.57).

To determine whether 1/f characteristics change before and after LOR, we fit each patient's spectral power to c/f α (**Figure 3A**). Overall, 71% of the patients showed an increase in the value of α from pre-LOR to post-LOR (**Figure 3B**). This difference in α before and after LOR differed significantly from the null hypothesis of no change (median change = 0.17, p < 10−3 , **Figure 3C**).

#### Characterization of Dynamical Attractors

We began by testing whether significant differences could be observed in CD between 3D attractors plotted at an 8 ms embedding delay. We observed a similar flattening of the attractor and more ellipsoidal shapes after LOR (**Figure 4A**). We quantified the pre- to post-LOR changes in 3D attractors using CD. We did not observe a significant difference using this measure (**Figure 4B**). We also calculated the CD for 5 dimensional attractors. Here again we did

Percentage of power in individual frequency bands are reported as medians [25, 75 percentiles]. Significance values reported are from Wilcoxon signed rank tests, corrected for multiple comparisons. The ∗ indicates the percentage of power change was significant at p < 0.05 corrected.

not observe a significant difference pre- and post-LOR (**Figure 4C**).

We then tested our ERR with the 8 ms attractors before and after LOR. We found a significant difference between pre- and post-LOR (**Figure 4D**, p = 0.04 corrected).

To test the impact of the embedding delay on the attractor calculations we chose multiple embeddings delays between the smallest possible delay (4 ms, shifting the EEG by 1 point) and the largest (set by the first zero-crossing of the autocorrelation). We calculated the first zero-crossing for both the pre- and post-LOR period to see if there was a difference. We did not find a significant difference between the pre- and post-LOR values (pre-LOR 1800 ms [1392, 2856], post-LOR 2176 ms [1664, 2480], p = 0.35 corrected). Given these results, we decided to calculate CD using the following embedding delays: 4, 8, 12, 52, 100, 500, 1000, 1500, 2000, and 2500 ms. The impact of embedding delay on attractor shape is shown in **Figure 5A**. We did not observe a significant difference at any delay for 3D or 5D CD calculations, nor did we observe a significant correlation between CD and embedding delay (**Figure 5B**).

We performed the same test of different embedding delays with our ERR analysis. The ERR showed a significant positive relationship with the difference in post-LOR–pre-LOR conditions (**Figure 5C**, p = 0.04, Spearman correlation). This relationship likely was driven by a reduction in the ERR between the pre- and post-LOR states at short embedding delays: only embedding delays of only 4 and 8 ms showed a significant difference (p4 = 0.017 corrected, p8 = 0.04 corrected).

#### Multiscale Entropy

A MSE analysis between the pre-LOR and post-LOR conditions revealed a scale-dependent change in sample entropy (**Figure 6A**). Complexity decreased at short scale factors, showing a statistically significant trend toward increasing at medium scale factors (rspearman = 0.45, p < 10−3 percentile permutation test). The difference in sample entropy between post and pre conditions then showed a decreasing relationship with a further increase in scale factor (rspearman = -0.31, p < 10−3 percentile permutation test). The initial decrease in complexity

at a scale factor of 1, which was not significantly different by median (p = 0.48), likely resulted from sharp decreases in MSE among several of the participants (**Figure 6B**).

#### Correlations Between EEG Measures and Effect Size of EEG Measures

observed that 71% of the patients show an increase in α after LOR.

Spearman correlation revealed no significant relationship in our ERR or MSE changes before and after LOR and with patient age or BMI (**Table 2**). Spearman correlation between our ERR (calculated at the shortest delay, 4 ms) and MSE changes before and after LOR revealed that they are correlated with each other (**Table 3**, r = 0.54, p = 0.0034 corrected). Additionally, a significant correlation between ERR change and change in percentage of gamma activity before and after LOR was observed (r = 0.67, p = 0.00012 corrected). Significant correlations are also observed between changes in MSE and changes in delta (r = -0.60, p = 0.00084 corrected), alpha (r = 0.57, p = 0.0019 corrected), beta (r = 0.74, p = 10−5 corrected) and gamma power (r = 0.64, p = 0.00034 corrected). A summary of all of the measured values before and after LOR can be found in **Table 3**.

for 6 patients. We quantified the difference in attractors using both correlation dimension and our ellipse radius ratio (ERR). (B) We did not observe a significant difference in correlation dimension in 3 dimensions, nor did we observe a significant difference in correlation dimension in 5 dimensions (C). (D) We did observe a significant difference in our ERR measure from before to after LOR (p = 0.04 corrected). The ∗ indicates significant difference between pre-LOR and post-LOR measures at p < 0.05 corrected.

Cohen's D values revealed medium effect sizes for both complexity measures: sample entropy and ERR (**Table 4**, DMSE = −0.55, DERR = −0.69). Cohen's D values were also sizeable for percentage of power differences before and after LOR in alpha and gamma frequency bands (**Table 4**, D1/f = 0.68, DAlpha = 0.61, DGamma = −0.55). The rest of the frequency bands had small effect sizes (**Table 4**).

#### DISCUSSION

We found spectral results similar to those reported for propofol anesthesia (Gugino et al., 2001; John et al., 2001; Purdon et al., 2013, 2015a). Specifically, we saw increases in the percentage of alpha activity and decreases in the percentage of gamma activity with LOR. Conversely, we did not observe a significant difference in the percentage of slow (0.1–1 Hz) frequency from before to after LOR, nor did we observe a significant difference in spectral edge frequency before and after LOR. This is likely due to the overall reduced amplitudes of EEG signals in our elderly patients (Purdon et al., 2015b), which reduces the magnitude of the spectral changes. In addition, patients were sedated heavily before they lost consciousness, so EEG changes were subtle. Since we were doing a retrospective data analysis, we were not able to control the initiation of maintenance anesthesia, so some of our patients were on sevoflurane in oxygen with 50– 60% nitrous during the post-LOR timepoints. This may have obscured a change in spectral edge frequency as nitrous oxide

FIGURE 5 | Attractor shape changes are embedding delay dependent. (A) Attractors from PT02 shown at 4, 52, 1000, and 2500 ms. The shape change in the attractor is not observed at higher embedding delays. We tested whether our attractor characterization analyses correlated with embedding delay. We computed correlation dimension and ellipse radius ratio (ERR) at 10 different delays from 4 to 2500 ms. (B) No significant changes in correlation dimension were observed from before to after LOR for any embedding delay nor was correlation dimension correlated with the embedding delay. Here we show results from the 3D correlation dimension calculation. (C) We observed a significant correlation between the ERR and embedding delay (rspearman = 0.17, p = 0.04) and a significant difference in the ERR between pre- and post-LOR for 4 ms and 8 ms delays (p4 = 0.017 corrected, p8 = 0.04 corrected). The ∗ indicates a significant difference between the measured values before and after loss of response at that embedding delay at p < 0.05 corrected.

maintains higher frequencies (Rampil et al., 1998). Sevoflurane, on the other hand, causes similar changes in EEG activity compared to propofol (Akeju et al., 2014; Purdon et al., 2015a). Overall, these nuances between common clinical anesthetics highlight the importance of developing new tools to better distinguish anesthetic states using EEG (Eagleman and Drover, 2018).

One interesting spectral analysis not previously applied to anesthesia EEG is the calculation of the 1/f frequency scaling. We chose this measure as it distinguishes different brain states

TABLE 2 | Summary of correlations between observed changes before and after loss of response (LOR) and patient demographics.

The Spearman correlation between changes in the ellipse radius ratio (ERR, calculated from the shortest delay 4 ms) and multi-scale entropy (MSE) values before and after LOR (post – pre) with age and body-mass index (BMI) of patients used in the analysis is reported here. No correlations were significant.

TABLE 3 | Summary of correlations between measured changes before and after loss of response (LOR).

The Spearman correlation between changes in the ellipse radius ratio (ERR) and multi-scale entropy (MSE) values before and after LOR (post-pre) is reported here. Additionally, Spearman correlation between these values and the changes in percentage of power in individual frequency bands (delta, theta, alpha, beta, and gamma) before and after LOR (post–pre) are reported as well. The ∗ indicates the correlation between the two parameters was significant at p < 0.05 corrected.

TABLE 4 | Summary of Cohen's D for measured values before and after loss of response (LOR).

We calculated Cohen's D (sign convention: post-pre) to quantify the effect sizes of our EEG measures. We included all measures that showed significant differences before and after LOR, and the standard frequency bands for comparison. ERR value shown here is for the 4 ms delay.

(Bédard et al., 2006) and ages (Voytek et al., 2015). To our knowledge, we have demonstrated the first observation of a change in 1/f frequency scaling in EEGs in an anesthesia protocol. 1/f frequency scaling was sensitive to before and after LOR.

Previous studies reported differences in CD with anesthetic depth (Widman et al., 2000; Walling and Hicks, 2006), but our study differs in several important aspects. In these earlier studies, patients were anesthetized with sevoflurane to deep levels of anesthesia, but in our study, patients were anesthetized with propofol and fentanyl; further, we limited our analysis to 20 s before and after LOR. Additionally, we used a unique patient population consisting only of geriatric patients with betaadrenergic blockade. We determined that our CD results were not based on embedding delay or differences between 3 and 5 dimensions.

There are several possibilities as to why we did not find significant differences in CD before and after LOR. For instance, the attractor might be better resolved with higher sampling frequency. Additionally, brain activity might be better represented in an even higher dimensionality embedding. Changes in CD also may not be observable due to agerelated changes in complexity (Pierce et al., 2000; Müller and Lindenberger, 2012; Sleimen-Malkoun et al., 2014). The other possibility is that CD might not be sufficiently sensitive to detect the changes that occur during before and after LOR.

We did observe a change in attractor shape similar to what has been previously described (Watt and Hameroff, 1988; Walling and Hicks, 2006; MacIver and Bland, 2014; Eagleman et al., 2018). Additionally, our phase-space analysis based on the geometry of attractors showed significant differences before and after LOR at very short delays (shifting the EEG signal by 1 or 2 points), as previously reported (Eagleman et al., 2018). However, our results demonstrate variability in our population. To determine the source of this variation, better control and measurement of the anesthetics administered in a prospective study is needed. We noticed the magnitude of our results were reduced in the current study compared to our recently published results (Eagleman et al., 2018). This may be due to the current patient population being more sensitive to anesthetics and thus more sedated at our pre-LOR timepoint. We also tested the impact of the embedding delay on our analysis to explore whether changes in EEG signals before and after LOR existed at longer timescales. We found that shorter delays better distinguished before and after LOR and results were no longer significant when the signal was delayed by three points. Visual inspection of the attractor shapes supports this result. Further work is needed to elucidate whether calculations performed in real-time can classify anesthetic depth adequately.

In addition, we tested whether another complexity measure, MSE, could distinguish before- and after-LOR timepoints. We noticed that MSE values appear to converge within a smaller range in the post-LOR period. The distribution suggests that there is one group of patients that show a large decrease in this EEG measure with LOR, while another group does not. The differences between the groups were outside the scope of the current study; however, we plan to examine these differences in the future in a well-controlled prospective study. Additionally, previous reports that showed high correlation of MSE with existing anesthetic depth measures (BIS index and expert anesthesiologist assessment) used EEG data from the entirety of the surgery (induction through recovery), which exposed patients to deeper anesthetic levels than in our study

(Liu et al., 2015). Thus, further testing of this analysis on more clinically relevant timepoints is needed.

Since several EEG measures were applied here to an anesthetic dataset including some that have not been explored much previously (ERR), we calculated the Spearman correlation between the changes in EEG measures before and after LOR. We did this for the shortest delay of ERR, MSE, and the percentage of all the individual frequency bands we included. We found a significant correlation with the ERR and MSE values. This is interesting because it suggests that the ERR change, and thus attractor shape change, may be related to changes in the complexity of the signal (revealed by reduction in MSE with anesthesia onset). Additionally, a significant correlation between ERR and MSE with the percentage of gamma change was observed. These relationships are expected as both reductions in gamma activity and reductions in entropy with anesthetic administration have been previously reported (John et al., 2001; Li et al., 2010). The change in MSE is also correlated with changes in the percentages of delta, alpha, and beta power. This again is expected given the ability of MSE to capture all of the spectral changes observed with anesthesia onset (Li et al., 2010).

We tested the effect size of the EEG measures using Cohen's D. None of these measures had large effect sizes (> 1) indicating the challenge with detecting subtle EEG changes in geriatric patients. ERR and 1/f frequency scaling had the largest effect sizes, which were medium in magnitude. This indicates that supplementing complexity measures may improve geriatric patient monitoring. However, further testing of these analyses on more clinically relevant time points and on full EEG traces is needed to test this idea.

As with any study involving EEG, muscle contamination is a potentially important source of artifact. Three of the authors visually inspected all of the EEG clips as well as the processed EEG spectrum and spectrograms to ensure they were free of artifacts. Activity from facial and neck muscles can appear 20 Hz and above, and thus into the frequency ranges we used for analyses (Shackman et al., 2009; Claus et al., 2012; Muthukumaraswamy, 2013); however, as we have reasoned previously (Eagleman et al., 2018), it is important not to throw out higher frequency activity, as it plays an important role in brain-state dynamics (Muthukumaraswamy, 2013), especially in judging anesthetic depth (Sleigh et al., 2001). The measures that we have used have been tested in similar experimental paradigms on intracranial recordings free from EMG contamination (Bédard et al., 2006; MacIver and Bland, 2014; Voytek et al., 2015).

Since these results are only from retrospective analyses our work is limited in several ways. We were not able to control the timing of drug delivery (such as delivery of muscle relaxants); several anesthesiologists with potentially diverse clinical practices were involved, and consciousness measures were restricted to the first loss of response to verbal commands. However, all participating anesthesiologists were instructed to administer anesthesia and a small number of adjuvant agents as per strict protocol guidelines (Drover et al., 2011). Additionally, the protocol of anesthetic and adjuvant agent administration in this retrospective dataset is aligned with current clinical practice. Thus, our results are relevant to current practices of balanced anesthesia administration. Additionally, we tested whether any of our measured results were correlated with patient age or BMI and found no significant correlations. Future prospective work will include several measures to better titrate our analysis to anesthetic action. Collection of blood samples or exhaled vapor can help us correlate results with anesthetic delivery more accurately. Whenever possible, future work should control the delivery of muscle relaxants and the initiation of maintenance anesthetics to separate out the effects of individual anesthetic and adjuvant agents on our measures. Additionally, our results need to be tested on more clinically relevant timepoints, and alongside spectral measures on full EEG traces instead of clips, to better prepare our analyses for clinical application.

Nonetheless, we have observed significant differences before and after LOR using several techniques in a traditionally hardto-monitor patient group. Future work will discern if these results are useful supplemental tools to better guide physicians in monitoring anesthetic depth in sensitive patient populations. Development of better EEG analysis techniques will hopefully encourage the wide adoption of EEG monitoring and improve the standard of care.

### CONCLUSION

We found that frontal spectral changes before and after LOR in geriatric patients were limited to the alpha and gamma ranges. Further, we showed that 1/f frequency scaling differed before and after LOR. We tested the ability of several measures from nonlinear dynamics, including CD, MSE, and a geometric characterization of time-delayed embeddings, to distinguish LOR timepoints. Among these, MSE and the geometric characterization showed significant differences and had comparable or greater effect sizes to standard frequency measures. In the future, these results may enable the development of better methods of quantifying anesthetic depth in geriatric patients as they are able to significantly discriminate between the subtle EEG changes that occur before and after loss of response.

### AUTHOR CONTRIBUTIONS

MM, DD, and SE conceived the study. DD acquired the data. SE, DD, and CD organized anesthesia records for analysis. SE and DV analyzed the data. MM, NO, and MC supervised the data analysis. SE and DV wrote the manuscript. All authors contributed to the intellectual content with rounds of review and approved the final version of the manuscript for publication.

### FUNDING

This work was supported by a seed grant from the Department of Anesthesiology, Perioperative and Pain Medicine at Stanford University (SE, CD, DD, NO, and MM), the Anesthesia Training Grant in Biomedical Research (NIH T32 GM 089626-09, SE), and by the UCLA QCBio Collaboratory Postdoctoral Fellowship directed by Matteo Pellegrini (DV and MC).

#### REFERENCES

fnins-12-00645 September 18, 2018 Time: 16:50 # 12

consciousness from propofol. Proc. Natl. Acad. Sci. U.S.A. 110, E1142–E1151. doi: 10.1073/pnas.1221180110

**Conflict of Interest Statement:** DD is a consultant for Masimo Inc. The original study was conducted with financial support from Hospira, Inc.

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

Copyright © 2018 Eagleman, Vaughn, Drover, Drover, Cohen, Ouellette and MacIver. 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.

# Dynamic Complexity of Spontaneous BOLD Activity in Alzheimer's Disease and Mild Cognitive Impairment Using Multiscale Entropy Analysis

#### Edited by:

Danny J. J. Wang, University of Southern California, United States

#### Reviewed by:

Kay Jann, University of Southern California, United States Albert Yang, Harvard Medical School, United States Robert X. Smith, Washington University in St. Louis, United States Hengyi Rao, University of Pennsylvania, United States

\*Correspondence: Jie Xiang<EMAIL_ADDRESS>
†These authors are co-first authors

#### Specialty section:

This article was submitted to Brain Imaging Methods, a section of the journal Frontiers in Neuroscience

Received: 27 April 2018 Accepted: 07 September 2018 Published: 01 October 2018

#### Citation:

Niu Y, Wang B, Zhou M, Xue J, Shapour H, Cao R, Cui X, Wu J and Xiang J (2018) Dynamic Complexity of Spontaneous BOLD Activity in Alzheimer's Disease and Mild Cognitive Impairment Using Multiscale Entropy Analysis. Front. Neurosci. 12:677. doi: 10.3389/fnins.2018.00677 Yan Niu1† , Bin Wang1,2† , Mengni Zhou1 , Jiayue Xue1 , Habib Shapour1 , Rui Cao1 , Xiaohong Cui1 , Jinglong Wu3,4 and Jie Xiang1 \*

1 College of Information and Computer, Taiyuan University of Technology, Taiyuan, China, 2 Department of Radiology, First Hospital of Shanxi Medical University, Taiyuan, China, 3 Key Laboratory of Biomimetic Robots and Systems, Ministry of Education, Beijing Institute of Technology, Beijing, China, 4 Graduate School of Natural Science and Technology, Okayama University, Okayama, Japan

Alzheimer's disease (AD) is characterized by progressive deterioration of brain function among elderly people. Studies revealed aberrant correlations in spontaneous blood oxygen level-dependent (BOLD) signals in resting-state functional magnetic resonance imaging (rs-fMRI) over a wide range of temporal scales. However, the study of the temporal dynamics of BOLD signals in subjects with AD and mild cognitive impairment (MCI) remains largely unexplored. Multiscale entropy (MSE) analysis is a method for estimating the complexity of finite time series over multiple time scales. In this research, we applied MSE analysis to investigate the abnormal complexity of BOLD signals using the rs-fMRI data from the Alzheimer's disease neuroimaging initiative (ADNI) database. There were 30 normal controls (NCs), 33 early MCI (EMCI), 32 late MCI (LMCI), and 29 AD patients. Following preprocessing of the BOLD signals, whole-brain MSE maps across six time scales were generated using the Complexity Toolbox. One-way analysis of variance (ANOVA) analysis on the MSE maps of four groups revealed significant differences in the thalamus, insula, lingual gyrus and inferior occipital gyrus, superior frontal gyrus and olfactory cortex, supramarginal gyrus, superior temporal gyrus, and middle temporal gyrus on multiple time scales. Compared with the NC group, MCI and AD patients had significant reductions in the complexity of BOLD signals and AD patients demonstrated lower complexity than that of the MCI subjects. Additionally, the complexity of BOLD signals from the regions of interest (ROIs) was found to be significantly associated with cognitive decline in patient groups on multiple time scales. Consequently, the complexity or MSE of BOLD signals may provide an imaging biomarker of cognitive impairments in MCI and AD.

Keywords: multiscale entropy, Alzheimer's disease, mild cognitive impairment, blood oxygen level-dependent signals, dynamic complexity

#### INTRODUCTION

fnins-12-00677 September 27, 2018 Time: 19:7 # 2

Functional connectivity (FC) of spontaneous blood oxygen leveldependent (BOLD) signals in functional magnetic resonance imaging (fMRI) has become an important tool for probing brain function changes in normal aging and neurodegenerative diseases. However, relatively few studies have investigated the temporal dynamics of BOLD signals and its relations with pathologic changes in neurophysiology (Sporns et al., 2000; Friston et al., 2003; Wu et al., 2012). As the most complex organ of the human body, the human brain regulates multifaceted actions with billions of neurons and synapses (Fox et al., 2007). Therefore, the BOLD signals possess complex temporal fluctuations, which could be imitated by nonlinear dynamical processes (Soltysik et al., 2004; Stephan et al., 2008; Yan et al., 2017).

Over the past few years, several statistical methods have been applied to quantify the temporal dynamics of physiological systems. A widely used non-linear statistical method is sample entropy (SE) proposed by Richman and Moorman (2000). SE improved approximate entropy (ApEn) proposed by Pincus (1991), by resolving the problem of erratic results due to vector self-matching. Many studies evidenced the effectiveness of SE in the complexity analysis of time series data of biological systems (Sokunbi et al., 2013, 2014). However, recent studies found that neural signals in the brain possess correlations over a wide range of temporal and spatial scales, stemming from long-range interactions (Costa et al., 2005; Peng et al., 2009; Morabito et al., 2012). Therefore, SE may not be adequate to fully capture the complexity of neural signals by only calculating signal entropy on a single scale.

The multiscale entropy (MSE) was proposed (Costa et al., 2002) to investigate the dynamic complexity of a time series across multiple temporal scales. Several studies have demonstrated the efficacy of MSE for quantifying the complexity of BOLD signals in aging (Yang et al., 2013; Smith et al., 2014). Yang et al. (2013) employed MSE analysis to investigate the complexity of BOLD signals between the younger and older groups, and found significant decreases in MSE in older subjects. Smith et al. (2014) explored the effect of healthy aging on the entropy of resting-state fMRI (rs-fMRI) using MSE analysis, and the results revealed enhanced contrast between healthy young and aged volunteers at longer time scales. However, the dynamic complexity of BOLD signals in neurodegenerative diseases across multiple temporal scales remains largely unexplored.

Alzheimer's disease (AD) is a neurodegenerative disease characterized by progressive deterioration of cognitive and behavioral function (Ballard et al., 2011). Mild cognitive impairment (MCI) is a neurological disorder occurring before the onset of early AD as an intermediate stage at a high risk of developing AD (Petersen et al., 1999; Belleville et al., 2011). A few studies found decreased complexity of BOLD signals in AD by using single-scale entropy analysis (Liu et al., 2013; Wang et al., 2017). However, the complexity alterations of BOLD signals in MCI and AD patients across multiple time scales remain unclear.

We obtained BOLD rs-fMRI data from the Alzheimer's disease neuroimaging initiative (ADNI1 ) database, including 30 normal control (NC), 33 early MCI (EMCI), 32 late MCI (LMCI), and 29 AD subjects. MSE maps of the four groups across multiple time scales were calculated and the clusters with significant MSE differences were identified. We then examined the relationships between MSE values and scores of cognitive assessments on all time scales. Finally, we investigated the relationship between MSE and gray matter volume (GMV) on all time scales.

#### MATERIALS AND METHODS

#### Participants

A total of 124 subjects were selected from ADNI-2 database, including 30 NC subjects (aged 74.18 ± 5.96 years; 19 females; education: 16.8 ± 2.0 years), 33 EMCI subjects (aged 72.01 ± 5.87 years; 16 females; education: 15.5 ± 2.4 years), 32 LMCI subjects (aged 72.57 ± 8.16 years; 13 females; education: 16.5 ± 2.1 years), and 29 AD subjects (aged 72.33 ± 7.26 years; 18 females; education: 16 ± 2.7 years). For each subject, there were cognitive assessments including Mini-Mental State Examination (MMSE), Clinical Dementia Rating (CDR), and Functional Activities Questionnaire (FAQ). **Table 1** summarizes the demographic and clinical characteristics of the participants.

#### Data Acquisition and Data Processing

All subjects went through resting-state BOLD fMRI scans with their eyes closed on a 3.0 T scanner (Philips Medical Systems)

1http://adni.loni.usc.edu/

Data are given as the mean ± standard deviation (SD).

The p-values were obtained by one-way ANOVA.

using the following parameters: repetition time (TR) = 3000 ms; echo time (TE) = 30 ms; slice thickness = 3.3 mm; flip angle = 80◦ ; slice number = 48, and 140 time points.

Resting-state fMRI data were preprocessed using Statistical Parametric Mapping (SPM122 ), Data Processing and Analysis for (Resting-State) Brain Imaging (DPABI; Yan et al., 2016) and the rs-fMRI Data Analysis Toolkit (REST 1.8; Song et al., 2011) packages. The following steps were performed: removing the first 10 time points; slice-timing correction; image realignment; normalization to the Montreal Neurological Institute (MNI) space (resampled into 3 mm × 3 mm × 3 mm voxels). The linear trends of time courses were removed, and the effect of nuisance covariates was removed by signal regression using the global signal, the motion parameters, the cerebrospinal fluid (CSF) and white matter (WM) signals. Temporal filtering (0.01 Hz < f < 0.2 Hz) was applied. Finally, each voxel time series was standardized to a mean of zero and standard deviation of unity.

The analysis of the GMV was performed according to the voxel-based morphometry (VBM) protocol using DPABI. The VBM procedure involves the segmentation of the original anatomic MRI images in gray matter (GM), WM, and CSF tissues, followed by GM image normalization to templates in stereotactic space to acquire optimized normalization parameters, which were applied to the raw images. Finally, GM images were smoothed using a 6-mm full-width at half-maximum (FWHM) Gaussian kernel.

#### MSE Theory

Multiscale entropy is based on the theory of SE over a range of scales and consists of two steps (Costa et al., 2002).

(1) The coarse-graining procedure of time series represents the system dynamics on different scale factors. Given time series {xi, i = 1, 2, . . . , N}, for the time scale l, the coarse-grained time series {y l } is calculated as follows:

$$\chi\_j^l = \frac{1}{l} \sum\_{i=(j-1)l+1}^{jl} \varkappa\_i, 1 \le j \le N/l \tag{1}$$

The length of new time series is N/l. For scale 1, the new time series corresponds to the original time series.

(2) The SE for each coarse-grained time series is calculated.

Sample entropy (Richman and Moorman, 2000) is calculated as:

$$SE(m, r, N) = -\ln \frac{P\_{m+1}(r)}{P\_m(r)}\tag{2}$$

where m is the sequence length of time points to be compared, r is the radius of similarity, N is the length of the time series, and P is the probability that points falling within r.

Multiscale entropy consists of a set of SE values under multiple time scales, which reflects the complexity of time series on multiple scales. MSE can be used to compare the complexity of different time series, based on the specific trend of SE changes with scales (e.g., complex time series show constant entropy over

2http://www.fil.ion.ucl.ac.uk/spm

various time scales, while random noise shows a marked decrease in entropy at longer time scales; Wang et al., 2018).

#### MSE Calculation

We used the Complexity Toolbox3 [Laboratory of Functional MRI Technology (LOFT), Department of Neurology, University of Southern California] to calculate MSE of rs-fMRI data.

Three parameter values were set for the calculation of MSE, including pattern length m, distance threshold r, and time scale l. The point to be made is that the r value is generally correlated with the standard deviation of the original time series (Lu et al., 2015). Various theoretical and clinical applications have indicated that, m = 1 or 2 and r = 0.1–0.35 of the standard deviation of the original sequence, provides reasonable statistical validity for calculating SE (Richman and Moorman, 2000). Because no rigorous standard exists for choosing the parameters to calculate SE, prior studies on SE analysis of biomedical signals have shown inconsistent selection of parameters.  In addition, different parameters were also used in the studies of EEG, including m = 2 and r = 0.15 (Catarino et al., 2011), m = 2 and r = 0.25 (Xiang et al., 2015), m = 1 and r = 0.25 (Escudero et al., 2006). In this study, MSE was calculated for each BOLD time series based on different parameter pairs: (m = 2, r = 0.15), (m = 2, r = 0.25), (m = 2, r = 0.30), (m = 2, r = 0.35), (m = 1, r = 0.25), and (m = 1, r = 0.35) across the range of scales from 1 to 6.

#### Statistical Analyses

For every time scale, one-way analysis of variance (ANOVA) was used to assess differences in MSE maps of BOLD signals among four groups (NC, EMCI, LMCI, and AD) using REST 1.8. For multiple comparison corrections, a stringent statistical significance level was employed by setting a voxelwise threshold of p < 0.001 and a cluster threshold of p < 0.05 with a Gaussian random field (GRF) correction among four groups after adjusting for age, sex, and education differences.

Then, fivefold cross-validation was used for regions of interest (ROIs) analyses. We divided the data into five independent subsets. For each fold, we used one subset for selective analysis after using other four subsets for selection (ANOVA). According to the peak MNI coordinates (X Y Z), we extracted the average MSE and GMV by using DPABI toolbox to define ROIs and the radius of the spheres at all scales (8 mm). For each ROI, differences on MSE values among four groups at all scales were compared using ANOVA using Statistical Package for Social Sciences (SPSS 20.0) software. Bonferroni's post hoc pairwise test on ANOVA was performed.

Spearman's correlation was used to assess the relationship between MSE and MMSE, FAQ, CDR, and GMV for four groups using SPSS 20.0 software.

3http://loft-lab.org/index-5.html

## RESULTS

### Demographic and Clinical Characteristics

fnins-12-00677 September 27, 2018 Time: 19:7 # 4

**Table 1** summarizes the demographic and clinical characteristics of the participants. The p-values were obtained by one-way ANOVA. The results indicated no difference in age, sex, and education across four groups. Significant differences (p < 0.001) among four groups were found on the MMSE, FAQ, and CDR scores.

### Parameter Selection for MSE Calculation

The comparison was made by calculating MSE using six different parameter combinations (m, r). All subjects' MSE maps were calculated across time scales from 1 to 6. We performed the one-way ANOVA on MSE maps of four groups on every time scale. Based on the final results, the findings using m = 2 and r = 0.35 as the optimal parameter were mainly reported in this study. Previous study has shown that the accuracy of the calculation results is least dependent on the sequence length N when m = 2 (Smith et al., 2014). Other results are presented in **Supplementary Data Sheet 1**. For m = 2 and r = 0.15, four clusters were significantly different among the four groups across multiple time scales on scale 2, scale 3, scale 4, and scale 6 (**Supplementary Table S1** and **Supplementary Figure S1**). Five clusters were found on scale 2, scale 4, and scale 6 when m = 2 and r = 0.25 (**Supplementary Table S2** and **Supplementary Figure S2**). For m = 2, r = 0.30 and m = 2, r = 0.35, similar results were found and nine clusters showed significant differences on scale 2, scale 4, scale 5, and scale 6 (**Supplementary Table S3** and **Supplementary Figure S3**). For m = 1, r = 0.25 and m = 1, r = 0.35, however, only one consistent cluster was found (left middle occipital gyrus) on scale 1 and no difference was found on the rest scales (**Supplementary Tables S4, S5** and **Supplementary Figures S4, S5**).

### Significant Differences on MSE Among the Four Groups

Using m = 2 and r = 0.35, the result is presented in **Figure 1**. The detailed information is summarized in **Table 2**. Significant differences (p < 0.001, GRF correction) were found on the MSE maps among the four groups on scale 2, scale 4, scale 5, and scale 6. We found no significant difference on scale 1 and scale 3. On scale 2, one cluster was found: right thalamus (THA.R). On scale 4, one cluster was found: left superior frontal gyrus (SFGdor.L). We found two clusters on scale 5: right lingual gyrus (LING.R) and right insula (INS.R). For the scale 6, five clusters were found: right superior temporal gyrus (STG.R), left middle temporal gyrus (MTG.L), right olfactory cortex (OLF.R), left inferior occipital gyrus (IOG.L), and right supramarginal gyrus (SMG.R).

We also extracted the mean MSE of whole brain (WB), GM, WM, and CSF using the corresponding masks on all time scales. Then, one-way ANOVA was performed to examine the differences among the four groups. The result is presented in **Supplementary Table S6**. Only GM showed a trend of entropy difference (F = 2.283, p = 0.083) among four groups on scale 6. **Figure 2** shows the mean entropy curve of GM across the scale of 1–6 for four groups as well as the differences between each pair of the four groups on scale 6 (p < 0.05, two-sample t-test, uncorrected).

### Time Scales Analysis on MSE From Scale 1 to Scale 6

We extracted the average MSE of 9 ROIs over multiple time scales. **Figure 3** displays the MSE curve across scale 1 to 6 among four groups (NC, EMCI, LMCI, and AD) for nine ROIs. Each group exhibited a drop in SE values with increasing scale. SE values on scale 1 showed no difference among the four groups for nine ROIs. For scale 2, there were two ROIs (THA.R and OLF.R) showing significant differences among four groups. SFGdor.L, INS.R, and OLF.R showed significant differences on scale 3. For scale 4, there were six ROIs (SFGdor.L, LING.R, INS.R, MTG.L, OLF.R, and SMG.R) showing significant differences among four groups. There were four ROIs showing significant differences (LING.R, INS.R, and IOG.L) on scale 5. On scale 6, there were six ROIs showing significant differences (THA.R, STG.R, MTG.L, OLF.R, IOG.L, and SMG.R) among four groups. Specifically, OLF.R showed significant differences on four time scales (scale 2, scale 3, scale 4, and scale 6).

### Comparison of MSE Among the Four Groups

Multiscale entropy values of nine ROIs at all scales were compared among the four groups (NC, EMCI, LMCI, and AD) using ANOVA, and for MSE of ROIs with significant differences among the four groups, Bonferroni's post hoc pairwise test on ANOVA was performed. **Figure 4** displays the comparison of MSE of nine ROIs at different scales between any two groups. The results showed that, compared with NC subjects, patient groups demonstrated reduced complexity. Specifically, the AD group showed lower complexity than the NC group for all ROIs.

Compared with the NC group, the EMCI subjects had significantly reduced MSE of BOLD signals in INS.R. The LMCI subjects showed significantly decreased MSE in eight of the nine ROIs except THA.R. Compared with the EMCI group, the LMCI group showed decreased MSE in SFGdor.L while the AD group showed decreased MSE in three ROIs (THA.R, MTG.L, and OLF.R). In addition, THA.R had lower complexity in the AD group than that in the LMCI group.

### Relationships Between MSE and Clinical Measurements

We performed Spearman's correlations between MSE values and the clinical measurements (MMSE, FAQ, and CDR) in patient groups (MCI and AD). After corrections for multiple comparisons, significant correlations were found.

**Figure 5** shows the scatter plots between MSE values of BOLD signals and clinical measurement scores (MMSE, FAQ, and CDR) in patient groups in the significantly correlated brain regions. On scale 2, MMSE was positively correlated with the complexity of BOLD signals in THA.R (r = 0.354, p = 0.006). SFGdor.L

TABLE 2 | Characteristics of the brain regions those were significantly different among the four groups across multiple time scales.

The location coordinates are those of the peak significance in each region (p < 0.001, GRF corrected).

exhibited the significant positive correlation (r = 0.293, p = 0.030) between the MSE values and MMSE scores on scale 4. The four ROIs (INS.R, STG.R, IOG.L, and SMG.R) exhibited significant positive correlations (r > 0.283, p < 0.048) between MSE and MMSE scores on scale 6. Some trend correlations were also found (p < 0.05, uncorrected) and the results are shown in **Supplementary Table S7**.

As shown in **Figure 5**, the MSE values of THA.R exhibited significant negative correlations (r = −0.344, p = 0.006) with the FAQ scores on scale 2. OLF.R exhibited significant correlations (r = −0.291, p = 0.042) between FAQ scores and MSE values of BOLD signals in patient groups on scale 6. SFGdor.L, LING.R, INS.R, and IOG.L exhibited trend correlations (p < 0.05, uncorrected) between FAQ scores and MSE values on multiple time scales (**Supplementary Table S8**).

After corrections for multiple comparisons, MSE values of THA.R exhibited the significant negative correlation with the CDR scores on scale 2 (r = −0.303, p = 0.024) and scale 3 (r = −0.286, p = 0.042). LING.R showed the significant negative correlation between CDR scores and MSE values on scale 5 (r = −0.331, p = 0.012) and MSE values of SMG.R were negatively correlated with CDR scores on scale 6 (r = −0.312, p = 0.018; **Figure 5**). In addition, SFGdor.L, INS.R, STG.R, OLF.R, and IOG.L exhibited trend correlations (p < 0.05, uncorrected) between CDR scores and MSE values on multiple time scales. **Supplementary Table S9** summarizes the correlation results between CDR scores and MSE values of BOLD signals in patient groups.

We also performed Spearman's correlations between MSE values and the clinical measurements (MMSE, FAQ, and CDR)

FIGURE 2 | (A) MSE curve across scale factor 1–6 in gray matter (GM) for four groups. Each point represents group average SE. (B) Mean SE values of GM in the NC, EMCI, LMCI, and AD subjects on scale 6. Significant differences between each pair of the four groups (p < 0.05, two-sample t-test, uncorrected) are indicated. The error bar represents the standard error of MSE within the group. ∗ indicates p < 0.05.

for each group (NC, EMCI, LMCI, and AD). After corrections for multiple comparisons, for the NC group, SMG.R exhibited the significant negative correlation (r = −0.516, p = 0.048) between MSE values and FAQ scores on scale 5. No correlation was found between MSE values and MMSE, FAQ, and CDR scores in the EMCI group on all scales. For the LMCI group, MSE values of STG.R were positively correlated (r = 0.512, p = 0.030) with MMSE scores and MSE values of MTG.R and IOG.L were

negatively correlated (r < −0.485, p < 0.048) with FAQ scores. STG.R exhibited the significant negative correlation(r = -0.488, p = 0.048 between MSE values and CDR scores in the AD group. Some trend correlations were also found (p < 0.05, uncorrected) and the results are shown in **Supplementary Tables S11**–**S14**.

#### Relationships Between MSE and GMV

We extracted the average GMV values of nine ROIs for four groups. Then, we explored the relationships between the MSE and the GMV in patient groups. After corrections for multiple comparisons, no significant correlation was found between the complexity of BOLD signals and GMV values. But the LING.R exhibited trend positive correlations (r > 0.209, p < 0.043, uncorrected) between the MSE and the GMV in patient groups on four time scales (scale 3, scale 4, scale 5, and scale 6). STG.R exhibited a positive correlation (r = 0.203, p = 0.050, uncorrected) between the MSE and the GMV on scale 6 and MTG.L showed a positive correlation (r = 0.235, p = 0.023, uncorrected) on scale 5. The results are presented in **Supplementary Table S10**.

Correlation analyses for each group (NC, EMCI, LMCI, and AD) were also performed between the MSE and the GMV on all time scales. After corrections for multiple comparisons, no significant correlation was found between the complexity of BOLD signals and GMV values in the NC, EMCI, and LMCI groups. Only SFGdor.L showed the significant positive correlations in the AD group on scale 6. Some brain regions exhibited trend positive correlations (p < 0.05, uncorrected) between MSE values and GMV values in each group and the results are presented in **Supplementary Tables S11**–**S14**.

### DISCUSSION

In this study, we employed MSE analysis to assess the complexity of BOLD activity in AD and MCI patients from scale 1 to 6. We discovered that the spontaneous BOLD signals of nine clusters had significant differences among four groups on four time scales. The significant MSE differences were mainly detected in the occipital, frontal, temporal, limbic, and parietal lobes, which were significantly correlated with clinical measurements in patient groups from scale 2 to 6. These results suggest that the complexity analyses using MSE of BOLD signals can provide information on the temporal dynamics of neural signals across multiple scales that are relevant to the cognitive impairments in MCI and AD.

### The MSE Differences Among Four Groups on Multiple Time Scales

This study found that MSE of BOLD activity exhibited significant contrasts among four groups on four time scales (scale 2, scale 4, scale 5, and scale 6; **Figure 1**), mainly distributed in the occipital lobe (IOG.L and LING.R), frontal lobe (SFGdor.L and OLF.R), parietal lobe (SMG.R), temporal lobe (STG.R and MTG.L), limbic lobe (INS.R), and the subcortical region (THA.R). In the MSE analysis for nine ROIs over all time scales, we found that three ROIs (SFGdor.L, INS.R, and OLF.R) had significant differences on scale 3 (**Figure 3**). This means that more useful information was found on multiple time scales. This is consistent with previous reports using MSE analysis on rs-fMRI and EEG signals that detected differences in entropy on multiple time scales (Mizuno et al., 2010; Liu et al., 2013; Yang et al., 2013; Mcbride et al., 2014; Smith et al., 2014; Michalopoulos and Bourbakis, 2017). Particularly, as **Figure 1** demonstrates, five clusters showed significant differences on scale 6. More significant differences were found among the four groups with increasing scales. As can be seen from **Figure 3**, six ROIs showed prominent differences among the four groups on scale 4 and 6. Based on the mechanism of MSE analysis, at the shortest scale, the entropy is dominated by the high

frequency fluctuations from random noise (Wang et al., 2018). By filtering out these random fluctuations, the contrast in entropy becomes larger at longer time scales (Smith et al., 2014).

In this study, each of the nine ROIs was observed on a single scale. In the process of calculating MSE, the key step is to coarse-grain the time series to reflect the system dynamics on different time scales, which means that, MSE mainly calculates the complexity of high frequencies at low scales, while the complexity of low frequencies is calculated at large scales. Our results showed that different brain regions displayed differences at different frequencies. Consisted with our result, Wang et al. (2018) investigated the neurophysiological underpinnings of complexity (MSE) of fMRI signals and their relations to FC and the results showed that the associations between MSE and FC were dependent on the temporal scales or frequencies It has been proposed that each frequency band is generated by different mechanisms and relates to different physiological functions, higher frequency oscillations are confined to a small neuronal space, whereas lower frequencies may reflect longrange interactions (Buzsáki and Draguhn, 2004; Zuo et al., 2010). More recently, studies on rs-fMRI have hypothesized that frequency-dependent effects in different brain regions which reflect different synaptic and functional characteristics that are affected by the progression of cognitive impairment (He et al., 2010; Yu et al., 2013; Wang et al., 2016; Zhou et al., 2016). Hence, we propose that the observed complexity changes on different time scales might represent different region or networkdependent neuropathophysiological mechanisms in MCI and AD.

We also analyzed the complexity of WB, GM, WM, and CSF on all time scales. Only GM showed a trend of entropy difference (F = 2.283, p = 0.083) among four groups on scale 6. Many studies on the complexity analysis of rs-fMRI data found global complexity differences in aging (Yang et al., 2013; Sokunbi et al., 2015) and AD (Liu et al., 2013; Wang et al., 2017). Smith et al. (2014) found greater age-related decline in average GM complexity of rs-fMRI at longer time scales, and Liu et al. (2013) found mean complexity of rs-fMRI in GM and WM decreased with normal aging. Thus, the complexity of global brain activity may decrease with age. For AD-related cognitive decline, Liu et al. (2013) found that mean ApEn of GM showed a significant positive correlation with MMSE scores in the cohort of familial AD subjects, and Wang et al. (2017) found significant differences (p < 0.05) in permutation entropy (PE) of GM and WM across the four groups of ADNI data. Possibly due to differences in data samples and complexity analysis methods, the MSE analysis used in this study was only able to reveal a trend of entropy difference (p = 0.083) among four groups in GM as well as decreased complexity in the AD and LMCI groups compared to that of the NC group at the statistical threshold of uncorrected p < 0.05. In contrast, we found highly significant MSE differences (p < 0.001, GRF corrected) in several brain regions on multiple time scales. This is not surprising as the pathological process of AD first affects the network of temporal, frontal, and parietal regions before progressing to the whole GM and brain level. Different complexity analyses may have different sensitivities in detecting global and regional changes of neural complexity with AD progression. This question awaits further investigation in future studies.

### Decreased Complexity and Cognitive Decline in Patient Groups

Using the post hoc pairwise test on ANOVA, reduced complexity in the AD group was detected in all ROIs compared with the NC group (p < 0.05, Bonferroni corrected). In addition, MSE also showed strong sensitivity in differentiating NC from EMCI (one ROI), NC from LMCI (eight ROIs), EMCI from LMCI (one ROI), EMCI from AD (three ROIs), and LMCI from AD (one ROI). As can be seen from **Figure 4**, the complexity of BOLD signals in most ROIs showed a gradually decline from NC to EMCI to LMCI and to AD. Previous complexity studies of fMRI signals also consistently reported reduced complexity in AD patients compared to matched control subjects (Liu et al., 2013; Wang et al., 2017). Liu et al. (2013) reported decreased complexities in STG, MTG, and SMG in familial AD. Some of brain regions, such as SFGdor, MTG, and IOG, were also reported in our previous study using PE method to analyze the complexity of the same ADNI dataset (Wang et al., 2017).

We performed correlation analyses between complexity of BOLD signals in these brain regions with significant MSE differences and cognitive function scores (MMSE, FAQ, and CDR). These three clinical measurements provide quantitative assessments of cognitive function and are widely used (Ciesielska et al., 2016; Kaur et al., 2016; Kim et al., 2017). Higher scores of MMSE indicate higher aptitude of cognition; low functional performance is related to higher FAQ scores and the presence of dementia is indicated by higher CDR scores. Our correlation results showed that the average MSE of some brain regions was significantly positively correlated with the MMSE scores and significantly negatively correlated with FAQ scores and CDR scores in patient groups (p < 0.05, corrected). This means that lower MMSE and higher FAQ and CDR scores were observed in MCI and AD patients who exhibited lower MSE in some brain regions. Particularly, THA.R exhibited significant correlations between MSE values and three clinical measurement scores (MMSE, FAQ, and CDR) on scale 2. The MSE values of SMG.R showed significant correlations with the MMSE and CDR scores on scale 6. Previous fMRI studies suggested that THA and SMG are closely related to cognitive dysfunction in healthy aging and AD (Mevel et al., 2011; Yang et al., 2013; Xiaoying et al., 2014; Raczek et al., 2017). Studies found that activity in THA is associated with spatial working memory and memory processing (Jankowski et al., 2013; Saalmann and Kastner, 2015; Štillová et al., 2015; Hallock et al., 2016), and SMG is mainly involved in language perception, phonological processing and verbal working memory and processing (Hartwigsen et al., 2010; Kheradmand et al., 2013; Deschamps et al., 2014). In addition, as can be seen from **Figure 5**, significant correlations between MSE and cognitive measurements were dependent on the temporal scales. For example, we observed THA.R showed associations between complexity and MMSE at high temporal frequencies, and SMG.R exhibited significant correlations at low

temporal frequencies. The results showed that different brain regions displayed correlations at different frequencies and once again corroborated the MSE theory that high and low temporal frequencies may represent region or network-dependent different neuropathophysiological mechanisms (Buzsáki and Draguhn, 2004; Zuo et al., 2010).

### Potential Physiological Underpinnings of Altered Complexity in Patient Groups

It has been suggested that physiological diseases are associated with a loss of complexity in healthy systems (Lipsitz, 2004; Pincus, 2010). AD is a neurodegenerative disorder characterized by dementia and cognitive decline (Querfurth and Laferla, 2010). The brain regions that we found to have reduced complexity play important roles for cognitive functions. For example, the lingual gyrus is believed to play a role in the analysis of logical conditions and encoding visual memories. The superior temporal gyrus is involved in social cognition processes and middle temporal gyrus is mainly involved in the recognition of known faces and episodic memory (Bigler et al., 2007; Acheson and Hagoort, 2013). Some fMRI experiments have found proof that the superior frontal gyrus is involved in self-awareness, sensory system, and social cognitive processes (Goldberg et al., 2006). The altered complexity of these brain regions in patient groups may be associated with deterioration of brain function in these important networks.

Further, AD is characterized by the presence of neuritic plaques and neurofibrillary tangles, accompanied by widespread cortical neuronal loss, and loss of connections between brain systems (Sankari, 2010). This may degrade cortical and subcortical connections, leading to cognitive and behavioral disturbances. Many studies have reported that disrupted FC in the AD group in THA, SFGdor, INS, STG, MTG, IOG, and SMG (Zhang et al., 2009; Wang et al., 2010; Dennis and Thompson, 2014). Thus, this degeneration of both local and long-range connections disrupts the functional coherence of brain activation, decreasing the complexity of spontaneous brain activity.

In addition, we examined the relationships between MSE and GMV in patient groups. After corrections for multiple comparisons, no significant correlation was found between the complexity of BOLD signals and GMV. But LING.R, STG.R, and MTG.L exhibited trend positive correlations (p < 0.05, uncorrected) between the MSE and the GMV in patient groups. Many studies have also reported GM atrophy in these brain regions in MCI and AD (Busatto et al., 2003; Karas et al., 2004; Guo et al., 2010; Möller et al., 2013). In our previous study, we also found that the complexity of these brain regions was related to GMV and was associated with glucose metabolism (Wang et al., 2017). More pathologies of AD may lead to lower complexity of brain regions still requires further study.

### Comparison of SE, PE, MSE, and Multiscale PE

Sample entropy solved the problem of vector self-matching in the ApEn defined by the Heaviside function and has been widely used (Pincus, 1991; Richman and Moorman, 2000). PE is different from SE, as PE calculates the probability of a symbolic sequence of points in the phase space and the entropy value in the form of Shannon information entropy (Bandt and Pompe, 2002). Many researchers prefer to use SE and PE to study the complexity of time series and obtain findings on a single scale (Sokunbi et al., 2013; Berger et al., 2017; Wang et al., 2017; Aktaruzzaman, 2018). Compared with PE and SE, MSE and multiscale PE (MPE) investigate the dynamic complexity of time series data across multiple temporal scales, not only at the original time scale of 1 (Costa et al., 2002; Aziz and Arif, 2005; Ouyang et al., 2013).

In this study, we found significant complexity differences among four groups on multiple temporal scales, especially on longer time scales, due to MSE's capability to average out short time scale fluctuations (Smith et al., 2014; Yan et al., 2017). Thus, compared with SE, researchers performed MSE for complexity analysis obtained richer and more comprehensive information in aging and neurodegenerative diseases (Humeauheurtier, 2016; Shang, 2017). Our previous research investigated the abnormal complexity of BOLD signals in MCI and AD patients using PE analysis (Wang et al., 2017). Then, we also applied MPE to the same dataset, but no significant difference was found on longer scales (p < 0.005, GRF correction). Some studies demonstrated that PE had better anti-noise performance and thus, compared with SE, we got supplementary information in detecting differences among four groups on scale 1 (Bandt and Pompe, 2002; Nicolaou and Georgiou, 2012; Wang et al., 2017). The coarse-grained procedures in MPE with large scale factors may result in short data length, while PE requires more time points to contain more states of the reconstructed sequence (Bandt and Pompe, 2002). This may be the reason that we did not detect the significant PE difference at longer time scales. As a consequence, for our dataset, MPE had better performance at short time scales, while MSE could provide more information on multiple time scales. In the future, we will perform and compare SE and PE analysis across multiple time scales on more rs-fMRI datasets to further our understanding on this issue.

#### Limitation

A limitation of this study is the short BOLD time series used for MSE analysis which may lead to potentially erratic entropy estimation (Costa et al., 2002; Yang et al., 2013). In this study, we performed the parameter selection for MSE calculation by using 6 different parameter pairs based on previous studies, not all of the parameter pairs. The results of m = 2 and r = 0.35 were mainly reported in this study. The selection of parameters may be related to particular datasets, and different datasets may have different optimal parameters.

## CONCLUSION

Multiscale entropy is a powerful tool to quantify the nonlinear information of a time series over multiple time scales through the SE algorithm. This study applied MSE analysis to investigate the abnormal complexity of BOLD signals across multiple time scales in MCI and AD patients. Enhanced MSE differences were detected among four groups which were significantly correlated with clinical assessments in patient groups at multiple temporal scales. The MCI and AD patients demonstrated lower complexity than normal controls and AD patients showed lower complexity than MCI. These findings indicate that MSE of spontaneous BOLD signals may provide an imaging marker of cognitive impairment in MCI and AD.

#### ALZHEIMER'S DISEASE NEUROIMAGING INITIATIVE

fnins-12-00677 September 27, 2018 Time: 19:7 # 11

The data used in preparation of this article were obtained from the Alzheimer's disease neuroimaging initiative (ADNI) database adni.loni.usc.edu. As such, the investigators within the ADNI contributed to the design and implementation of the ADNI and/or provided data but did not participate in analysis or writing of this report. A complete listing of the ADNI investigators can be found at http://adni.loni.usc.edu/wpcontent/uploads/how\_ to\_apply/ADNI\_Acknowledgement\_List.pdf.

#### AUTHOR CONTRIBUTIONS

YN performed the experiment and completed the manuscript. BW, MZ, JX, HS, RC, and XC provided suggestions for this study. JX provided the guidance throughout the study.

#### FUNDING

This study was supported by the National Natural Science Foundation of China (61503272, 61305142, 61741212, and 61373101), the Natural Science Foundation of Shanxi Province<PHONE_NUMBER> and 201601D202042), a project funded by the China Postdoctoral Science Foundation (2016M601287), and the Shanxi Provincial Foundation for Returned Scholars, China (2016-037).

### REFERENCES

Data collection and sharing for this project was funded by the Alzheimer's disease neuroimaging initiative (ADNI) (National Institutes of Health Grant U01 AG024904) and DOD ADNI (Department of Defense award number W81XWH-12-2-0012). ADNI was funded by the National Institute on Aging and the National Institute of Biomedical Imaging and Bioengineering and through generous contributions from the following: AbbVie, Alzheimer's Association; Alzheimer's Drug Discovery Foundation; Araclon Biotech; BioClinica, Inc.;Biogen; Bristol-Myers Squibb Company; CereSpir, Inc.; Cogstate; Eisai Inc.; Elan Pharmaceuticals, Inc.; Eli Lilly and Company; EuroImmun; F. Hoffmann-La Roche Ltd and its affiliated company Genentech, Inc.; Fujirebio; GE Healthcare; IXICO Ltd.; Janssen Alzheimer Immunotherapy Research and Development, LLC.; Johnson & Johnson Pharmaceutical Research and Development LLC.; Lumosity; Lundbeck; Merck & Co., Inc.; Meso Scale Diagnostics, LLC.; NeuroRx Research; Neurotrack Technologies; Novartis Pharmaceuticals Corporation; Pfizer Inc.; Piramal Imaging; Servier; Takeda Pharmaceutical Company; and Transition Therapeutics. The Canadian Institutes of Health Research provides funds to support ADNI clinical sites in Canada. Private sector contributions are facilitated by the Foundation for the National Institutes of Health (www.fnih.org). The grantee organization is the Northern California Institute for Research and Education, and the study is coordinated by the Alzheimer's Therapeutic Research Institute at the University of Southern California. ADNI data are disseminated by the Laboratory for Neuro Imaging at the University of Southern California.

### SUPPLEMENTARY MATERIAL

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

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Niu, Wang, Zhou, Xue, Shapour, Cao, Cui, Wu and Xiang. 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.

# Topological Pattern Recognition of Severe Alzheimer's Disease via Regularized Supervised Learning of EEG Complexity

Miaolin Fan1 , Albert C. Yang2,3, Jong-Ling Fuh4,5 and Chun-An Chou1 \*

1 Department of Mechanical and Industrial Engineering, Northeastern University, Boston, MA, United States, 2 Division of Interdisciplinary Medicine and Biotechnology, Beth Israel Deaconess Medical Center, Harvard Medical School, Boston, MA, United States, 3 Institute of Brain Sciences, National Yang-Ming University, Taipei, Taiwan, 4 Neurological Institute, Taipei Veterans General Hospital, Taipei, Taiwan, 5 School of Medicine, National Yang-Ming University, Taipei, Taiwan

Alzheimer's disease (AD) is a progressive brain disorder with gradual memory loss that

#### Edited by:

Laura Marzetti, Università degli Studi G. d'Annunzio Chieti e Pescara, Italy

#### Reviewed by:

Elzbieta Olejarczyk, Institute of Biocybernetics and Biomedical Engineering (PAN), Poland Filippo Zappasodi, Università degli Studi G. d'Annunzio Chieti e Pescara, Italy

> \*Correspondence: Chun-An Chou<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Brain Imaging Methods, a section of the journal Frontiers in Neuroscience

Received: 28 April 2018 Accepted: 12 September 2018 Published: 04 October 2018

#### Citation:

Fan M, Yang AC, Fuh J-L and Chou C-A (2018) Topological Pattern Recognition of Severe Alzheimer's Disease via Regularized Supervised Learning of EEG Complexity. Front. Neurosci. 12:685. doi: 10.3389/fnins.2018.00685 correlates to cognitive deficits in the elderly population. Recent studies have shown the potentials of machine learning algorithms to identify biomarkers and functional brain activity patterns across various AD stages using electroencephalography (EEG). In this study, we aim to discover the altered spatio-temporal patterns of EEG complexity associated with AD pathology in different severity levels. We employed the multiscale entropy (MSE), a complexity measure of time series signals, as the biomarkers to characterize the nonlinear complexity at multiple temporal scales. Two regularized logistic regression methods were applied to extracted MSE features to capture the topographic pattern of MSEs of AD cohorts compared to healthy baseline. Furthermore, canonical correlation analysis was performed to evaluate the multivariate correlation between EEG complexity and cognitive dysfunction measured by the Neuropsychiatric Inventory scores. 123 participants were recruited and each participant was examined in three sessions (length = 10 seconds) to collect resting-state EEG signals. MSE features were extracted across 20 time scale factors with pre-determined parameters (m = 2, r = 0.15). The results showed that comparing to logistic regression model, the regularized learning methods performed better for discriminating severe AD cohort from normal control, very mild and mild cohorts (test accuracy ∼ 80%), as well as for selecting significant biomarkers arcoss the brain regions. It was found that temporal and occipitoparietal brain regions were more discriminative in regard to classifying severe AD cohort vs. normal controls, but more diverse and distributed patterns of EEG complexity in the brain were exhibited across individuals in early stages of AD.

Keywords: Alzheimer's disease, EEG, complexity analysis, pattern recognition, LASSO

## 1. INTRODUCTION

Alzheimer's disease (AD) is a neurodegenerative disorder characterized by progressive loss of memory and cognitive dysfunctions. Despite of many efforts, the pathological mechanism of AD progression still remains unsettled. In recent decades, the emerging field of interdisciplinary studies between computational cognitive and data sciences has enabled data-driven knowledge discovery systems for investigating multivariate patterns based on largescale, complex brain data. More specifically, advances of machine learning techniques have contributed to the clinical science by not only improving the automated diagnostic/predictive tools, but also enhancing the understanding of pathological mechanism underlying AD progression. In the past years, there were studies to demonstrate the capability of machine learning algorithms in addressing the sophisticated patterns using various brain data, e.g., electroencephalography (EEG) and magnetic resonance imaging (MRI). Trambaiolli et al. (2011) identified the bipolar peaks of EEG signals as biomarkers for differentiating AD, mild cognitive impairments (MCI) and early dementia patients. Casanova et al. (2011) found that most informative voxels in structural MRI data locate in the gray and white matter tissues, which can discriminate patients from cognitive normal subjects accurately using large-scale regularization. Other studies encouraged the utilization of an integrative EEG biomarkers derived from various sources in order to provide predictive models with diverse and comprehensive information (Poil et al., 2013; Triggiani et al., 2017).

Among modern neuroimaging modalities, EEG as a noninvasive, inexpensive technique has drawn extensive attentions for investigating nonlinear dynamics of neuronal brain functions. It was reported that AD progression can be characterized by the reduced complexity in EEG signals, which is hypothesized to be related to the loss of neurons and possible connectivity caused by pathological aging process. A recent and comprehensive review is refered to Dauwels et al. (2010). In this study, we used Multiscale Entropy (MSE) for estimating the nonlinear complexity of EEG signals across multiple temporal scales (Costa et al., 2002). Previous studies investigated MSE as a measure of complexity for understanding AD pathology using univariate (Escudero et al., 2006; Park et al., 2007) and multivariate EEG dynamics (Labate et al., 2013). It was reported that the decreased complexity in short-time scale and increased complexity in long-time scale distinguish AD patients from normal controls (Mizuno et al., 2010; Yang et al., 2013). A recent study (Azami et al., 2017) also indicated the potentials of the second-order MSE features for characterizing EEG changes with AD progression. Moreover, correlation was found between MSE features from various brain regions and multiple neuropsychiatric symptoms, particularly in temporal and occipitoparietal electrodes (Yang et al., 2013). Since the previous study only assessed the bivariate correlation, we extend to investigate the relationships in a multivariate feature space by applying canonical correlation analysis (CCA) (Hotelling, 1936). CCA is a multivariate technique that is capable to capture multiple causes and effects to further investigate the relationship between MSE and neuropsychiatric symptoms.

One of the most challenging tasks for understanding AD pathology is to characterize the biomarkers and associated patterns that differentiate different AD severity levels. Considering the pathological aging of the brain is a highly heterogeneous process, the generalizability in many existing research studies is limited by the small sample size, large individual variability, and high-dimensional data structure. While most state-of-the-art machine learning algorithms suffered from over-fitting data and produced poor generalized prediction results, regularized learning methods attempt to address this over-fitting issue by adding a regularization term (called L1-norm or L2-norm) to the cost function. Least Absolute Shrinkage and Selection Operator (LASSO) (Tibshirani, 1996) is a classic method that builds a regression model of correlating input variables (MSE features in this study) to the prediction outcome (severe AD or not) while posing a penalty on the number of non-zero coefficients of input variables (L1-norm feature selection). Later, Elastic net (ENet) was proposed by combining L1-norm and L2-norm for the purpose of addressing several drawbacks of LASSO, including the group effect among input variables in addition to feature selection. The flexibility and variability of regularization methods allow one to develop variants for specific purposes (Tibshirani et al., 2005; Bach, 2008). Specifically, the interpretability/stability of feature selection is desired for providing scientific insights, since consistent feature selection across different samples and individuals is more likely to suggest a meaningful pattern (Fan and Chou, 2016). Stability selection (Meinshausen and Bühlmann, 2010) is thus proposed based on the combination of feature selection method and repeated subsampling. For the cost of computational resources, the stability selection aims to provide a statistical control on the error rate of feature selection in a sparse dataset.

Based on the general concept of stability selection approach, the present study intends to provide a stability-based feature selection and identify important EEG biomarkers using the frequency of selection across multiple replicates in cross validation. The objective of our study is two fold. On one hand, we are interested in characterizing the functional brain activities with varying temporal scales that best discriminate severity levels of AD groups and normal controls based on EEG complexity. On the other hand, we aim to profile the topolographic map of EEG biomarkers for various AD severity and investigate the multivariate correlation patterns to cognitive dysfunctions.

### 2. MATERIALS AND METHODS

### 2.1. Participants

One hundred and twenty-three participants were recruited from the Dementia Clinic at the Neurological Institute, Taipei Veterans General Hospital in Taiwan. The diagnosis for AD was based on the criteria of the National Institute of Neurological and Communicative Disorders and the Stroke/Alzheimer's Disease and Related Disorders Association (McKhann et al., 1984). All patients had received neurological examinations, laboratory tests, EEG monitoring, and neuroimaging evaluation during the diagnostic process. Our study was approved by the Institutional Review Board of Taipei Veterans General Hospital to conduct retrospective analysis of the patients' clinical and EEG data. We excluded patients who had other conditions that caused secondary dementia, such as vascular dementia, Parkinson's disease, hypothyroidism, vitamin B12 deficiency, syphilis, and prior history of major psychiatric illness (e.g., major depression, bipolar disorder, or schizophrenia). The participants were categorized into four groups according to their severity of dementia, assessed by the Clinical Dementia Rating (CDR) scale (Morris, 1993). In the following sections, we refer to these groups as HC (healthy control; N = 15), AD1 (very mild, CDR = 0.5; N = 15), AD2 (mild, CDR = 1; N = 69), and AD3 (moderate to severe, CDR = 2; N = 24).

#### 2.2. EEG Data Acquisition and Pre-processing

A routine EEG recordings were performed on all participants (Nicolet EEG, Natus Medical, Incorporated, San Carlos, CA, USA) in the EEG examination room at the Neurological Institute of Taipei Veterans General Hospital. The EEG recording protocol began with a 5-min habituation to the examining environment, followed by three consecutive sessions of 10–20 s with the eyes closed and then open, and a session of photo stimulation, while only the eye closed data was used in the present study. The recordings were performed using the international 10–20 system of 19 electrodes (Fp1, Fp2, F7, F3, Fz, F4, F8, T3, C3, Cz, C4, T4, T5, P3, Pz, P4, T6, O1, and O2) with linked ear reference, 256 Hz sampling rate and filtered at 0.05 Hz high-pass, 70 hZ low-pass and notch filter of 60 Hz, and impedance below 3 k. Vigilance was monitored by the EEG technician, who alerted patients when signs of drowsiness appeared in the tracings. Vertical eyeball movement was detected from electrodes placed above and below the right eye, while the horizontal eyeball movement was detected from electrodes placed at the left outer canthus. EEG signals were preprocessed to remove the linear trend and visually inspected to ensure there were no eye movement artifacts. The EEG signals were exported in European Data Format and were processed using MATLAB 2016b (Mathworks, Inc.).

#### 2.3. Multiscale Entropy Analysis (MSE)

In this study, we employed MSE (Costa et al., 2002) to measure the nonlinear complexity of EEG signal. Let us consider a singlechannel EEG signals with length = N, denoted by {x1, x2, ....xN}. MSE provides an estimate of the sample entropy over multiple time scales in two steps: (1) the construction of coarse-grained time series based on various scale factors, denoted by τ , and (2) the estimation of sample entropy for each time scale. In the first step, the range of τ need to be pre-defined as a set of increasing integers starting from 1. (i.e., [1, 2, ... T]). For each possible value of τ , the corresponding coarse-grained time series yj(τ ) is obtained by applying a non-overlapping sliding window with length = τ and taking the average of all values in each window, represented by the following equation (1 ≤ j ≤ N/τ ):

$$\chi\_j(\tau) = \frac{1}{\tau} \sum\_{i=(j-1)\tau+1}^{j\tau} \chi\_i. \tag{1}$$

If we denote M as the largest integer such that M ≤ N/τ , the coarse-grained time series is then rewritten as {y1(τ ), y2(τ ), ..., yj(τ ), ...yM(τ )}.

In the second step, the sample entropy (Richman and Moorman, 2000) is calculated for each coarse-grained time series as a function of τ . To calculate the sample entropy for a time series with length = M, two parameters need to be determined: the pattern length m and the similarity criterion r. Within the coarse-grained time series {y1(τ ), y2(τ ), ..., yj(τ ), ...yM(τ )}, we denote a vector of pattern length = m as Ym(k) = {yk (τ ), yk+1 (τ ), ..., yk+m−1 (τ )}. Accordingly, the total number of pairs of vectors that satisfy D(Ym(k), Ym(l)) < r(k 6= l) is denoted by Nm. The sample entropy I(τ ) for this time series with parameters τ and r is defined as:

$$I(\mathbf{r}, r) = -\log \frac{N\_{m+1}}{N\_m}.\tag{2}$$

In this study, we use m = 2 and r = 0.15, and the range of scale factors is [1, 20] by following our previous work (Yang et al., 2013). **Figure 1** shows the averaged raw EEG signals, spectral power and MSE scores across all groups; a cross-over is observed in the MSE curves with the increasing scale factors. In short-time scales (≤ 8), lower MSE features are observed from the severe AD group comparing to normal controls, but in long-time scales (> 8) an opposite pattern is observed.

#### 2.4. Hybrid Machine Learning Model for Classification and Biomarker Identification

The objective of applying machine learning model to analyze the MSE features of EEG signals is two-fold: first, we intend to discriminate between control group and AD groups (AD1, AD2, and AD3) by performing a binary classification task in a one-to-one manner (exhaust all the possible combination of pairs). Second, we aim to examine the multivariate correlation patterns between MSE features and dementia symptoms rated by clinicians based on The Neuropsychiatric Inventory (NPI) (Cummings et al., 1994). After extracting MSE features from 19-channel EEG device using 20 scale factors, 380 (= 19 × 20) dimensions were obtained for the feature space. The machine learning model may be over-fitted in training with the relatively less samples on this high dimensional feature space. Therefore, regularization learning methods are employed to perform classification tasks between different AD/HC groups while reducing the dimensionality of trained model. A logistic regression (LR) model is trained and fitted with a penalization on the number of features with non-zero coefficients. As a result, an automatic feature selection is performed by forcing some features to yield zero coefficients. In the following subsections, we present two classic types of regularized LR models. Furthermore, we implement canonical correlation analysis, a unsupervised learning method, for inferring the correlations among two sets of variables.

#### 2.4.1. L1-Norm and L2-Norm Regularized Learning Methods

The original form of LASSO is a linear regression model with a penalty term that controls the number of non-zero coefficients for all variables. In a classification problem, LASSO is reformulated with the cost function of LR, which is rewritten as the following problem (Tibshirani, 1996; Friedman et al., 2001):

$$\max\_{\beta\_0, \beta} \left\{ \sum\_{i=1}^N \left[ \mathbf{y}\_i (\beta\_0 + \boldsymbol{\beta}^T \mathbf{x}\_i) - \log(1 + \boldsymbol{\varepsilon}^{\beta\_0 + \boldsymbol{\beta}^T \mathbf{x}\_i}) \right] - \lambda \sum\_{j=1}^p |\beta\_j| \right\},\tag{3}$$

FIGURE 1 | An illustration of general distribution for each group, including (A) raw EEG signals, (B) spectral powers, and (C) MSE on channel F8 for all groups. These curves show that MSE curves are more distinguishable than EEG signals and spectral powers in overall, and the trend with increasing scale factors in the MSE curve of each AD group is different.

where f(xi) = β0 + β T xi and yi are the prediction and target class for the ith sample respectively. Pp j=1 |βj | is also known as L1-norm penalty that controls the shrinkage with corresponding parameter λ selected via nested cross-validation.

However, LASSO attempts to addresses cluster information of correlated variables, which is referred to as grouping effect. It only selects one and drops the other variables when fitted with a group of related variables (Zou and Hastie, 2005). In this study, this grouping effect is observed among MSE features extracted from the same electrode; however, we may want to keep multiple correlated MSE variables in our model in order to characterize the correlation in spatial patterns of functional brain activity. Therefore, we used ENet, a variation of LASSO, to account for this grouping effect (Zou and Hastie, 2005). Similar to Equation (3), ENet is formulated with a penalty term but in a different format:

$$\lambda \sum\_{j=1}^{p} \left[ (1 - \alpha) \|\beta\_j\| + \alpha |\beta\_j| \right],\tag{4}$$

where α is a trade-off parameter that controls the balance between L1-norm and L2-norm. As α approaches 1, the sparsity of solution will increase such that α = 1 is equivalent to LASSO. On the other hand, α = 0 is equivalent to ridge regression. As α approaches 0, the algorithm tends to encourage group selection of correlated features and stabilize the solution path. In our study, we choose the α = 0.7 for ENet as a empirical choice.

#### 2.4.2. CCA Between MSE and Cognitive Declines

In our study, we used CCA for analyzing the multivariate correlation patterns between MSE features and cognitive decline symptoms related to dementia. The NPI scores includes 12 symptoms: delusions (DEL), hallucinations (HAL), agitation (AG), dysphoria (DEP), anxiety (ANX), apathy (APA), irritability (IRR), euphoria (EUP), disinhibition (DIS), aberrant motor behavior (ABE), night-time behavior disturbances (NIG), and appetite and eating abnormalities (APP). CCA (Hotelling, 1936) is a multivariate analysis approach for finding the relationship

between two sets of variables, **X** and **Y**, with the objective to maximize the Pearson correlation based on projections on new subspaces of **X** and **Y**. **Figure 2** illustrates the cencept. The new feature space is constructed by canonical variables set **U** and **V**, which correspond to original MSE and symptoms rating scales. CCA is formulated as follows:

features as set X and the scores of 12 symptoms from NPI scale as set Y.

$$\underset{\mathbf{u}\in\mathbb{R}^{\mathcal{P}},\mathbf{v}\in\mathbb{R}^{\mathcal{q}}}{\arg\max}\quad\frac{\mathbf{u}^{T}\mathbf{X}^{T}\mathbf{Y}\mathbf{v}}{\sqrt{(\mathbf{u}^{T}\mathbf{X}^{T}\mathbf{X}\mathbf{u})(\mathbf{v}^{T}\mathbf{Y}^{T}\mathbf{Y}\mathbf{v})}},\tag{5}$$

where **X** is a n × p matrix that represents n samples in pdimensional space; **Y** is a n×q matrix that represents n samples in q-dimensional space; **X** and **Y** are two sets of paired variables that correspond to n samples. This problem is solved as a generalized eigen-decomposition problem.

#### 2.4.3. Model Validation and Biomarker Identification

The evaluation of overall performance uses the following three metrics: (1) accuracy indicates the ratio of correctly classified patients in the entire sample; (2) sensitivity indicates the ratio of correctly identified AD patients; and (3) specificity indicates the ratio of correctly identified normal controls, defined as follows:

$$\text{Accuracy} = \frac{\text{TN} + \text{TP}}{\text{TN} + \text{TP} + \text{FP} + \text{FN}},\tag{6}$$

$$\text{Sensitivity} = \frac{\text{TP}}{\text{TP} + \text{FN}},\tag{7}$$

$$\text{Specificity} = \frac{\text{TN}}{\text{TN} + \text{FP}},\tag{8}$$

where TP = true positive, TN = true negative, FP = false positive, and FN = false negative. In particular, normal control group is treated as the negative class in the classification task of this study. If two groups are both AD patients, the less severe group is defined as the negative class. We reported the accuracy of both training and test set to show the potential risks of overfitting, indicated by the gap between training and testing accuracy.

The classifier will be impacted by the imbalanced data during training phase, and the trained model is usually more biased to the majority class. The Receiver Operating Characteristic (ROC) analysis is thus employed for performance evaluation. The area under ROC curve (or AUC) is used as an alternative metric without bias from the selection of threshold parameter (e.g., cut-point) in binary classification of logistic regression.

In addition, we use a leave-one-subject-out cross-validation design to minimize the bias introduced by sample variability. That is, the generalization error is estimated by leaving out samples collected from in the three sessions of one participant for testing and training the model on remaining samples. Validation repeats for all participants as testing samples. Furthermore, the importance of EEG biomarkers was assessed by overall selection frequency in all iterations.

#### 3. RESULTS

#### 3.1. Classification for AD Severity

**Table 1** presents the classification performances of three algorithms. AD groups are considered as the target class. ENet classifier with α = 0.7 (Enet 0.7) yields the best accuracy for classification tasks of HC vs. AD2 and AD1 vs. AD2, and LASSO classifier performs better in discriminating HC vs. AD1, AD1 vs. AD3, and AD2 vs. AD3. Neither model is able to classify AD1 vs. AD2 given the low specificity, although the AUC achieved ∼ 0.7. From the feature selection perspective, grouping effect is accounted for in ENet, which allows for multiple selection among correlated MSE features. This property, considering the high correlation among EEG biomarkers, may better describe the topological patterns for brain activity. Finally, LR with no regularization performed 100% accurate for the training tasks, but the model has poor generalizability because of low test set accuracy and AUC, which indicates the over-fitting issue. All the above results show that the regularized learning methods provide insights about EEG biomarkers with lower risks of over-fitting than LR models.

#### 3.2. Multivariate Correlation Between MSE and Cognitive Declines

The structure coefficients in canonical variables for all channels and symptoms are presented in **Figure 3**. These structure coefficients can be interpreted as the loadings of each original variables (MSE features and cognitive declines) projected into the canonical space. In **Figure 3**, the left panel shows the coefficients of symptoms and right panel shows the absolute values of coefficients for MSE features across all channels. These figures describe how the MSE features and cognitive symptoms contributed to all canonical variables, which suggests a multivariate correlation pattern between clinician's rating and functional brain activity. Our study focus on canonical variables 1–6, since they have higher coefficients of MSE features. For example, in canonical variable 1, the combination of symptoms IRR, DIS, ANX and ABE is associated with channels P3, O1, O2 and central electrodes in short-term complexity, but associated with the frontal area in long-term complexity. In canonical variable 2, the combination of DIS, DEL and APA is associated with central-frontal region. In canonical variable 3, the combination of symptoms DEP, ANX, AG, APA and APP with is associated with frontal region. Canonical variables 4 and 5 present a similar correlation pattern between symptoms ANX, EUP, and APP, and frontal region, but with different signs (positive and negative). Canonical variable 6 presents a positive functional correlation between temporal regions with HAL and AG, but a negative correlation with DIS and IRR. We noted that most significant coefficients are assigned to low time factors (1– 4), while very few non-zero coefficients are distributed in frontal regions for higher (5–8) time factors. In addition, canonical variables 7–12 yield relatively small coefficients comparing to canonical variables 1–6.

#### 3.3. Topological Patterns of EEG Changes Associated With AD Severity

**Figures 4**, **5** display the frequency distribution of selected MSE features in all EEG channels across the brain regions. In the classification tasks of HC vs. AD1 and HC vs. AD2, the selected MSE features were concentrated in the low scale factors (1-4) and distributed diversely from frontal-central to temporal and occipital regions. In contrast, in the classification task of HC vs. AD3, a relatively consistent selection of channels was shown across subjects, mainly in channels T5, T6, O1, and O2.

#### 4. DISCUSSION

#### 4.1. Classification Results