Neonatal cortical activity organizes into transient network states that are affected by vigilance states and brain injury

Early neurodevelopment is critically dependent on the structure and dynamics of spontaneous neuronal activity; however, the natural organization of newborn cortical networks is poorly understood. Recent adult studies suggest that spontaneous cortical activity exhibits discrete network states with physiological correlates. Here, we studied newborn cortical activity during sleep using hidden Markov modeling to determine the presence of such discrete neonatal cortical states (NCS) in 107 newborn infants, with 47 of them presenting with a perinatal brain injury. Our results show that neonatal cortical activity organizes into four discrete NCSs that are present in both cardinal sleep states of a newborn infant, active and quiet sleep, respectively. These NCSs exhibit state-specific spectral and functional network characteristics. The sleep states exhibit different NCS dynamics, with quiet sleep presenting higher fronto-temporal activity and a stronger brain-wide neuronal coupling. Brain injury was associated with prolonged lifetimes of the transient NCSs, suggesting lowered dynamics, or flexibility, in the cortical networks. Taken together, the findings suggest that spontaneously occurring transient network states are already present at birth, with significant physiological and pathological correlates; this NCS analysis framework can be fully automatized, and it holds promise for offering an objective, global level measure of early brain function for benchmarking neurodevelopmental or clinical research.


Introduction
Cortical activity exhibits constellations where far apart brain areas are functionally connected to support global functions. These large-scale networks can be characterized by their respective spatiotemporal correlations in locally recorded brain activity (Hallett et al., 2020;Helwegen et al., 2023;Kobeleva et al., 2022;Yrjölä et al., 2022;Zhang et al., 2021). Such cortical activity networks are believed to support a wide range of higher brain functions (Gao et al., 2015;Haque et al., 2020;Rouhinen et al., 2020). At the same time, sleep is an archetypal vigilance state which emerges in the newborns at the first place and holds crucial information about early brain function Tham et al., 2017;Uchitel et al., 2022). Infant sleep architecture is dominated by interplay of two major sleep states: active sleep (AS) and quiet sleep (QS) which are precursors of rapid eye movement (REM) sleep and non-REM sleep Lokhandwala and Spencer, 2022). Recent studies have shown that already starting from preterm age, these sleep states (AS and QS) are associated with distinct functional brain networks for different intrinsic coupling modes (Tokariev et al., 2016b(Tokariev et al., , 2019aYrjölä et al., 2021). Experimental studies also reported distinct neurodevelopmental roles of these states (Knoop et al., 2021): AS is thought to support formation of sensorimotor networks (Del Rio-Bermudez et al., 2017;Del Rio-Bermudez and Blumberg, 2018) whereas QS is important for memory consolidation (Friedrich et al., 2020).
The cortical networks are often characterized as an average connectivity over longer periods of time, i.e. as if it was static functional connectivity (FC) (Hallett et al., 2020). This approach is methodologically convenient and straightforward, however it ignores the well-established rapid network dynamics at a sub-second timescale that are a hallmark of healthy cortical function (Hutchison et al., 2013;Preti et al., 2017;Shriki et al., 2013;Siebenhühner et al., 2016). Such rapidly alternating network reconfigurations can be seen by analyzing bivariate networks in short time windows, or by identifying otherwise independent, statistically discrete brain states (Allen et al., 2014;Khanna et al., 2015;Michel and Koenig, 2018;Vidaurre et al., 2016). Several recent studies have suggested analytic approaches to measure non-stationary changes in FC, typically in the timescales of seconds or even shorter (Hutchison et al., 2013;Preti et al., 2017). The most common strategy is to use short time windows where recurring brain states are identified in a straightforward manner with clustering algorithms (Allen et al., 2014;Preti et al., 2017). The results of the sliding window are readily challenged by biases that may arise from the selection of window length and shape, which should be optimized a priori to include all temporal scales of interest (Liuzzi et al., 2019;Preti et al., 2017;Shakil et al., 2016). A data-driven approach using hidden Markov modelling (HMM) was recently introduced to characterize brain activity as a dynamic sequence of discrete brain states (Baker et al., 2014), which are first identified from the data using relevant timescales; subsequently, cortical activity characteristics like functional connectivity and/or spectral content can be estimated from these brain states (Baker et al., 2014;Vidaurre et al., 2016Vidaurre et al., , 2017Vidaurre et al., , 2018a. While prior studies have shown the presence of neonatal large-scale cortical networks (Fransson et al., 2007;Gao et al., 2017;Molloy and Saygin, 2022;Tokariev et al., 2019a;Uchitel et al., 2022;Yrjölä et al., 2021Yrjölä et al., , 2022, the properties and physiological correlates of the rapid network dynamics in newborns are not known. The pathophysiological meaning of novel measures of brain function can be assessed from studying their effect on neurological conditions. Perinatal asphyxia (PA) is the globally most common perinatal incidence that may cause neonatal hypoxic-ischemic encephalopathy (HIE), a neurological condition requiring intensive care treatment and posing the infants at risk of neurodevelopmental consequences (Saugstad, 2011). Perinatal asphyxia refers to a temporary compromise in oxygenation and/or circulation near the time of birth, and its clinically most significant consequence is cerebral injury, associated with mild, moderate, or severe clinical signs of HIE; additionally, infant may have a PA exposure without developing clinical signs of HIE (Sarnat, 1976). Clinically, moderate and severe HIE have been the primary interest in developing treatment options while there has been considerable recent interest in understanding the impact of mild HIE and PA exposures on the neurodevelopment (Chalak et al., 2018;Conway et al., 2018;Finder et al., 2020;Murray et al., 2010). Neuroscientifically, this spectrum of conditions gives an opportunity to study relationships between pathophysiological exposures and their effects on brain function. Several effects of moderate and severe HIE are well documented in the literature (Dilena et al., 2021;Lacan et al., 2021;Nevalainen et al., 2017;Temko et al., 2015), however it is poorly understood how the cortical function, i.e., EEG activity, is affected by milder forms of PA exposure.
Here, we aimed to characterize network states and their dynamics in the newborn cortical activity. In particular, we i) estimated the spectral, temporal, and FC properties of discrete cortical states in infants, and ii) extracted the temporal sequence, or dynamics of these cortical states. In search of behavioral and clinical correlates of these newly described network states, we studied the effects of sleep states (behaviour) and neonatal brain injury (clinical). For the latter, we chose to assess EEG recordings during birth asphyxia, which is one of the most common clinical adversities (Chalak et al., 2018;Murray et al., 2010), with well documented effects on the EEG activity (Dilena et al., 2021;McLaren et al., 2019).

Subjects and recordings
During daytime sleep, EEG was recorded from 60 full-term, healthy controls (HC) and 47 hypoxic-ischemic encephalopathy (HIE) infants. These EEG recordings are parts of the data used in our previous studies (Tokariev et al., 2019a(Tokariev et al., , 2019bTuiskula et al., 2022). Further details of the HIE cohort are presented in Supplementary information (see Table S1). All EEG recordings were done at term-equivalent age of 42 ± 1.5 (mean ± std) weeks for both groups at the Children's Central Hospital, Helsinki University Central Hospital. The gestational age of the healthy infants was 41.3 ± 2 (mean ± std) weeks, whereas the gestational age of HIE infants was 40.3 ± 1.3 (mean ± std) weeks. All the procedures in this study were approved by the Ethics Committee of the Helsinki University Central, and informed written consent was received from a guardian.
All the EEG data was initially collected as prospective study for other purposes. The different numbers of electrodes were due to collation of EEG recordings from many different earlier studies in our lab, with some recordings carried out in the hospital EEG unit, and others performed in the research center related to the hospital; all done by the same staff and overall protocols. The present study is a re-use of previously collected datasets (Tokariev et al., 2019b(Tokariev et al., , 2019aTuiskula et al., 2022), hence we consider it retrospective.

EEG preprocessing
All EEG data was visually inspected, and 3-min-long artifact-free EEG segments for each infant were retained for further analyses. The identified artifacts that were excluded belonged to the following categories: excessive electromyography (EMG), mechanical, and instrumental noise. The EEG signals were then band-pass filtered between 0.4 and 22 Hz using a combination of seventh-order low-pass and high-pass Butterworth filters to preserve the physiologically relevant range (Vanhatalo et al., 2005). Filters were applied in both forward and reverse directions to prevent phase lags in recordings. All EEG segments were then downsampled to 100 Hz and re-referenced to the common-average montage (Fig. 1a).

Source reconstruction
Source reconstruction was performed as described elsewhere (Tokariev et al., 2019b). Briefly, the source space was composed of 8014 dipoles orthogonal to the cortical surface. Magnetic resonance imaging (MRI) data from a full-term newborn infant was used to generate a three-shell head model, which included scalp and inner/outer skull surfaces. A magnetic resonance image (MRI) data was obtained using a Philips 3T scanner at Helsinki University Central Hospital (Finland) from a healthy full-term infant. Each slice had a resolution of 240×256 pixels (pixel resolution was 1×1 mm), and the slice thickness was 0.9 mm. From the full image stack, 176 slices covering the cranium were segmented manually into 5 different compartments (scalp, skull, CSF, brain, eyes) by a clinician using FSL software (Odabaee et al., 2014;Sairanen et al., 2017). The tissue conductivities were set to: 0.43 S/m for the scalp, 1.79 S/m for the intracranial volume, and 0.2 S/m for the skull (Despotovic et al., 2013;Odabaee et al., 2014). The forward operator was generated using the symmetric boundary element method (Gramfort et al., 2010), and the inverse operator was computed with dynamic statistical parametric mapping approach (Dale et al., 2000) as they both are implemented in the Brainstorm software (Tadel et al., 2011). Finally, cortical sources were clustered into 58 parcels and cortical signals representing the neural activity of each parcel were calculated as the weighted mean of source signals belonging to that parcel. This parcellation scheme was designed specifically to analyze infant EEG data and is a compromise between spatial resolution and reliability of reconstruction of cortical parcel signals from standard low-density clinical EEG montages. Parcels were symmetric across hemispheres, roundish in shape and comprised 125 ± 26 (mean ± standard deviation) sources each. Weighting of individual sources within the host parcels was aimed to prioritize contribution of the sources with high fidelity (Tokariev et al., 2019b). Parcels were assigned to one of the anatomical brain regions, including frontal, central, occipital, and temporal, based on their anatomical location (Fig. 1a).

Identification of newborn cortical states
It is assumed that the observable data in each time series is generated by switching between several latent states. As shown recently (Quinn et al., 2018;Vidaurre et al., 2018b), HMM is a potentially useful tool for discovering such hidden states in the data. The fundamental assumption of HMM is that the probability of a state being active at the next time point depends solely on the current state. Each state in HMM has its own observational model that represents the link between that state and the observed data (Rabiner, 1986). In order to provide a condensed data representation, the HMM first assigns a state to each point (Rabiner, 1986), and these states can be represented by several measures of the data, including temporal, spatial, and spectral information (Vidaurre et al., 2018a(Vidaurre et al., , 2017. Most importantly, an HMM inferred at the group level will support the extraction of dynamic patterns at the level of individual patients or even trials (Vidaurre et al., 2017), which has made HMM a widely used method in dynamic network analysis (Ahrends et al., 2022;Bai et al., 2021;Gascoyne et al., 2021;Higgins et al., 2021;Hirschmann et al., 2020;Tao et al., 2022;Zhang et al., 2021Zhang et al., , 2022. In the present study, we used time series of 58 cortical parcels, which were reconstructed from the scalp-recorded EEG signals (see above). The HMM approach assumes that observations in each state are drawn from a probabilistic Gaussian model with zero mean and a covariance matrix in the time-delay embedded (TDE) space (Vidaurre et al., 2018b). The HMM in the TDE space (TDE-HMM) is constructed by adding time lagged versions of the data to the original data, resulting for each segment in an extended data matrix with a size of timepoints × (parcels × lags). Next, a principal component analysis is applied to reduce the dimensionality by only including a predefined number of principal components (PCs), transforming the data matrix into a new one with the size of timepoints × PCs. In the current study, time lags varied between -5 and 5 (in samples), corresponding to windows of 100 ms. Only 12 PCs were kept for further analysis, which explained 80% of the variance in the data. Before applying TDE-HMM, the cortical activities in HC cohort were concatenated across subjects and sleep states followed by a normalization (mean = 0, variance = 1). All the computations related to TDE-HMM were done using HMM-MAR toolbox (https://github. com/OHBA-analysis/HMM-MAR) implemented in MATLAB.

State-wise spectral and connectivity assessments
TDE-HMM inference outputs time courses for the brain being in one of the characteristic states (Fig. 1b), which could be used as a mask to compute the following state-specific characteristics: (i) the spectral profiles of activity in each cortical parcel, and (ii) the phasesynchrony between parcel signals using the cross-spectral density.
The phasesynchrony between different brain regions was measured by means of the lagged coherency (Pascual-Marqui, 2007). This measure ignores the spurious connections (Palva et al., 2018;Palva and Palva, 2012) originated from the volume conduction. Moreover, each connectivity matrix, whose entities are the lagged coherence values between all possible pairs of brain regions, was corrected and the same subset of unreliable connections was removed in the parcel space (Tokariev et al., 2019b). The rationale behind this is that for standard clinical low-density EEG caps (in our case 19 channels), the reconstruction of the whole-brain source activities is challenging (Tokariev et al., 2016a) as some cortical areas are in a 'blind zone' for the scalp sensors. The unreliable connections were defined based on the extensive simulations that test head model with defined number of cortical parcels (N = 58) and scalp EEG electrodes (N = 19). More specifically, we generated many cortical networks (reference data) and compared them to their copies after forward-inverse modelling using head model (reconstructed data). Pairwise connections that systematically were poorly reconstructed in simulations (about 32%), were rejected from empirical data (see Tokariev et al., 2019b for details). Thus, unreliable interactions were removed for further analysis.
Cross-spectral density was estimated by means of a state-wise multitaper approach (Vidaurre et al., 2018b). The spectral properties are estimated using several tapers (Thomson, 1982) and Fourier transform. The estimation by this method reduces the variance and bias in comparison with conventional Fourier-transform-based methods (Babadi and Brown, 2014). Here, discrete prolate Slepian sequences with the time-bandwidth product 2, the number of tapers 3, and the length of window 5 seconds were used as tapers of the multitaper method. These were chosen to enable the method to perfectly capture the low-frequency content (~ 0.4 Hz).
To visualize the strongest patterns for each NCS (on the connectivity ring plots) we used Gaussian mixture model (GMM) with two Gaussian curves (see Supplementary for more details). This approach was selected to avoid arbitrary thresholding (Vidaurre et al., 2018b).

Defining the number of newborn cortical states
The robustness of HHM inference output is sensitive to the predefined number of NCS. Thus, to find the solution, we repeated the analysis with various number of states (from 3 to 10) and compared results for stability. Briefly, for each number of NCS, we extracted NCS time courses, by HHM inference, five times on the whole dataset (60 segments for active sleep, 60 segments for quiet sleep in the healthy group; in total 120 concatenated segments) and calculated the correlation between the extracted NCS time courses as a measure of reproducibility. The maximum correlation value was reached with four NCSs, and its robustness to data selection was confirmed by split-half testing (see Supplementary Fig. S1).

Temporal features
State time courses can also be used to look at the temporal dynamics of each NCS. Three common temporal features are fractional occupancy (FO), lifetime (LT), and transition probabilities (Vidaurre et al., 2018b).
FO is a measure that shows the percentage of the total time covered by a specific NCS and LT is the average time during which an NCS is active. The transition probabilities indicate the probability of the transition from one NCS to another one. Together, these metrics summarize the brain dynamics as captured by the HMM (Baker et al., 2014). By splitting the HMM output per segment, we calculated these metrics for each infant individually, during AS and QS.

Complexity analysis of newborn cortical state transitions
The state time courses can be described as a sequence of ordinal patterns. The relative occurrence of each of these patterns can be calculated by permutation entropy (PE) (Bandt and Pompe, 2002;Olofsen et al., 2008). PE is a measure to determine the complexity of a sequence. To calculate PE values, the sequence of brain states was reduced to its transitioning sequence (e.g., sequence 1112233341144 is reduced to 123414) because the between-state transitions are of interest, not self-state transitions. Since the length of transitioning sequences is different among EEG segments and to avoid the bias caused by the length, only the first 400 transitions were considered for each EEG segment. This value equals the shortest sequence length among all EEG segments.
The permutation entropy can be calculated as: where p(π i ) is the probability of the occurrence of i th ordinal pattern denoted by π and D indicates the length of ordinal patterns.

Statistical test
Most statistical tests in this study compare metrics of different NCSs with each other or compare NCS measures between clinical groups (sleep state or brain injury status). These tests were done using Wilcoxon's rank sum test. In the case of multiple comparisons, we corrected the p-values using the Benjamini-Hochberg method (Benjamini, 1995).

Results
As described above and in the supplementary material (see Supplementary Fig. S1), four cortical activity states were found to optimize the key analytical targets: robustness and stability. We also tested the significance of the HMM states using surrogates (see Supplementary  Fig. S2). In the following, we will describe the characteristics and clinical correlates of these four cortical states.

Spectral features of the newborn cortical states
The spectral analysis of the state time courses showed that all NCSs have narrow-band spectral profiles (see Supplementary Fig. S3), enabling us to extract spectral features of NCSs by means of the statewise cross spectral estimation. The spectral features in newborn cortical states (NCS) over the frequency range 0.4-22 Hz were found to organize into complementary pairs of states 1 vs. 2 and states 3 vs. 4, respectively (Fig. 2a). The relative power (compared to all combined other states) was high in NCS1 and low in NCS2 in the frontotemporal regions (p < 0.001) (Fig. 2a). NCS3 was characterized by a higher occipital and lower frontotemporal activity, while NCS4 was characterized by a reduced activity in most posterior areas with only higher power in a few frontal cortices proximal to midline (p < 0.001). The average frequency spectra showed a clear 1/f-like decline in all NCSs; however, NCSs were mutually different with the highest power across the whole frequency range in NCS1, in turn NCS2 had the lowest power (Fig. 2b). A measure of bivariate cortico-cortical interactions, the average coherence, was comparable in shape between all states, but NCS1 showed the highest coherence strength (Fig. 2c). A scatter plot combining coherence and power showed clearly discrete clustering of NCS1 and NCS2, whereas NCS3 and NCS4 were partly overlapping (Fig. 2d). The more detailed spectral topographies of the NCSs are shown in the supplementary material (see Supplementary Fig. S4).

Large-scale cortico-cortical connectivity in the newborn cortical states
The large-scale cortico-cortical connectivity by phasesynchrony showed robustly different patterns between the NCSs (Fig. 3). The NCS1 exhibited a significantly higher phase-coupling strength compared to other states throughout all cortical regions (p < 0.001), involving up to 19-44% of connections linking specific anatomical regions. The other NCSs showed lower phase-coupling strength (p < 0.001) and distinct spatial fingerprints with the most prominent differences seen in frontal, fronto-occipital, and fronto-temporal connectivity. These networks comprised all anatomical regions but with predominant of the frontal cortices in all states. The connection lengths were found to be different between states (Fig. 4): most connections in NCS1 and NCS2 are midrange (~ 4.5 cm), NCS3 is dominated by short-range connections (< 3 cm), and NCS4 presents with a somewhat bimodal distribution including both short-range (< 3 cm) and long-range connections (> 6 cm). Further detailed properties of the state-and frequency-specific networks are shown in the Supplementary Figs. S5 and S6.

Sleep effects on the newborn cortical states
The NCSs consist of brief segments that together form sequences, or temporal dynamics which can be studied as NCS time courses; in turn, they may be studied for spontaneous behavior that possibly correlates to subjects' behavior and/or vigilance states, such as sleep in the newborn infant context (André et al., 2010;Omidvarnia et al., 2015;Tokariev et al., 2016b).
We first inspected the temporal state properties, fractional occupancy and mean lifetime and found that quiet sleep (QS) and active sleep (AS) are mostly comparable within each NCS. The only differences were found in NCS1, where fractional occupancy was significantly higher in quiet sleep compared to active sleep (p < 0.001), while mean lifetime was also greater NCS in quiet sleep compared to active sleep (p < 0.001) (Fig. 5).
Next, we wanted to analyze how the sleep states affect the temporal sequences of NCSs. The transition probabilities between NCSs (no selftransitions were considered, e.g., 1→1) were computed and compared Fig. 2. Spectral topographies of the identified newborn cortical states (NCSs). The power maps (a) are in relation to the average across states and across the frequency range. All maps were thresholded in such a way that only 50% of the cortices with the strongest absolute value are shown. The averaged power (b) and phasecoupling across brain regions (c) are shown as a function of frequency. Phase-coupling versus power is shown in panel (d). Each dot represents a brain region and phase-coupling is calculated as the summation of the phase-lag index value between those regions with the rest of the regions. (L: left hemisphere, R: right hemisphere).
between AS and QS (Wilcoxon's ranksum test). Transitions 3→1 and 2↔4 was observed more during AS than QS (p < 0.001), while QS exhibited significantly higher number of transitions 2→3 and 4→1 (p < 0.001). Notably, there were no direct transitions between NCS1 and NCS2 during either sleep states (Fig. 6). Additionally, we also examined the longer NCS sequences by estimating complexity (permutation entropy), which indicated that state transitions exhibit higher complexity in QS than in AS (p < 0.05) (Fig. 6).

Effects of a neonatal brain injury on newborn cortical states
The above findings have characterized discrete NCSs and how they are linked to physiological neurobehaviour, infant's sleep. Finally, we studied whether the reported NCSs are affected by clinically occurring medical adversities, such as neonatal brain injury. To this end, we examined NCSs in a cohort of infants suffering from mild to moderate degree of hypoxic-ischemic encephalopathy (HIE) as a cause of birth asphyxia. It is the most common form of neonatal neurological adversity leading to treatment in the neonatal intensive care units with a significant risk of neurodevelopmental compromise (Douglas-Escobar and Weiss, 2015;Kurinczuk et al., 2010). The healthy control (HC) group that was analyzed in detail above was used as a reference.
Temporal characteristics of brain states showed multiple significant differences in HIE group compared to HC (Wilcoxon rank sum test): During active sleep, the HIE infants showed an increased fractional occupancy in NCS2 (p < 0.001) and a decrease in NCS3 (p < 0.001). During   Fig. 3. Cortico-cortical connectivity of the identified newborn cortical states (NCS). Network strength is shown relative to the mean strength across other states (red is over the average and blue is below the average). The connectivity rings show the strongest connections after thresholding of the full networks with Gaussian mixture model approach (see Methods section and Supplementary Fig. S7). The bottom row shows the involvement index for each anatomical brain region in the network. For example, NCS1 is supported mainly by frontotemporal regions while NCS2 is supported by the frontal region. (L: left hemisphere, R: right hemisphere, F: frontal, C: central, T: temporal, O: occipital).

Fig. 4.
Connection-length differences between the identified newborn cortical states (NCS). The probability distribution of the Euclidean distance between the parcel centroids in the identified NCS shows that they are different in terms of the connection length. quiet sleep, the HIE infants showed a longer mean lifetime in all NCSs (p < 0.001 for all); likewise, the mean lifetime of NCS1 and NCS2 was longer during active sleep in the HIE infants compared to the healthy group (p < 0.001) (Fig. 7). There were also significant effects by HIE on the temporal state dynamics (Fig. 8), in particular in the complexity of state transitions in quiet sleep; the longer sequencies of state transitions patterns showed a lower permutation entropy in the HIE infants compared to the healthy infants (p < 0.01). Transition probabilities showed no clear effects from HIE, with only the transition probability 3→4 being lower in the HIE infants during quiet sleep (p < 0.001).

Discussion
We show here that the newborn cortical activity organizes spontaneously into discrete newborn cortical states (NCS) that exhibit their respective characteristics in terms of spectral content, large-scale network connectivity, and temporal sequences. Moreover, we show that NCSs link differently to vigilance states as well as to a common neurological adversity, hypoxic-ischemic encephalopathy due to birth asphyxia. The findings are fully in line with the prior studies in adults showing that cortical activity can be characterized as a sequence of discrete network states (Baker et al., 2014;Stevner et al., 2019;Vidaurre et al., 2018bVidaurre et al., , 2018aVidaurre et al., , 2017. Here, we extend earlier knowledge by showing that discrete cortical states are robustly present already at birth and that they are nested within qualitatively different forms of neonatal cortical activity associated with distinct vigilance states (Vanhatalo and Kaila, 2006;Wallois et al., 2021); we also show for the first time that such cortical states may correlate to both physiological and pathophysiological conditions in the newborn brain.
The complex and rapid global dynamics in human cortical function call for methods that can reliably capture such changes in a way that also links to neurobehavioral read-outs. It has been shown recently that cortical activity can be described as a sequence of microstates, i.e. segments with distinct scalp potential distributions (Khanna et al., 2015;Michel and Koenig, 2018) where cortical connections may remain somewhat stable (Hutchison et al., 2013;Preti et al., 2017). Recently, many methods have been proposed to capture these dynamics ranging from EEG microstates (Khanna et al., 2015;Michel and Koenig, 2018) to hidden Markov modeling (Quinn et al., 2018;Vidaurre et al., 2016Vidaurre et al., , 2018b. While microstate analysis only considers scalp potential fields, the presently used HMM-based approach identifies brain states based on their changes in the cortical connections. A large number of studies have reported cortical network dynamics in relation to cognition, pathology, behavior (Fosque et al., 2022;Hutchison et al., 2013;Murphy et al., 2020;Preti et al., 2017;Sanchez-Alonso et al., 2021;Siebenhühner et al., 2016), and age (Coquelet et al., 2020;França et al., 2022;Ma et al., 2020;Wen et al., 2020). Our study extends these data by assessing dynamics in spontaneously occurring cortical network states in the neonatal brain.
Our findings confirm prior studies on sleep-related differences in the static networks (Tokariev et al., 2019a(Tokariev et al., , 2019b(Tokariev et al., , 2016; however, we identified global networks that are present in both vigilance states, yet showing robust sleep state differences akin to the prior observations in the adult brain (Stevner et al., 2019). The deeper sleep, QS, was found to exhibit pronounced amounts of NCS1 that is characterized by frontotemporally prominent power and a solid brain-wide neuronal coupling. This is directly compatible with previous work on static networks reporting that QS links to a broadly increased cortical connectivity (Tokariev et al., 2019a) and power (Tokariev et al., 2016b). The longer LT in QS implies a higher stability, being in line with the observations from scalp-potential-based microstates (Khazaei et al., 2021). These observations together are compatible with an idea that the rapidly changing cortical dynamics related to information processing is generally decreased in the deeper sleep (see also Tononi and Massimini, 2008).
We also observed the effects of perinatal brain damage on the cortical network states. AS periods exhibit pronounced amounts of NCS2, which is characterized by a lower brain-wide neuronal coupling. That observation is entirely in line with the previous study showing HIErelated reduction in cortical coupling (McLaren et al., 2019). A general finding across sleep states in the HIE cohort was that all NCSs exhibited longer lifetimes, suggesting higher stability or conversely, lower flexibility compared to the healthy infants. A reduced flexibility or dynamics has been suggested as a general property in many brain pathologies (Hallett et al., 2020;Kobeleva et al., 2022;Polverino et al., 2022;Van Schependom et al., 2019).
This study also has some methodological limitations that should be considered. First, like other clustering-based approaches, HMM analysis does not yield a unique number of states. Fewer states may lead to a toosimplified interpretation of the data, whereas more states may divide one mental process into many steps and contain sporadic states (Borst and Anderson, 2015;Quinn et al., 2018). Second, HMM assumes states are mutually exclusive (Vidaurre et al., 2018b). The probabilistic state assignments can explain every time point, but only one state is considered to activate at each time point. This assumption may limit expressible patterns to HMM states. EEG microstate analysis, which divides scalp maps into microstates, uses the same assumptions (Michel and Koenig, 2018). A recent microstate analysis study suggested that EEG microstates are spatially and temporally continuous rather than discrete neuronal population activations (Mishra et al., 2020). Third, the Markovian assumption states that the current cortical state is only reliant on the cortical state at the previous time point. Consequently, the dependence between state occurrences decays exponentially with time and does not reach very temporally distant data; this is known as temporally-distant dependency. There is an evidence that the brain exhibits long-term dependencies (Sikka et al., 2020). Fourth, due to low-density recording montage, we had to remove about 40% of unreliable phase-locking connections from the analysis. Lastly, although the common HIE group in this study consisted of three clinical subgroups  . The significant differences between HC (healthy control) and HIE (hypoxic-ischemic encephalopathy) groups are reported by p-value and z-statistic (Wilcoxon's ranksum test) in relevant plots. In each distribution, medians are depicted with a horizontal line and interquartile ranges with darker shading. Fig. 8. Effect of asphyxia on the transition between newborn cortical states (NCS) during AS and QS. The permutation entropy as a function of pattern length (top panel, HC (healthy control): in green, HIE (hypoxic-ischemic encephalopathy): in purple) shows that the transitions during QS are more complex for the HC group for long patterns (p < 0.01). The asterisks show significant differences between HC and HIE groups. The significant transitions (bottom panel) during QS and AS were identified by means of the Wilcoxon's ranksum test (p < 0.001). The green arrow shows that the transition probability is greater in the HC group. The dashed arrows show non-significant transitions.
(mild, grade I, and grade II) due to their small samples, reliable statistical comparison to test for the effects of HIE severity on the dynamics of the cortical networks was not possible. For the clinical evaluation of the proposed pipeline, it can be further tested for its: i) diagnostic value by discriminating HIE severity, and ii) prognostic value by linking NCSs dynamics to neurodevelopment.
In conclusion, we show that the newborn cortical network activity organizes spontaneously into rapidly alternating states (NCS) with their respective, stable spectral and connectivity patterns. These NCSs are found in both cardinal sleep states, but with different dynamics, indicating their likely physiological significance in supporting sleep functions. A structural insult, such as perinatal brain injury, was shown to affect the NCS dynamics. The proposed framework can be fully automatized by adding a reliable artifact removal method (Kumaravel et al., 2022;O'Sullivan et al., 2023;Webb et al., 2021) and it holds promise for offering an objective, global-level measure of early brain function for benchmarking neurodevelopmental or clinical research.

Data and code availability
The datasets presented in this article are not publicly accessible, but connectivity matrices for all subjects are available upon request to Anton Tokariev, the article's senior author. anton.tokariev@helsinki.fi.
All the computations related to TDE-HMM were done using HMM-MAR toolbox (https://github.com/OHBA-analysis/HMM-MAR) implemented in MATLAB. This toolbox was cited in the main text properly.

Declaration of Competing Interest
No conflict of interest exists. All authors approve the submission.

Data availability
The authors do not have permission to share data.