Neurophysiologic Complexity in Children Increases with Developmental Age and Is Reduced by General Anesthesia

This was a prospective, single-center, cross-sectional, observational study assessing cortical complexity in children undergoing general anesthesia for elective outpatient surgery. The study was approved by the University of Michigan Medical School Institutional Review Board (Ann Arbor, Michigan; approval No. HUM00142298). After careful discussion, written informed consent by parents/guardians and verbal or written assent by pediatric patients were obtained before study enrollment. All study operations were conducted at C.S. Mott Children’s Hospital, Michigan Medicine, University of Michigan. Recruitment took place from November 2018 to March 2020. This study adheres to the Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) guidelines.²⁴  Statistical analysis plans for the primary and secondary outcomes and meaningful effect size were defined a priori.

Pediatric patients aged 8 to 16 yr old, with American Society of Anesthesiologists physical status I or II, and scheduled for outpatient elective surgery with a halogenated ether as the primary anesthetic were eligible for study enrollment. Exclusion criteria included a patient history of seizure disorder, developmental delay, neurologic disease, current use of stimulant medications (e.g., amphetamine, dextroamphetamine), surgery above the neck (which might preclude neurophysiologic monitoring), history or suspicion of a difficult airway, physical characteristics that prevent electrode contact with scalp, enrollment in conflicting research protocol, or where English was not the primary language.

The goal of this study was to determine the changes in cortical complexity across age in the preoperative, baseline state of consciousness and through the perioperative period. As such, no protocol was implemented for altering patient care. Clinical teams provided care as indicated, based on standard perioperative protocols, and were blinded to the electroencephalogram data to prevent additional sources of bias.

The electroencephalogram was recorded from 16 Ag/AgCl scalp electrodes using a wireless electroencephalographic system (mobile 128 system, Cognionics Inc., USA) and applied based on the international 10-20 system (Supplemental Digital Content, fig. 1, https://links.lww.com/ALN/C668). Head circumference was measured to ensure proper cap size (EASYCAP, Germany). Cap placement was based on electrode position Cz, which was localized to half of the distance between the nasion and inion and the preauricular notch measurements. Data recordings were sampled at 500 samples/s and referenced to the mastoid. Electrode impedances were continuously monitored and maintained at less than 100 kΩ per manufacturer recommendations. The raw electroencephalogram signals were exported into MATLAB (version 2019b; MathWorks, Inc., USA) and downsampled to 250 Hz. The 60-Hz power-line interference, if present, was removed using a multitaper regression technique and Thomas F-statistics implemented in CleanLine plugin for EEGLAB toolbox.²⁵ 

Electroencephalogram epochs selected for data analysis are shown in figure 1. Baseline (n = 50) electroencephalogram data (range, 3.27 to 5.97 min) were recorded in the preoperative eyes-closed resting state. Premedication electroencephalogram data for all subjects were extracted from a 2-min segment that was recorded within 5 min before and as close as possible to the induction of general anesthesia. If premedication was administered, the selected segment was at least 2 min after intravenous administration (n = 18) or 20 min after oral administration (n = 3). Maintenance (n = 49) electroencephalogram data (5 min) were recorded during the maintenance of general anesthesia and approximately halfway between surgical incision and cessation of the anesthetic maintenance agent. The specific time epoch chosen was based on the following additional criteria: constant age-adjusted minimum alveolar concentration value (greater than 0.7 and less than ±0.1% change) and electroencephalogram data suitable for analysis (i.e., free from artifact by visual inspection). During this maintenance period, the mean age-adjusted minimum alveolar concentration value was 1.26 ± 0.35. Recovery (n = 45) electroencephalogram data (range, 1.18 to 5.75 min) were recorded in a postoperative, eyes-closed resting state after clinical recovery determined by achieving a University of Michigan Sedation Scale score of 0 to 1.²⁶  For the maintenance phase, one patient was excluded because of a minimum alveolar concentration below selection criteria, and five patients in the recovery phase were excluded because of data loss (e.g., associated with emergence delirium).

For each epoch, the electroencephalogram signals were preprocessed as previously described.²⁷  First, bad channels and noisy time segments with obvious artifacts were rejected by visual inspection; the data after artifact removal had an average length (mean ± SD) of 4.88 ± 0.39 min with 15 to 16 channels for the baseline, 5 ± 0 min with 14 to 16 channels for the maintenance, and 3.06 ± 1.21 min with 11 to 16 channels for the recovery epochs. Second, the signals were detrended using a local linear regression method with a 10-s window at a 5-s step size in the Chronux analysis toolbox,²⁸  and low-pass-filtered at 50 Hz using the eegfiltnew function in the EEGLAB toolbox.²⁵  Third, the signals underwent independent component analysis using the extended-Infomax algorithm in the EEGLAB toolbox.²⁵  Independent components representing cardiac, eye, muscle, or other transient artifacts were identified and removed using visual inspection of the time-domain waveform, power spectrum, and spatial scalp topography. The number of independent components removed were (median [25th, 75th]) 3 [2, 4] for the baseline epochs, 0 [0, 0] for the maintenance epochs, and 4 [2, 5] for the recovery epochs.

The power spectrogram was estimated using the multitaper method in the Chronux analysis toolbox,²⁸  with a window length of 4 s with 50% overlap, a time-bandwidth product of 2, and a number of tapers of 3; the estimates were then averaged over all available windows to obtain the averaged power spectrum for each available channel. The normalized power spectrum was further calculated as the absolute power spectrum divided by the total power at 0 to 50 Hz. Electroencephalogram power was calculated for delta (0.5 to 4 Hz), theta (4 to 7 Hz), alpha (7 to 13 Hz), beta (13 to 25 Hz), and gamma (25 to 50 Hz) from both absolute and normalized power spectra. The topographic maps of group-level spectral power across all subjects were constructed for each frequency band using the topoplot function in the EEGLAB toolbox.²⁵ 

We used the Lempel–Ziv algorithm to determine the complexity of cortical dynamics across different states. Lempel–Ziv complexity is a method of symbolic sequence analysis¹⁸  that serves as a surrogate measure of temporal and spatiotemporal complexity of brain activity.^{14–16,21,23}  In this study, we assessed both spatiotemporal complexity across multiple channels and temporal complexity in individual channels.

The spatiotemporal complexity was measured as previously described.²³  Specifically, the instantaneous amplitude was estimated by applying the Hilbert transform, which was then segmented into 4-s windows with 50% overlap (additional analysis was performed to test the effect of window length; Supplemental Digital Content, fig. 2, https://links.lww.com/ALN/C668). The data were then converted to a binary value using the mean value as the threshold for each channel. The data window was then converted into a binary matrix in which rows represent channels and columns represent time points. The complexity of the matrix was assessed by spatiotemporal Lempel–Ziv algorithm, which reflects the number of different spatial patterns across time. If the matrix is random, the spatiotemporal Lempel–Ziv complexity tends to be high; if channels behave similarly (or identically), Lempel–Ziv complexity is low. Because the Lempel–Ziv complexity value for a sequence of fixed length is maximal if it is entirely random, we normalized the raw spatiotemporal Lempel–Ziv complexity by the mean of those from n = 50 surrogate data generated by randomly shuffling the original spatial order for each time point; thus, the resultant spatiotemporal Lempel–Ziv complexity values range from 0 to 1. To test whether the Lempel–Ziv complexity is accounted for by spectral changes, we normalized the spatiotemporal Lempel–Ziv complexity by the mean of those from a surrogate time series (n = 50) generated through phase randomization that preserves the spectral profiles of the signal for each channel. If the complexity change is entirely due to spectral changes, the difference in the Lempel–Ziv complexity values across the states will be completely preserved in the Lempel–Ziv complexity values from the surrogate data; thus, the normalized spatiotemporal Lempel–Ziv complexity will be close to 1 and equal across the states. If the change in complexity is not due to spectral changes, the alterations of Lempel–Ziv complexity values from the surrogate data will be distinct from those of spatiotemporal Lempel–Ziv complexity, and the normalized spatiotemporal Lempel–Ziv complexity will reflect the signal diversity beyond the spectral changes. Last, the spatiotemporal Lempel–Ziv complexity and normalized spatiotemporal Lempel–Ziv complexity values were averaged across all available windows as the final estimate for each state and participant.

To measure individual-channel temporal complexity, the number of different temporal patterns in each individual channel was analyzed. These data were then normalized by the mean of those from n = 50 surrogate data generated by randomly shuffling the temporal order for each channel to obtain temporal Lempel–Ziv complexity. We further tested whether the difference in the Lempel–Ziv complexity across states is due to spectral changes by comparing the temporal Lempel–Ziv complexity with the mean of those from n = 50 surrogate time series generated through phase randomization. The temporal Lempel–Ziv complexity and normalized temporal Lempel–Ziv complexity values were averaged across all available windows as the final estimate for each channel, state, and participant. The topographic maps of group-level temporal complexity across all subjects were constructed using the topoplot function in the EEGLAB toolbox.²⁵  For each epoch and subject, the mean temporal complexity was further obtained by averaging the temporal Lempel–Ziv complexity (and normalized temporal Lempel–Ziv complexity) in prefrontal (Fp1, Fp2), frontal (F5, F6, Fz), central (C3, C4, Cz), parietal (P5, P6, Pz), and occipital (O1, O2) regions, as well as all further channels across the scalp.

To assess the general relationship between brain network maturation and complexity, a large-scale functional brain network was simulated using a coupled Stuart–Landau model implemented in a neuroanatomically informed scaffolding derived from human diffusion tensor imaging. Complexity values of oscillations in the simulated brain network models were compared to assess the effect of developmental changes in the network hub structure. We chose the Stuart–Landau model because it can replicate the oscillatory dynamics of different types of brain signals.^29–32  The coupled Stuart–Landau model is defined as follows:

Here, the complex variable $zj(t)$ determined a state of the node (brain region) j at time t, $j=1,2,…,N$⁠. The anatomical structure A was acquired from group-averaged diffusion tensor imaging with n = 82 nodes.³³  The A_jk is determined by the connection weight between brain regions j and k. We modulated A_jk to simulate brain network maturation based on previous studies demonstrating that hub structure is associated with the developmental age.^9–11,34  Therefore, we additionally used weak (one tenth lower in connection strength compared with the ten strongest hubs) and strong (five times larger in connection strengths compared with the ten strongest hubs) hub structures to simulate the effect of brain network maturation. The brain network model simulation described below was performed with three different brain anatomical structures. The dynamics of the oscillator settle on a limit cycle if $λj>0$ and on a stable focus if $λj<0$⁠. We modulated the $λj=λ$ from –3 to 3 with $δλ=0.2$⁠. The $ωj$ is an initial angular natural frequency of each oscillator j. To simplify the model, we used a Gaussian distribution for natural frequency with a mean frequency of 10 Hz and SD of 0.3 Hz to simulate the narrow bandwidth of human electroencephalogram activity in the eyes-closed resting state.^30–32,35  We controlled a coupling term $Kjk=K$ between oscillators j and k from 0 to 0.5 with $δK=0.01$⁠, which determines the global connection strength among brain regions. To make the model more realistic, we introduced a time delay between brain regions, $τjk=Djk/s$⁠, with the average speed of axons in brain areas, s = 7 ms³⁶  and the distance D_jk between brain regions. The brain region j receives input from connected region k after the time delay $τjk$⁠. The model results are not qualitatively different if the time delay is smaller than a quarter of the period of oscillation.³¹  A Gaussian white noise $ξj(t)$ for each region was added with the SD $β=0.005$⁠. We numerically solved the differential equations of the Stuart–Landau model using the Stratonovich–Heun method with 1,000 discretization steps. The first 10 s of the generated signals were discarded, and the last 50 s were used for the analysis of each simulation. Each brain region generates its own spontaneous oscillatory dynamics in a network at each bifurcation parameter $λ$ and coupling strength K for one simulation. The simulation was repeated 50 times with different frequency configurations to obtain statistical robustness.

Among the simulated brain signals at various bifurcation parameter $λ$ and coupling strength K, we selected brain states at certain parameter sets that can represent conscious states. A variance of the level of global synchronization in a network, termed the pair correlation function, was calculated, and the state at certain coupling strength K with the largest pair correlation function was chosen as the state that can represent the conscious state for each bifurcation parameter $λ$⁠.³²  The instantaneous global synchronization level $r(t)$ at time t was calculated using phase difference $Δθjk(t)=θj(t)−θk(t)$ at each coupling strength K.

Here $r(t)=1$ if all phases are equal, but $r(t)$ is nearly 0 if all phases are randomly distributed.

Global pulsatile stimuli to the whole brain network were induced at the states with the largest pair correlation function to observe the complexity from the response to the stimuli. The coupled Stuart–Landau model with the stimulation term $u(t)$ is as follows:

Here p is the strength of the stimulus during a period $T= t2−t1$⁠. We fixed $p=30$ and set duration of the stimulus as T = 100 ms. We induced the stimulus at 10 different random timings for one iteration. Each stimulus was applied independently to generated signals within one frequency configuration.

A significant response was calculated by comparing the instantaneous amplitude values $I(t)$ before and after the stimuli. For each iteration, the $I(t)$ of each node j after stimuli was normalized by the mean and SD of the baseline amplitude values of node j. Baseline values were obtained by using a total of 100 s, consisting of 10 trials of a 10-s prestimulus segment for each iteration. We considered the one tail $(1−α)∗100th$ quantile with $α=0.05$ as a significantly increased amplitude. A perturbation response (represented in equations as PR) of node j at time t was defined in a binary fashion: $PRj(t)=1$⁠, for the significantly increased amplitude for node j, and $PRj(t)=0$⁠, otherwise. The complexity was calculated by measuring Lempel–Ziv complexity $LZc(t)$ of the $PR(t)$ over 1 s after the stimuli. The Lempel–Ziv complexity was calculated for brain network models with low, intermediate, and high hub structure connectivity.

To assess the effects of premedication on cortical complexity and spectral properties, the resting electroencephalogram data during the preinduction period were analyzed. The electroencephalogram signals were preprocessed as described for the primary analysis and, after preprocessing, 12 to 16 channels remained for n = 50 subjects. The spectral power and cortical complexity measures were calculated as described for the complexity analysis.

Statistical analyses were performed using MATLAB, and the data were tested for normality of distribution by Lilliefors corrected Kolmogorov–Smirnov tests. Because the null hypothesis of normality of distribution was rejected in some of the data sets (P < 0.05), the two-sided Wilcoxon signed rank test was used to compare the electroencephalogram measure (spectral power and cortical complexity) across baseline, maintenance, and recovery. With Bonferroni correction, P < 0.017 (0.05/3; 3 pairs) was considered statistically significant. Spearman correlation was used to investigate the relationship between age and each electroencephalogram measure. Univariate linear regressions, with age as the independent variable and complexity measures as the dependent exposure variables, were used to determine the slope [95% CI] of the associations, with R² used to assess goodness-of-fit. For the effect of premedication, the two-sided Wilcoxon rank sum test was used to compare the electroencephalogram measure (median [25th, 75th]) between subjects without (n = 29) and with premedication (n = 21). A Kruskal–Wallis test was used to calculate the statistical differences (median [25th, 75th]) across simulated brain networks, and P < 0.05 was considered statistically significant.

A sample size of 41 achieved 80% power to detect a change in slope between complexity and age from 0.000 under the null hypothesis to 0.011 under the alternative hypothesis when the statistical hypothesis is two-sided, the significance level is 0.050, the variance of age is 1.0816, the variance of complexity is 0.0009, the error variance of complexity is 0.5, and the correlation between observations within an individual is 0.25, assuming a compound symmetry correlation structure. This was a conservative estimate based on a previous study of complexity in healthy control patients of similar ages.¹⁹ 

In total, 175 children were screened for study eligibility. Of the 175 who met the study criteria, 36 children declined participation, 52 were excluded because of operative time change or time constraints, 37 because of enrollment in another study, 21 because of research staff availability, 10 because of cancellation or no show for surgery, 3 because of changes in anesthetic plan, and 2 because of technical complications, leaving 50 participants who completed the study and were included for analysis. The participant demographic, anesthetic, and surgical characteristics are presented in table 1. Five participants had incomplete data for race.

Age-related changes in spectral properties are shown in figure 2. Baseline normalized gamma power in the frontal and parietal region and beta in the parietal region were correlated with age (gamma: frontal r = 0.30, P = 0.034, and parietal r = 0.31, P = 0.027; beta: parietal r = 0.40, P = 0.003), whereas baseline theta was inversely correlated with age (r = –0.45, P = 0.001 and r = –0.57, P < 0.001 for frontal and parietal regions, respectively). Total power decreased with age, with a peak observed at 8 yr, and subsequent decreases were present during the baseline (slope [95% CI], –0.63 [–0.91, –0.35]; P < 0.001; R² = 0.29), maintenance (–0.65 [–0.94, –0.35]; P < 0.001; R² = 0.28), and recovery phases (–1.2 [–1.6, –0.79]; P < 0.001; R² = 0.43); Supplemental Digital Content, fig. 3, https://links.lww.com/ALN/C668). Total power also tended to decrease across regions and frequency bands with increasing age (Supplemental Digital Content, fig. 3, https://links.lww.com/ALN/C668). During maintenance, there was an inverse correlation for normalized theta power in the parietal region (r = –0.32, P = 0.027). During recovery, there was an inverse correlation for normalized delta power r = –0.44, P = 0.002 and r = –0.35, P = 0.020 for frontal and parietal regions, respectively). Conversely, higher (beta and gamma) frequencies were positively correlated during recovery (beta: frontal r = 0.45, P = 0.002 and parietal r = 0.57, P < 0.001; gamma: frontal r = 0.51, P < 0.001 and parietal r = 0.53, P < 0.001; fig. 2, A and B).

Group-level changes in spectral properties are shown in figure 3 (and Supplemental Digital Content, fig. 4, https://links.lww.com/ALN/C668). Overall, during general anesthesia, there were global increases in delta power and alpha anteriorization and decreases in the higher frequency band power (beta and gamma; fig. 3, A and B). In terms of total power, general anesthesia induced an overall increase in power (except for gamma), which returned to baseline level after recovery for most frequencies (except alpha; Supplemental Digital Content, fig. 4, https://links.lww.com/ALN/C668).

In the preanesthetic baseline state of consciousness, there was a positive correlation with age and spatiotemporal Lempel–Ziv complexity (r = 0.41, P = 0.003), as well as age and averaged temporal Lempel–Ziv complexity (r = 0.39, P = 0.005). Likewise, there was a significant linear association with age for both measures (slope [95% CI], 0.010 [0.004, 0.016], P < 0.001, R² = 0.20; and 0.007 [0.003, 0.011], P < 0.001, R² = 0.20, respectively; fig. 4, A and B). When these data were normalized to account for spectral effects, there was no correlation for the normalized spatiotemporal Lempel–Ziv complexity or normalized averaged temporal Lempel–Ziv complexity (r = 0.02, P = 0.916 and r = 0.08, P = 0.572). Likewise, the linear association was not statistically significant (slope [95% CI], 0.001 [–0.004, 0.005], P = 0.737, R² < 0.01, and 0.002 [–0.001, 0.004], P = 0.234, R² = 0.03, respectively; fig. 4, C and D). Regional analysis of temporal Lempel–Ziv complexity and normalized temporal Lempel–Ziv complexity mirrored the findings above for the whole brain average (data not shown). The results were qualitatively similar when we applied robustness testing by varying the duration of window analysis (Supplemental Digital Content, fig. 2A, https://links.lww.com/ALN/C668).

Because our study design presumed but did not measure brain network development in the cohort of participants, we employed a principled approach to assess complexity in developing networks, as defined by networks that had progressively increased hub structures. After a simulated perturbation, complexity was measured in computational brain models informed by human neuroanatomy and compared across weak, intermediate, and strong hub structures (as quantitatively defined under “Materials and Methods”). Lempel–Ziv complexity significantly increased across simulated brain models with increasing strength of connectivity in the network hub structure (P < 0.001; fig. 5).

Midazolam premedication was not associated with a change in cortical complexity when measured across the entire brain (fig. 6A). However, regional increases in complexity were found for both the temporal Lempel–Ziv complexity and normalized temporal Lempel–Ziv complexity in the central region (median [25th, 75th], 0.384 [0.352, 0.418] vs. 0.434 [0.414, 0.452], P = 0.002; 0.915 [0.898, 0.933] vs. 0.945 [0.915, 0.961], P = 0.014, respectively) and parietal region (0.377 [0.365, 0.401] vs. 0.420 [0.412, 0.456], P < 0.001; 0.923 [0.900, 0.939] vs. 0.938 [0.915, 0.958], P = 0.031, respectively; fig. 6, B and C). Midazolam premedication was associated with frequency- and region-specific changes in normalized electroencephalogram power (Supplemental Digital Content, fig. 5, https://links.lww.com/ALN/C668). In all brain regions, there was a decrease in theta and increase in beta power, as well as a decrease in frontal and prefrontal alpha power (P < 0.05).

During the stable maintenance phase of general anesthesia, cortical complexity decreased from baseline values for both complexity measures tested, with the greatest decrease from baseline in spatiotemporal Lempel–Ziv complexity (median [25th, 75th] for baseline 0.660 [0.620, 0.690] vs. maintenance 0.459 [0.402, 0.527], P < 0.001; fig. 7; Supplemental Digital Content, fig. 6, https://links.lww.com/ALN/C668). Changes in complexity were independent of spectral power during the maintenance phase of general anesthesia. During recovery, spatiotemporal Lempel–Ziv complexity exceeded the baseline level (0.704 [0.642, 0.745], P < 0.001); however, after controlling for spectral changes, the normalized spatiotemporal Lempel–Ziv complexity during recovery was significantly decreased compared with baseline (baseline 0.913 [0.887, 0.923], vs. recovery 0.873 [0.840, 0.902], P < 0.001). Similar results were found for the average temporal Lempel–Ziv complexity analysis, with the exception of the normalized average temporal Lempel–Ziv complexity when the recovery period returned to baseline levels (0.960 [0.949, 0.966] vs. 0.957 [0.945, 0.965], P = 0.404, respectively). Regional changes in normalized temporal Lempel–Ziv complexity results were similar to data from the whole brain average (Supplemental Digital Content, fig. 7, https://links.lww.com/ALN/C668). Overall, these results were consistent when we varied the analysis window duration (Supplemental Digital Content, fig. 2B, https://links.lww.com/ALN/C668).

In a post hoc analysis of changes in spatiotemporal Lempel–Ziv complexity during the maintenance epoch with (n = 11) or without (n = 38) the use of nitrous oxide, there was no difference except that the normalized spatiotemporal Lempel–Ziv complexity was higher in the nitrous oxide group (P = 0.045, Wilcoxon rank sum test; data not shown). In an additional post hoc analysis of the association of spatiotemporal Lempel–Ziv complexity and the average age-adjusted minimum alveolar concentration (maintenance epoch) with or without the use of nitrous oxide, we found a negative correlation (with nitrous oxide r = –0.64, P = 0.040, without r = –0.36, P = 0.028; Spearman correlation; data not shown), but this was not significant after controlling for spectral changes (with nitrous oxide r = –0.37, P = 0.261, without r = –0.02, P = 0.916; similar results were found for average temporal Lempel–Ziv complexity; data not shown).

This study of children undergoing surgical anesthesia tested the hypothesis that electroencephalogram cortical complexity would increase with developmental age and decrease with general anesthesia. Age was positively correlated with cortical complexity during baseline recordings, and further analyses revealed that this was attributable to spectral changes. We supported this empirical finding through a principled approach by investigating complexity in simulated brain networks representative of development. This analysis demonstrated increased complexity with increasing connectivity strength of the network hub structure. During anesthetic state transitions, we found that cortical complexity decreased during the maintenance phase of general anesthesia compared with the eyes-closed baseline. Furthermore, after recovery of consciousness, normalized spatiotemporal complexity remained reduced compared with baseline. Overall, we found that age and anesthetic-mediated perturbations in the level of consciousness were associated with changes in cortical complexity, with age-related changes likely resulting from spectral changes and, as suggested by our simulation results, evolving functional architecture.

The hypothesis that cortical complexity would increase as a function of age was informed by the developmental maturation occurring through this period. Structural and functional changes, such as increased white matter density, long-range connections,^1,37  and network connectivity,^1,38  could yield a less reducible and more “complex” brain network. Additionally, studies of childhood and adolescent brain network development have demonstrated increased differentiation, specialization, and organization of network function,^1,9,38  potentially leading to a greater diversity of neural oscillatory activity and measurable using complexity algorithms. Indeed, previous studies of resting state magnetoencephalogram data analyzed using the Lempel–Ziv complexity algorithm demonstrated a significant increase in complexity during childhood and adolescence that continued throughout the lifespan, peaking in the sixth decade.^19,20  Our results are consistent with this, but the effect is attributable to changes in spectral properties, which were not controlled for in the studies by Fernández et al.^19,20  The lack of an age effect in complexity beyond spectral properties is somewhat surprising given the magnitude of brain network change occurring during this period but is possibly accounted for by several factors. First, we did not control for developmental heterogeneity. It also remains possible that the normalized measure of complexity is not sensitive enough to measure a developmental difference but requires significantly larger perturbations (e.g., anesthetic state transitions). The robust finding of significant decreases during general anesthesia, even when controlling for spectral changes, supports this interpretation. It should also be considered that despite the large developmental changes occurring in children 8 to 16 yr old, it remains possible that this age range may not capture the structural and functional changes required for discernable differences in normalized complexity. Additionally, a developing network may give rise to both spectral and complexity changes that temporally align.³⁹  Our computational brain model findings that networks with stronger hub structures yield a more complex network further support the empirical age-related results. Moreover, previous work from our group has shown that network hub structure is disrupted with propofol administration and consistent with reductions in complexity with general anesthesia.¹² 

Complexity analysis using the Lempel–Ziv algorithm yielded dynamic results. Preoperative midazolam administration increased complexity in the central and parietal regions, whereas complexity was globally reduced during general anesthesia. Previous studies from our group demonstrated a similar increase in complexity with subanesthetic ketamine administration and reduction during ketamine anesthesia.²³  Although the mechanism for increased complexity during subanesthetic states remains unknown, it has been shown that γ-aminobutyric acid type A receptor agonists, such as propofol, can induce paradoxical excitation.⁴⁰  Additionally, benzodiazepines have been shown to increase beta activity in the central and parietal regions consistent with the increased beta activity we observed.^41,42  In support, a post hoc analysis of temporal complexity and beta power demonstrated a positive correlation and after controlling for spectral changes, was mitigated in the no-premedication group, but remained significant for the midazolam group (Supplemental Digital Content, fig. 8, https://links.lww.com/ALN/C668). This suggests an association between complexity and beta power, which could be a concurrent effect of midazolam on the electroencephalogram. However, a recent study of subanesthetic nitrous oxide administration demonstrated decreased electroencephalographic complexity (Shannon complexity).¹³  The effect of various subanesthetic medications on cortical complexity in children remains an open question. In addition, there may be region specific changes (fig. 7B), but further investigation including high-density recordings would be required to address this fully.

Leading theories of consciousness posit that the diversity of information in both space and time is important for conscious processing.^43,44  The biologic interpretation of decreases in cortical complexity may represent reduced information sharing or processing. In this study, decreases in cortical complexity induced by general anesthesia were similar to what has been shown in animal models^16,45  and adult patients.^14,17,21,23  The significant decrease during general anesthesia of cortical complexity, independent of spectral changes, suggests that measures of complexity may be helpful in the development of real-time monitoring modalities across the adult and pediatric populations. It is important to note that complexity measures in this study reflect the level of consciousness despite the age-related changes in power occurring during this period of massive neurodevelopment. When developing monitoring systems, measures of the conscious state need to be robust, have high fidelity, and be accurate across the population. In our study, there was a decrease in all children with both spatiotemporal and averaged temporal complexity; however, when the data were normalized, nine children in the spatiotemporal and seven in the averaged temporal complexity analysis groups had actual increases in complexity during the maintenance of general anesthesia (Supplemental Digital Content, fig. 6, https://links.lww.com/ALN/C668). Additional post hoc analysis of age, anesthetic dose, or channel rejection did not account for this increase (data not shown). Recent work by our group in an animal model artificially dissociated the observed level of consciousness from various electroencephalographic measures, including complexity.¹⁶  However, as noted in this study, just because state-related and neurophysiologic dynamic changes can be dissociated in the laboratory setting does not mean that they are not correlated during spontaneous physiologic, pharmacologic, or pathologic state transitions. These findings suggest that a single quantitative surrogate of consciousness using normalized complexity analysis may be insufficient to reliably determine the level of consciousness in all patients and future developmental studies should also address combinatorial surrogate measures. In addition, the finding of decreased normalized spatiotemporal complexity after clinical recovery is not surprising considering that functional brain network recovery extends beyond the initial return of consciousness.^46,47 

In an effort to confirm that our findings were consistent with the spectral characteristics across development shown in previous works, we analyzed both the absolute and normalized power of the electroencephalogram in this cohort. Similar to previous studies of school-age and adolescent children, we observed an age-dependent decrease in absolute electroencephalogram power over most brain regions during general anesthesia.^5,8  In addition, when the power data were normalized, regional differences were aligned with developmental structural and functional changes occurring during this period.^1,48  Specifically, there was a general shift toward faster oscillatory activity (i.e., decreased theta and increased beta and gamma power).^3,4  When we analyzed the spectral data during general anesthesia, we found increased absolute power, increased delta power, decreased gamma power, and anteriorization of the alpha activity, consistent with what has been previously shown during this developmental period.^5,49,50 

There are several limitations of this investigation. First, this was an observational study without a protocol for a standard anesthetic regimen, and although this was pragmatic, we are unable to dissect the role of individual anesthetic agents. We performed a post hoc analysis to investigate the role of nitrous oxide and minimum alveolar concentration values influence on complexity. Although we are cautious in interpreting these data because of the small subsamples, the findings were consistent with the overall findings. Additionally, patient characteristics that may be present and attributable to heightened preoperative anxiety could also influence complexity and spectral properties. Further, sex differences have been shown to affect cortical complexity with development, and future larger-scale studies should address this.²⁰  When considering the Lempel–Ziv complexity measure, there are several notable strengths: it is a nonparametric measure, can be applied to relatively short data segments, and can be used with nonstationary data. However, because the measure requires binarization of the original signal, there exists the possibility for the loss of information during the analysis as well as influence by slow waves in the electroencephalogram (Supplemental Digital Content, fig. 9, https://links.lww.com/ALN/C668). Additionally, this study was conducted in the perioperative environment and limited by a low electroencephalogram (16 channel) density and, in some instances, a loss of information from an electroencephalogram channel that has the potential to influence the assessment of cortically complexity. However, when we repeated our analysis excluding participants with channel rejection, the results were similar (data not shown).

In summary, cortical complexity increased with developmental age and decreased during general anesthesia. The age-related increase in cortical complexity could be accounted for by spectral properties but, based on modeling data, is consistent with a maturing brain network. By contrast, reductions in complexity during general anesthesia extended beyond spectral properties. Overall, these findings contribute to the growing body of literature on cortical complexity and how this potential neural correlate of consciousness is affected by age and anesthetic state transitions.

The authors thank the children and their families for their participation in this study.

Supported by Foundation for Anesthesia Education and Research (Schaumburg, Illinois) grant No. 20-PAF00823 and by the Department of Anesthesiology, University of Michigan Medical School, Ann Arbor, Michigan.

Dr. Kaplan has received research funding from Lundbeck Biopharmaceutical Company (Copenhagen, Denmark). The other authors declare no competing interests.