The intimate relationship between coalescent generators in very premature human newborn brains: Quantifying the coupling of nested endogenous oscillations.

Temporal theta slow-wave activity (TTA-SW) in premature infants is a specific neurobiomarker of the early neurodevelopment of perisylvian networks observed as early as 24 weeks of gestational age (wGA). It is present at the turning point between non-sensory driven spontaneous networks and cortical network functioning. Despite its clinical importance, the underlying mechanisms responsible for this spontaneous nested activity and its functional role have not yet been determined. The coupling between neural oscillations at different timescales is a key feature of ongoing neural activity, the characteristics of which are determined by the network structure and dynamics. The underlying mechanisms of cross-frequency coupling (CFC) are associated with several putative functions in adults. In order to show that this generic mechanism is already in place early in the course of development, we analyzed electroencephalography recordings from sleeping preterm newborns (24-27 wGA). Employing cross-frequency phase-amplitude coupling analyses, we found that TTAs were orchestrated by the SWs defined by a precise temporal relationship. Notably, TTAs were synchronized to the SW trough, and were suppressed during the SW peak. Spontaneous endogenous TTA-SWs constitute one of the very early signatures of the developing temporal neural networks with key functions, such as language and communication. The presence of a fine-tuned relationship between the slow activity and the TTA in premature neonates emphasizes the complexity and relative maturity of the intimate mechanisms that shape the CFC, the disruption of which can have severe neurodevelopmental consequences.

INTRODUCTION

The neural circuits responsible for processing sensory information are established early in development, through a combination of genetically determined guidance events and the activity‐dependent refinement of synaptic connections (Huberman, Feller, & Chapman, 2008; O'Leary, Chou, & Sahara, 2007; Tritsch et al., 2010). It has been suggested that activity‐dependent refinement relies on endogenous non‐sensory (Babola et al., 2018; Clause et al., 2014; Siegel, Heimel, Peters, & Lohmann, 2012; Tritsch, Yi, Gale, Glowatzki, & Bergles, 2007; Xu et al., 2011) and/or sensory‐driven oscillatory generators (M. T. Colonnese et al., 2010; M. T. Colonnese & Khazipov, 2010; Milh et al., 2006; Minlebaev, Ben‐Ari, & Khazipov, 2007; Wess, Isaiah, Watkins, & Kanold, 2017). These spontaneous endogenous activities are involved in neuronal guidance, synaptogenesis, and apoptosis. They also contribute to the functional optimization of neural networks (Ayoub & Kostovic, 2009; Kostović, Judaš, & Sedmak, 2011; Moore, Zhou, Jakovcevski, Zecevic, & Antic, 2011; Ulfig, Neudörfer, & Bohl, 2000) preparing the latter's processing of complex sensory information in a unique, learning interaction with the environment (Mahmoudzadeh et al., 2013; Mahmoudzadeh, Wallois, Kongolo, Goudjil, & Dehaene‐Lambertz, 2017). The functionality of these generators is related to specific temporal, spectral, and spatial characteristics, and a fine‐tuning between various oscillatory activities—the disruption of which might be dramatic for a developing system. In the present work, we examined the intimate temporal relationships between the activities of the very first generators (located in the perisylvian area) at the onset of cortical circuitry in very preterm infants.

In premature newborns, the characteristic signatures of some of these generators can be recorded using surface electroencephalography (EEG) (André et al., 2010; Dereymaeker et al., 2017; Lamblin et al., 1999; Vecchierini, d'Allest, & Verpillat, 2003). The generators display hierarchically nested oscillatory activities, consisting of a slow wave (SW) within which faster activities are nested (Kaminska et al., 2017; Routier et al., 2017; Vanhatalo et al., 2005). The earliest currently known EEG neurobiomarker (temporal theta slow‐wave activity [TTA‐SW]) has this characteristic feature. TTA‐SW can be recorded as early as 24 weeks of gestational age (wGA), with a peak in activity at around 27–28 wGA (André et al., 2010; Lamblin et al., 1999). Before 28 wGA, the main structural event is the accumulation of thalamic afferents in the superficial part of the subplate (Dubois, Kostovic, & Judas, 2015; Kostović & Jovanov‐Milošević, 2006), allowing synaptic interaction between thalamic afferent fibers and subplate neurons (Kostović & Judaš, 2010; Murata & Colonnese, 2019; Rakic, 1988). This constitutes a transitional period during which thalamic afferents “wait” in the subplate before being relocated to layer IV of the cortex; the afferents subsequently provide the developing cortex with direct access to the outside world (Dubois et al., 2015). From 28 wGA onward, this complex transient network (involving the subplate and the cortical plate) is able to process complex sensory information (Mahmoudzadeh et al., 2013; Mahmoudzadeh et al., 2017). This assumes that the previous steps (subtended by genetic fingerprints and/or spontaneous electrical activities) have functionally and structurally prepared the cortical auditory network in the perisylvian area, independently of any exogenous sensory information (Adebimpe, Routier, & Wallois, 2019; Mahmoudzadeh et al., 2013; Routier et al., 2017). The TTA‐SW provides an opportunity to investigate the perisylvian area during this transitional period; that is, the initial functional wiring is being setup. TTA‐SWs are delicately characterized by the coalescence (from the Latin coalescere, meaning “to unite and grow with”) of two activities with different frequencies: an SW activity (0.5–2.5 Hz) within which a theta activity (TA) (4–7.5 Hz) is nested (André et al., 2010; Vecchierini et al., 2003). What, then is the functional role of TTA‐SWs? Are they sensory‐driven neural activities? What is the underlying mechanism? TTA‐SWs are located in the area of the planum temporal (Routier et al., 2017). The occurrence of TTA‐SW increases the connectivity over the perisylvian network (and notably between frontal, parietal, and temporal structures), and might be involved in the wiring of the auditory network (Adebimpe et al., 2019). Before 32 wGA, the probability of occurrence of TTA‐SW is not modified by auditory stimuli—suggesting that this activity is not sensory‐driven (Routier et al., 2017). A source localization analysis performed with high‐resolution (64‐channel) EEG suggested that TTA‐SWs might be located in the subplate (Routier et al., 2017). Since TTA‐SWs might be involved in the initial network wiring (Adebimpe et al., 2019; Routier et al., 2017), disorganization of this “well‐oiled horology” of interaction between the constituting activities is likely to impact neurodevelopmental skills at different steps—notably those relying on and necessitating fine functional tuning of the perisylvian networks. We hypothesized that the nesting (the cross‐frequency coupling [CFC], in other words) between the SW and the TTA is of functional importance and thus has profound implications for the developing neural circuits. To address the functional status of this early non‐sensory‐driven neurobiomarker, one has to define its temporal characteristics, that is, the coupling between the different nested activities that built up the TTA‐SW. In line with our hypothesis, recent studies in rodents have demonstrated the presence of bursts of nested oscillatory activities in the neonatal primary sensory cortices at later stages of development, where they might facilitate the refinement of cortical maps (Clause et al., 2014; Siegel et al., 2012; Tritsch et al., 2007; Xu et al., 2011). It is highly likely that over the course of development, the underlying mechanisms and function of these nested oscillations determine not only where they occur but also the temporal pattern of CFC. This observation holds in adults (Chacko et al., 2018; Helfrich, Mander, Jagust, Knight, & Walker, 2018; Mak‐McCully et al., 2017; Mölle, Bergmann, Marshall, & Born, 2011; Voytek et al., 2010), where intrinsic CFC constitutes a key feature of ongoing neural activity (Engel, Gerloff, Hilgetag, & Nolte, 2013; Hyafil, Giraud, Fontolan, & Gutkin, 2015), which in turn supposedly enhances combinatorial opportunities for encoding (Fell & Axmacher, 2011; Rasch & Born, 2013) and facilitates synaptic plasticity (Bergmann & Born, 2018; Buzsáki & Draguhn, 2004; Salimpour & Anderson, 2019). In the present work, we show that this generic mechanism is already in place early in the course of development, that is, before the onset of better‐known development‐related oscillatory activities such as spindle bursts (M. T. Colonnese & Phillips, 2018; Khazipov & Milh, 2018), delta brush (Kaminska et al., 2017) and neonatal prefrontal oscillations (Brockmann, Pöschel, Cichon, & Hanganu‐Opatz, 2011).

MATERIALS AND METHODS

Subjects

We analyzed 11‐electrode EEG recordings obtained during the routine clinical monitoring of 25 healthy premature newborns (Table S1; gestational age: 24 to 27 wGA). The EEGs were recorded in incubators at Amiens University Hospital's neonatal intensive care unit (Amiens, France). All infants had appropriate birth weight, size, and head circumference for their term age, an APGAR score greater than 6 at 5 min, and normal auditory and clinical neurological assessments. None were considered to be at risk of brain damage. In particular, the results of a neurological examination at the time of the recordings had to correspond to the corrected gestational age, with no history of abnormal movements. The gestational age (estimated from the date of the mother's last period and from ultrasound measurements during pregnancy) was in line with the degree of brain maturation (evaluated on the EEG). The brain imaging results (particularly transfontanellar ultrasound and standard EEG) had to be normal. In addition, the EEG evaluation in the follow up at the time they left the NICU, according to the EEG monitory guideline of the French Society of Clinical Neuroscience for very young premature, had to be normal. The main exclusion criteria were a corrected age other than 24–27 wGA, the presence of neurologic damage (such as intraventricular hemorrhage and parenchymal malformation), and abnormal EEG results.

EEG recordings

The routine EEG was recorded with Deltamed® amplifiers and 11 Ag/AgCl surface electrodes placed according to the International 10–20 system. An additional electrode at position AFz was used as a common earth. Recordings were acquired at a sampling rate of 512 Hz, a 0.3–100 Hz bandpass filter, and an additional notch at 50 Hz. Electrode impedance was kept below 5 kΩ. The electrocardiogram, deltoid electromyogram, and oxygen saturation level were also monitored. The EEG session was video‐recorded.

EEG signal pre‐processing

The EEG recordings were analyzed with MATLAB® software (The MathWorks, Inc., Natick, MA), using the FieldTrip (Oostenveld, Fries, Maris, & Schoffelen, 2011), EEGLAB (Delorme & Makeig, 2004) and CircStat (Berens, 2009) toolboxes and custom MATLAB functions and codes. The EEG signal was pre‐processed with an automated artifact rejection algorithm. For each participant and each electrode location, the sample was z‐scored using the mean ± SD of the absolute amplitude. A sample was recognized as an artifact if the z‐score exceeded a value of 5. The 3 s segments preceding and following the artifact were marked and discarded from further analysis.

A characteristic EEG feature at this stage of development is the presence of spontaneous bursts with long quiescent long inter‐burst intervals (IBIs) (André et al., 2010). To characterize the SW‐to‐TA CFC, we detected SW and TA events but (unless stated otherwise) limited our analyses to the events that occurred during the IBI. In order to detect the bursts, each time point was z‐scored (using the mean ± SD absolute amplitude for the participant and the electrode) after applying a 0.5–20 Hz two‐pass finite impulse response bandpass filter (order: 3 cycles of the low frequency cut‐off) to the pre‐processed data. Samples were marked as bursts if the z‐score exceeded a value of 3. The 3 s preceding the burst and the 6 s following the burst were marked and excluded from the following processing steps.

Detection of SW and TTA events

Slow waves and TTAs were detected separately for each participant and each channel, using the algorithms described in supplementary information (see also the detection algorithm's flowchart in Figure S1). The detection techniques were adopted from those previously used to detect sleep oscillations (namely SWs and spindles; Gonzalez et al., 2018; Mak‐McCully et al., 2017; Staresina et al., 2015), and then adjusted to match the characteristics of the SWs and TTAs observed here. The manual detection of events among randomly sampled raw data by an experienced clinical neurophysiologist (FW) agreed well with the automated TTA‐SW detection.

Time–frequency analysis

We performed TFR to test the nesting hypothesis of SW and TTA through looking for power modulation during these events. The TFRs were computed per event (using the FieldTrip toolbox) for frequencies between 0.5 and 20 Hz in steps of 0.05 Hz, using a sliding Hanning tapered window with a frequency‐dependent length that comprised a full number of cycles (always at least two). For group‐level statistics, two‐tailed paired‐samples t tests were used to check for significant peri‐event power changes relative to the pre‐event baseline (−5 to −3 s). To correct for multiple comparisons (SW‐ and TTA‐triggered: −3 to +3 s × 0.5–20 Hz), the cluster‐based permutation procedure implemented in the FieldTrip toolbox was applied (Maris & Oostenveld, 2007). The initial threshold for cluster definition was set to p < .01. Lastly, the final significance threshold for summed t values within clusters was set to p < .01.

Phase–amplitude coupling analysis

To evaluate the PAC between the phase of the low‐frequency oscillation (the modulating SW) and the amplitude of the high‐frequency oscillation (the modulated TTA), the raw EEG signal was filtered separately in the frequency range corresponding to the SW (0.5–2.5 Hz; two‐pass finite impulse response bandpass filter; order: 3 cycles of the low cut‐off frequency) and the TTA (4–7.5 Hz; two‐pass finite impulse response bandpass filter; order: 3 cycles of the low cut‐off frequency). Phase values were calculated for all samples of each extracted SW, and the corresponding TA amplitude variations were calculated using a Hilbert transform. Theta activity power‐time series were created using TFR bins, averaged across the respective frequencies (4–7.5 Hz), up‐sampled to the sampling frequency of 256 Hz, and filtered in the SW's frequency range. The synchronization index for the two phase value time series was then computed for each event epoch and channel. The synchronization index is a complex number, the radius of which indicates the strength of locking between SW and TA power fluctuations, and the angle of which represents the “preferred phase” of synchronization (Gonzalez et al., 2018; Staresina et al., 2015)—in other words, the SW's phase at which the TA's power is greatest over time. The SI is calculated aswhere N is the number of samples, is the phase value of the fluctuations in the TA power time series at sample t, and is the phase value of SW time series at sample t (Cohen, 2008). The synchronization index was calculated for the 2 s interval centered on the SW trough. Prior to the statistical analysis, the non‐uniformity of circular data without a specified mean direction at the single‐participant level was evaluated using the Rayleigh‐test (Berens, 2009). In order to evaluate the preferred PAC phase at the group level, the V‐test was used to check for non‐uniformity of circular data with a specified mean direction (Berens, 2009). In the V‐test, hypothesis H1 indicates that the population is clustered around a known mean direction (180°, in the current case), whereas H0 states that the population is either uniformly distributed or is clustered elsewhere (but not in the pre‐specified direction).

Event histogram analysis

To precisely investigate the temporal relationship between the TTA and the SW, event correlation histograms were built. Three types of histograms were calculated for the TTA, using 10 s windows with 5 s offsets and a bin size of 100 ms. The histograms were referenced to the SW's trough. Furthermore, (a) the largest trough, (b) the first peak, and (c) the last peak of each detected event were taken into consideration. The count in each bin was divided by the number of TTA events. The first event correlation histogram calculated TTA beginning with reference to the SW trough. The second event correlation histogram calculated TTA largest trough event with reference to the SW trough. The last event correlation histogram calculated the TTA end with reference to the SW trough. For statistical analysis and to check whether the event correlation histograms indicated significant temporal relationships, a randomization procedure was applied. The latter was similar to that used previously to check the temporal relationship between spindles and slow oscillation in sleep data (Mölle et al., 2011). We randomized the time points of all detected SWs and TTAs, and used them as reference points in the event correlation histograms. Next, we recalculated the event correlation histograms for the randomized data. This procedure was performed 500 times separately for each subject. After the recalculation of individual event correlation histograms, the results were collapsed across subjects and used for statistical comparisons. For group‐level statistics, two‐tailed paired‐samples t tests were used to compare each bin with the corresponding bin in the random condition. To correct for multiple comparisons (−2 to +2 s), the cluster‐based permutation procedure implemented in the FieldTrip toolbox was applied. The initial threshold for cluster definition was set to p < .05. Lastly, the final threshold for significance of the summed t value within clusters was set to p < .05.

Comodulogram analysis

We used the method developed by (Tort, Komorowski, Eichenbaum, & Kopell, 2010) to create a comodulogram and further evaluate PAC over a broad frequency range. For a given frequency pair, the raw EEG signal was filtered into the two corresponding frequencies (two‐pass finite impulse response bandpass filter, order: 3 cycles of the low frequency cut‐off), with lower and higher frequencies ranging from 1.5–10 Hz (0.5 Hz increments, with a 1 Hz filter bandwidth) and 3–20 Hz (1 Hz increments, with a 2 Hz filter bandwidth), respectively. A Hilbert transform was used to extract the time series for the lower (modulating) frequency phase and the higher (modulated) frequency amplitude. The lower‐frequency phase values were then binned into eighteen 20° bins, and the mean amplitude of the higher‐frequency corresponding to each bin was computed and normalized against the sum of all bins. Tort et al. (Tort et al., 2010) defined the modulation index (MI) with regard to the Kullback‐Lieber distance between the amplitude distribution P and a uniform distribution U, , where the Shannon entropy H of the distribution P is . The MI ranges from 0 to 1; a value of 0 shows that the mean amplitude is uniformly distributed over the phases, and a MIof 1 shows that the mean amplitude has a Dirac‐like distribution.

For the comodulogram analysis, epochs were defined as 12 s of continuous artifact‐free EEG data centered on the SW's trough, plus an additional 2 s of artifact‐free data for filter padding. To accommodate inter‐subject differences in the number of SW event during the IBI (mean number: 45.15 ± 4.93 after artifact rejection), we calculated the MI for 30 randomly selected concatenated epochs (thus ensuring that the index was calculated from the same amount of data). Data from six subjects were not included in this specific analysis because the number of SW events during the IBIs was below 30. For statistical analysis across subjects, surrogate PAC data were created by calculating the MI after permutation of the epoch numbers corresponding to amplitude and phase time series. This procedure was repeated 500 times for each frequency step. Lastly, the cluster‐based permutation procedure implemented in the FieldTrip toolbox was used to compare the original comodulogram with that corresponding to the surrogate data at the group level. The initial threshold for cluster definition was set to p < .01, and the final threshold for significance of the summed t value within clusters was set to p < .01.

Effect of auditory stimulation on the occurrence of TTA‐SW

To determine whether the TTA‐SW is sensory‐driven, we studied the impact of click stimuli on the occurrence of TTA‐SW in two additional groups of premature newborns (see Table S2 in supplementary information). In the first group, nine premature infants were tested by applying the protocol used previously to evidence the sensory‐driven nature of delta brushes over the temporal areas before 35 wGA (Chipaux et al., 2013). In the second group, six premature infants were tested using stimuli that mimic the fetal environment (i.e., where the womb acts as a low pass filter; Granier‐Deferre, Bassereau, Ribeiro, Jacquet, & DeCasper, 2011). Two loudspeakers were placed inside the incubator (Logitec VR, 20–20 kHz), 10 cm away from the infant's head. The click stimuli consisted of either a 4 ms square signal at 5000 Hz and70 dB (group 1) or a 400 ms square signal at 500 Hz and 70 dB (group 2). To avoid habituation a pseudo‐random sequence was applied every 20–25 s throughout the sleep period. The mean ± SD total number of clicks delivered per child was 115.44 ± 4.07 in group 1 and 120.33 ± 4.49 for group 2. Online synchronization between the stimulus delivery system and the EEG recording system was performed via a built‐in MATLAB procedure, which provided a marker for each click on the EEG acquisition device. The click stimuli period started after 15 min of background activity had been recorded.

In order to determine whether or not the TTA‐SW could be triggered by auditory stimuli, TTA‐SW events (including bursts and IBIs) were detected using the automatic approach described above. The TFRs were analyzed (as in the event‐locked analysis, to assess whether stimulation trials evoked an oscillatory response with reference to baseline (−3 to −5 s, with respect to the stimulus onset). Furthermore, the event histograms of the TTA‐SW events were computed using a 6 s window and a 75 ms bin, with respect to the stimulus onset. Each TTA‐SW event was defined with respect to the SW's trough in the filtered data (0.5–2.5 Hz). To test the statistical significance of (a) the oscillatory power changes in the TFR and (b) the TTA‐SW's probability of occurrence in the event‐related histograms, the statistical approaches were similar to those described above for the TFR and event histogram analyses.

RESULTS

To demonstrate the specific, highly tuned temporal relationship between SW and TTA, we characterized the phase–amplitude coupling (PAC, that is, the moment at which theta rhythms appear in the course of the SW) (Figure 1), on 54.91 ± 14.25 min EEG recordings of 25 premature neonates (24–27 wGA). To determine whether the SW phase systematically modulated the amplitude of the theta oscillations, we developed four complementary analytical approaches after selectively identifying the SW and TTA events, with a specific automatic detection algorithm. Firstly, to determine whether SW and TTA events were nested, we looked for power modulations over the course of the respective events. To this end, we performed an event‐locked analysis, and aligned the time–frequency representations (TFRs) of the peri‐event epochs with the center of SW and TTA events. Secondly, to investigate the temporal relationship between TTA and the SW in detail, we created event correlation histograms that were referenced to the negative half‐wave peaks of the SW. Thirdly, to extract the preferred phase of coupling, we analyzed the PAC. Fourthly, to further evaluate the selectivity of SW‐TTA coupling, we performed a comodulogram analysis and evaluated the PAC across a broader range of frequency pairs.

FIGURE 1

EEG recorded in a premature infant (age 25 wGA; band‐pass filter: 0.5–20 Hz). A TTA‐SW event is clearly seen over FP1‐T3 and T3‐O1

Furthermore, in order to verify that TTA‐SWs were not triggered by exogenous auditory inputs, we studied the responses of an additional 15 newborns (24–27 wGA) to click stimuli via two analytical approaches. Firstly, to determine whether there were significant power differences at various frequencies before vs. after stimulation, we analyzed the TFRs. Secondly, to test whether auditory stimulation would modulate the probability of the occurrence of TTA‐SW (relative to a random condition), the peri‐stimulus histogram was analyzed.

The results are presented below for the bipolar FP1‐T3 channel in the left hemisphere. To demonstrate the consistency of the results across the two hemispheres, all the proceeding results are presented over FP2‐T4 in the right hemisphere in the Supplementary Information.

Our automated event detection algorithm yielded an average of 98.22 ± 9.10 SWs and 81.43 ± 6.98 TTAs per infant over the left hemisphere. The mean event densities were 1.85 ± 0.15 per min for SWs and 1.52 ± 0.12 per min for TTAs. On average, 44.56 ± 3.78% of all TTAs occurred within 3 s of an SW's negative peak.

In preterm newborns, endogenous mechanisms induce cortical bursts of spontaneous high‐amplitude activities separated by long quiescent periods (i.e., IBIs) (André et al., 2010; Lamblin et al., 1999). In order to characterize the SW‐TTA CFC, the analyses presented below were (unless stated otherwise) limited to events that occurred during IBIs.

Nesting of SWs and TTAs, and the phase of coupling

The TTA power (4–7.5 Hz) is modulated during the SW event. First, TFR time‐locked to the SW trough (Figure 2a) demonstrates that TTA power increased during the descending phase of the SW, peaked during the SW trough (slightly preceding the SW trough), (zoomed in from ‐2 to 2 s and on 3–10 Hz to highlight the nesting of TTA power in the SW trough), decreased during the subsequent positive slope, and showed a lasting suppression during the SW peak. After thresholding, a significant cluster emerged around the SW trough (p < .01, corrected), relative to the pre‐event baseline interval of −5 to −3 s. The TFR locked to the TTA's deepest trough (Figure 2b) was associated with greater SW power (0.5–2.5 Hz) during the TTA event. A significant cluster emerged at the TTA event (p < .01, corrected); the SW frequency range was statistically significant from −1 to 1 s, relative to the pre‐event baseline. Furthermore, two clusters emerged 1–2 s away from the event centers. A PAC analysis (assessing the preferred modulation phase of TTA power by the SW) highlights the high temporal precision of the TTA power nesting within SW cycles (slightly before the SW trough) for all events (Figure 2c, Video S1) and for each participant (Figure 2d, V = 26.56 p < .001, average preferred phase = 175.21°). The Rayleigh test for non‐uniformity was significant (p < .05) for all subjects.

FIGURE 2

Event‐locked analysis of SW‐TTA relationship. (a) Grand average of SW trough‐locked TFR (zero at SW trough). The bottom panel shows the TFR power averaged over the theta range (4–7.5 Hz). (b) Average of TTA locked TFR (zero at the largest trough of TTA). In panels (a) and (b), the statistically significant change from the pre‐event baseline (−5 to −3 s, p < .01, cluster‐corrected) is indicated with a white contour. (c) Normalized histogram of the preferred modulation phase of SW to TTA pooled over all newborns. The average phase is indicated by a red line. (d) The unit circle of preferred phases of SW to TTA modulation in each newborn. The Rayleigh test for non‐uniformity was significant (p < .05) for all subjects. The preferred phases of SW to TTA modulation clustered significantly around the trough of SW across neonates, as shown by the red line (V = 26.56, p < .001, average preferred phase = 175.21°). (e) Grand average 0.5–20 Hz EEG trace across neonates, aligned to the largest TTA trough (trough, time 0). Intervals whose amplitude differed significantly from baseline (−5 to −3 s) are shaded (p < .01, corrected). The green arrow marks the location of the “blup.” (f) Grand average, 0.5–20 Hz EEG time‐locked to the first filtered positive theta peak. Blue and black lines indicate the top and bottom quartiles, respectively. The red arrow shows that specifically in the top quartile the first peak of the TTA is at the onset of the SW

To examine the general structure of the detected TTAs, the filtered EEG (0.5–20 Hz) signal was averaged time‐locked to all the detected TTAs. The superimposing of all the detected TTAs (locked to the deepest trough) revealed that the TTA event is nested around the trough of SW, which in turn displays an alternating pattern (as shown in Figure 2e). The amplitude of the signal (found to be statistically significant, relative to baseline) in the shaded regions highlighted two positive clusters and two negative clusters (p < .01, corrected). Next, we looked for the presence of oscillatory activities in the averaged filtered EEG signal that was time‐locked to the first peak of the filtered TTA. To this end, we sorted the events according to the TTA trough's amplitude and then averaged the signal across all events within a quartile that was time‐locked to the first peak in the band‐pass‐filtered TTA. This analysis revealed clear oscillations within the theta range, and showed two to three peaks for the top quartile; the peaks were less pronounced in the lowest quartile (Figure 2f). The data also showed that the first detected peak is located on the SW's descending slope. Inspection of Figure 2e,f shows a significant increase from the baseline preceding the onset of the descending slope. We chose to call this increase from baseline a “blup,” by analogy with the sound that boiling water makes. The TTA is nested within a SW; the first peak is located in the descending slope of the SW, and its deepest trough slightly precedes the SW's trough. The amplitude of the SW's subsequent rebound peak was also significant. The peak was followed by a small but statistically significant negative deflection.

On average, 23.12 ± 2.23% of the TTAs were located at between −200 and 200 ms, relative to the SW trough (t = 0). However, the result of the event histogram analysis (i.e., the proportion of TTAs associated with SWs) depends on the event detection criteria for SWs and TTAs. If (for example) the TTA detection threshold was set to 90% (and not 95%) of the root mean square value, the proportion of TTAs associated with SWs (based on our event detection criteria) increased from 27.34 ± 1.9% to 32.87 ± 2.98% for the left hemisphere. Event correlation histograms (built with randomized data across subjects) were used to check the reliability of these temporal relationships. Statistical analysis of these histograms indicated a consistent timing of the TTA during the SW cycles (Figure 3, with t = 0 at the SW trough). In agreement with the TFR and PAC analyses, the TTA occurred at the onset of the SW descending slope (p < .05, corrected, from −700 to ‐300 ms, compared to the random condition) and then showed a pronounced increase around the SW trough with a distinct maxima reached between −300 and +100 ms (p < .05, corrected, compared to the random condition). The histogram's mean skewness was ‐0.23 ± 0.14. We observed that the TTA events ended during the ascending slope of the SW and before the SW's rebound peak (p < .05, corrected, from +100 to +550 ms, compared to the random condition).

FIGURE 3

Event‐related histograms (normalized for the total number of counts) showing the timing of the beginning of TTA events in relation with the trough of the SW (t = 0) (a), the timing of the largest trough of the TTA events (b), and the timing of the end of TTA events (c). The data in the bins are quoted as the mean ± SE. The red line indicates the mean ± SE for the randomized condition. The results of subject‐level statistical analysis are shaded (p < .05, corrected)

The selectivity of SW‐TTA coupling

To further evaluate the selectivity of SW‐TTA coupling, we investigated PAC across a broader frequency range by applying the comodulogram method (Tort et al., 2010). The MI reflects the degree to which the amplitude of the higher (modulated) frequency varies as a function of the phase of the lower (modulating) frequency. In order to assess the statistical significance of the observed CFC, we compared the results with those generated with epoch‐shuffled surrogate data—the same data as for the original PAC analysis, and with exactly the same spectral power characteristics. This procedure highlighted a single positive cluster with an MI greater than the 90th percentile (p < .01, corrected, Figure 4). This cluster corresponded to the SW‐TTA PAC also observed in the event‐locked analysis. These PAC results not only corroborated the main findings from our event‐based analysis but also highlighted the preferred phase of the modulating frequency within which the higher frequency power is nested (see the inset in Figure 4).

FIGURE 4

Phase–amplitude coupling for TTA‐SW events. A single cluster showed a significant modulation index when comparing the original data with the epoch‐shuffled surrogate data (p < .01 in a two‐tailed, paired‐samples t test, and with a modulation index greater than the 90th percentile). The inset shows the probability distribution of average amplitudes of the modulated high‐frequency as the function of the phase of the modulating low frequency

The endogenous non‐auditory sensibility of the TTA‐SW before 28 wGA

In line with our previous results (Routier et al., 2017), click stimuli similar to those used to evoke delta‐brush activities in preterm newborns before 35 wGA (Chipaux et al., 2013; Kaminska et al., 2017) did not evoke TTA‐SW events between 24–27 wGA. TTA‐SW events for both bursts and IBI were considered in this analysis. A time‐locked analysis (relative to stimulation onset) did not reveal any significant changes in the power of the frequency bands, when compared with the pre‐stimulus baseline of −5 to −3 s (Figures S2 and S6). Furthermore, a correlation histogram for peri‐stimulus events did not reveal a significantly greater probability of TTA‐SW events than in the randomized condition. These findings suggest that in the premature newborns studied here, TTA‐SW is not auditory‐driven.

DISCUSSION

Our main results demonstrate, (a) a precise temporal relationship between TTAs and SWs in early stage of development, in a way that TTAs are orchestrated by SWs, and (b) that the TTA‐SW is not evoked by auditory stimulation, which emphasizes the endogenous nature of this early neurobiomarker, at least before 28 wGA. The strength of the study lies in the fact that the results were observed using complementary analyses.

TTAs were synchronized to the SW trough, and were suppressed during the SW peak. In event‐locked analyses, we identified SW and TTA events with established detection algorithms, and used the SW event centers to align the TFRs of faster neuronal oscillations and to assess the TTAs' probability of occurrence. This analysis revealed an increase in the theta‐range spectral power (4–7.5 Hz) on the descending slope of SW, with its peak preceding the trough of the SW. The event histogram analysis also highlighted the disappearance of TTAs at the ascending slope of SW events. We also used a comodulogram analysis to evaluate the PAC of an exhaustive range of frequency pairs. Together, the PAC and comodulogram analyses strongly corroborated the findings of the event‐locked analyses, and revealed the same phase‐to‐amplitude relationship (i.e., an increase in TTA power during the SW trough). Our present results provide evidence that at this early stage of development, SWs coordinate the occurrence of faster theta oscillations, and together they shape the earliest endogenous EEG biomarker of neurodevelopment (TTA‐SW) recorded from 24 wGA onward (André et al., 2010; Lamblin et al., 1999)—long before the thalamic afferents are relocated to the deep cortical layers (Dubois et al., 2015; Kostović & Judaš, 2010).

The SW component of the TTA‐SW has an alternating pattern. The TTA starts after the “blups” (the initial reversal), which might trigger the descending phase of the SW. In sleep, such initial slow component has been attributed, to the preceding slow waves, which were the object of the analysis (Neske, 2016). In the present study, we did not observe a train of SWs, single SWs arise in isolation (as observed in Figure 1 in the present article and figure 1 in Routier et al., 2017). The “blups” thus correspond to the first significant change in the EEG signal during the TTA‐SW event. The phase relationship is of interest as the maximum amplitude of the TTA slightly precedes the SW trough, suggesting that the beginning of the fading of theta activity is modulated by the SW. Lastly, the end of the TTA lies in the ascending part of the slow wave. Overall, the SW is likely to reflect the activity of generators that synchronize a train of theta oscillations confined to the descending part of the SW, providing a window of opportunity for the emergence of theta activity. Alternatively, TTA might initiate the SW ‐as has been suggested for the SW associated with sleep spindles in adults (Gonzalez et al., 2018). However, the event histograms observed in the present study do not support this hypothesis because (a) the theta activity appeared to be located in the descending slope of the SW, and (b) the beginning of a “blup” preceded the onset of the theta activity. However, a threshold effect on our detection of the onset of theta activity cannot be ruled out; these effects might complicate the detection of the event onset—particularly in surface EEG.

Phase–amplitude coupling is a potentially useful measure of coupling between neural oscillations on different timescales. The mechanisms underlying PAC have recently received a lot of attention in both experimental and theoretical studies. It has been suggested that PAC supports the encoding, storage, and retrieval of information (Bergmann & Born, 2018; Fell & Axmacher, 2011). In adults, PAC translates as precise temporal relationships between modulating and modulated frequencies in (for instance) the thalamocortical and hippocampal networks during sleep (Staresina et al., 2015), in the hippocampus during the operation of multi‐item working memory (Axmacher et al., 2010), and in the cortical networks during cognitive functions (Canolty et al., 2006; Chacko et al., 2018; Combrisson et al., 2017). It has been hypothesized that the phase of the slower oscillation generally reflects greater excitability among postsynaptic neurons, which in turn synchronizes the synaptic input (as reflected by an increase in the amplitude of the faster oscillation) (Bergmann & Born, 2018). In adults, the mechanisms underlying PAC lead to specific temporal patterns of coupling across multiple frequencies at both local and distal sites (Bergmann & Born, 2018; Engel et al., 2013; Hashemi, Dehnavi, Moghimi, & Ghorbani, 2019; Hyafil et al., 2015). This coupling can be modulated by structural changes (Helfrich et al., 2018; Salimpour & Anderson, 2019), and task‐related network dynamics (suggesting a functional role for the CFC) (Canolty & Knight, 2010; Combrisson et al., 2017; Engel et al., 2013; Haegens, Nácher, Luna, Romo, & Jensen, 2011; Tort, Komorowski, Manns, Kopell, & Eichenbaum, 2009). Ever since the 1950s, it has been accepted that fast oscillatory activities (at various frequencies <30 Hz) in newborns are superposed on an SW (e.g., delta brushes) from 28 wGA to full‐term (Dreyfus‐Brisac, 1962; Lamblin et al., 1999). More recently, slow activity transients have been described as ultra‐SWs within which faster, more complex activities are observed (Vanhatalo et al., 2005). By extrapolating the data from adults with caution, one can consider that distinct CFC patterns (as at least visually observed in different neurobiomarkers) reflect the mechanisms of the developing wirings and the functionality of the future neural networks. Therefore, considering the possible issues in premature dysfunctional networks, which might lead to neurodevelopmental disorders, detailed analysis of the temporal patterns of relationships across multiple frequencies in early neurobiomarkers of development is of both fundamental and clinical value. The precise phase–amplitude coupling we observed here for TTA‐SW events is rather similar to that reported for theta activity and spindles in the different “up” and “down” states of SWs in adult sleep (Gonzalez et al., 2018; Mak‐McCully et al., 2017; Mölle et al., 2011; Staresina et al., 2015). In fact, the phase–amplitude relationships between slow and fast activities are associated with specific functions, depending on when they occur during “up” or “down” phases. However, we do not mean to suggest that coupling in adults and coupling during TTA‐SW events in newborns necessarily share the same neural mechanism. However, we believe that the existence of this fine‐tuning from 24 wGA onward testifies to an essential functional interaction on the cellular level, the disruption of which is likely to impact the immature network's functional and structural status. The extent to which the PAC observed in premature newborns is a general feature and can be extrapolated to the SW in adult sleep must now be determined.

At this early stage of development, the source of the TTA‐SW is located in the subplate (Routier et al., 2017), which has a key role in (a) the establishment and functional maturation of thalamocortical connections through projections to the growing cortical plate (Dubois et al., 2015; Kostović & Judaš, 2010), (b) feedback projection to the thalamus (Murata & Colonnese, 2016), (c) the activity‐dependent development of cortical columns (Elias & Kriegstein, 2008; Elias, Wang, & Kriegstein, 2007; Rakic, 1988), and (d) radial neuronal migration (for a review, see Luhmann & Khazipov, 2018). Neuronal dysfunction of the subplate has been implicated in neurodevelopmental disorders because of a resulting hyper‐connectivity, which at this period of development is mainly related to the interactions between the subplate and the cortical plate (Nagode et al., 2017). Furthermore, spontaneous endogenous activities are involved in the underlying mechanisms that build up the neocortical network (for a review, see Kilb, 2012). Therefore, neurogenesis, cell death, cell migration, neuronal differentiation, and the development of functional synaptic pathways (i.e., the first visual and tonotopic maps, somatotopic arrangement, and the formation of barrels) (Minlebaev et al., 2007; Tritsch et al., 2007) require delicate interplay between various cellular populations at different levels in the neuronal network. In the perisylvian region, the precise pattern developed around TTA‐SW might be the hallmark of this mechanism in action before the thalamic afferents relocate to the cortical plate, consistent with the protomap hypothesis (Rakic, 1988), that is, prior to the functional and structural refinement of the network by complex sensory inputs arising from the external world. It is therefore likely that the fine‐tuned temporal relationship between TTAs and SWs is a marker of endogenous local interactions between the subplate, cortical plate, and thalamic afferents. Indeed, even if these pattern generators are located in the subplate, they cannot rely solely on subplate neurons—the orientation of which does not allow generation of a recordable bipolar electric field on the scalp (Luhmann & Khazipov, 2018).The latter requires synchronized dendritic activity by well‐oriented migrating pyramidal cells. Connectivity across the perisylvian network (notably between frontal, parietal, and temporal structures) increases during the occurrence of TTA‐SW—suggesting that the latter feature might also facilitate the development of the perisylvian network. It is known that from 28 wGA onward, the perisylvian network is functionally sufficiently developed to discriminate phonemes and voices (Mahmoudzadeh et al., 2013; Mahmoudzadeh et al., 2017).

The ongoing electrical activity of the brain is organized into functional resting state networks (Raichle, 2015). These networks are suggested to be present in term infants (Fransson, Åden, Blennow, & Lagercrantz, 2011) as well as preterm infants (Doria et al., 2010; van den Heuvel et al., 2015), but to be organized around hubs in sensory and motor cortex, and to mature in the course of development (He & Parikh, 2016). During the early development of neural circuits, the spontaneous activity in the primary sensory regions is determined mostly by transitory generators in the sensory periphery (e.g., spontaneous retinal and cochlear waves) (M. Colonnese & Khazipov, 2012). In this context, another hypothesis around the functional role of TTA‐SW might be its possible participation in the early development of the neural circuits corresponding to later processing of sensory information, noted in fMRI studies as resting state networks in pre‐terms with hubs in sensory and motor regions (Doria et al., 2010; van den Heuvel et al., 2015). Interestingly, spontaneous isolated nested oscillations are observed in different sensory and sensorimotor areas (over central, temporal, and parietal areas) at rather similar periods of development in preterm newborns between 24 and 28 wGA (personal data). These first oscillatory activities appear early in the course of development before the relocation of thalamic afferents in the cortical plate and are not modulated by exogenous stimulation, as we demonstrated in the case of TTA‐SW. The emergence time of these oscillations is before the emergence of the segregated and not‐yet developed higher order networks. It is therefore of interest to investigate to what extend these endogenous oscillatory activities can be considered as early neurobiomarkers of the development of multimodal sensory networks, which can later in the course of development lead to cortical sensory awareness (Padilla & Lagercrantz, 2020). This requires a multimodal analysis of the endogenous activates and a cohort which permits the follow‐up of the infants' neurodevelopment. The TTA‐SW's probability of occurrence does not increase with auditory stimulation of different frequencies—at least during this period of development (see also Routier et al., 2017). The lack of external sensory modulation of the TTA‐SW before 32 wGA does not rule out possible modulation by activities that develop spontaneously in the inner ear and travel along the sensory auditory tract toward the subplate (Tritsch et al., 2007; Tritsch et al., 2010). Spontaneous correlated activities have been also observed in other sensory systems. Activity‐dependent refinement of visual maps throughout the visual system before the onset of vision are believed to be guided by spontaneous retinal waves that convey spatiotemporal patterns (Ackman, Burbridge, & Crair, 2012), with a dynamic of interactions that changes with the progressive sensory experience as development progresses (Katz & Shatz, 1996). The retinal waves are shown to provide highly correlated patterns of spontaneous activity in the lateral geniculate nucleus and visual cortex (M. T. Colonnese & Khazipov, 2010; Weliky & Katz, 1999), which are proposed to be homologous to the slow activity EEG transients recorded in preterm human fetal occipital cortex (Tolonen, Palva, Andersson, & Vanhatalo, 2007; Vanhatalo et al., 2005). It has been shown that spontaneous activities in the cochlea triggers discrete bursts of action potentials in primary auditory neurons which in turn are required for the survival and refinement of tonotopic maps in auditory nuclei (Tritsch et al., 2007). A prudent analogy with the visual system may suggest that TTA‐SW is correlated with spontaneous cochlear activity before the onset of hearing, and independent of the auditory exogenous stimuli. These spontaneous activities thus differ from the delta‐spindles described in animals (Chipaux et al., 2013; M. T. Colonnese et al., 2010; Kaminska et al., 2017; Khazipov & Milh, 2018; Luhmann & Khazipov, 2018; Milh et al., 2006; Minlebaev et al., 2007) and the delta‐brushes observed in human newborns after 28 wGA (André et al., 2010), which are at least partly driven by exogenous sensory stimuli (Chipaux et al., 2013; Kaminska et al., 2017). It has been suggested that these delta oscillations (within which higher frequencies in the alpha and beta range are nested) are involved in amplifying and controlling sensory inputs via complex mechanisms involving the thalamus, the subplate, and the cortical plate (Murata & Colonnese, 2016). Even though the TTA‐SW and delta‐brush both display coalescence, these two EEG features differ with regard to the SW pattern (which is not a delta oscillation), the frequency of the nested oscillations (theta vs alpha/beta), and the phase of the SW within which the higher frequencies are nested (Watanabe, Hayakawa, & Okumura, 1999).Before 28 wGA, the functional synaptic contacts are rare in the cortex but the widespread arborization of the thalamocortical fibers enables synaptic interactions with the subplate neurons (Kostović & Judaš, 2010; Kostovic & Rakic, 1990; Penn & Shatz, 1999), which are the first cortical neurones being activated by sensory inputs (Wess et al., 2017). The cellular mechanisms underlying not sensory driven TTA‐SW activity (or similar localized nested oscillations at equivalent age) are not specifically described yet. However, using glass electrode recordings in the human fetal subplate neurons (17–23 wGA), it has been shown that although spontaneous unprovoked electrical activity at this period is moderately influenced by drugs that block neuronal network activity and synaptic transmission, it is more strongly influenced by drugs that interfere with the activity of connexin (Moore et al., 2014). These results suggest that before glutamatergic, GABAergic and glycinergic transmission are likely to dominate spontaneous activity (Dupont, Fourcaudot, Beekenkamp, Poulain, & Bossu, 2006), connexin hemichanels provide one of the plausible mechanisms for the generation of spontaneous unprovoked activity (Moore et al., 2014). Hemichannel mediated initiation of spontaneous activity has also been reported in the developing auditory system (Tritsch et al., 2007).

Preterm neonates are at high risk of neurodevelopmental disorders, for example; sensory motor, cognitive, and language related issues, Attention deficit hyperactivity disorder, Autism (for a review, Wallois, 2020) that might be associated with a disturbed neuronal wiring consistent with a disrupted functional organization (Padilla et al., 2017). To investigate the effect of age and the predictive value of the CFC a larger cohort including patients with abnormal outcome or specific cortical damage is necessary.

Spontaneous endogenous TTA‐SWs constitute one of the initial building blocks in the perisylvian neural networks' key functions, such as language and communication (Adebimpe et al., 2019; Routier et al., 2017). The presence of a fine‐tuned relationship between a slow activity (the SW) and an oscillator (the TTA) in very premature newborns emphasizes the complexity and relative maturity of the intimate mechanisms that shape the CFC. It is reasonable to suppose that a “grain of sand” in this well‐oiled horology would have long‐term neurodevelopmental consequences. A Better characterization of the phase–amplitude relationships for each EEG pattern recorded in premature newborns might enable researchers to define predictive neurodevelopmental criteria for the individualized management of this at‐risk population.

CONFLICT OF INTEREST

The authors declare no conflicts of interest.

Data S1: Supplementary Information

Click here for additional data file.

Video S1 Dynamic representation of TTA‐SW nested activities

Click here for additional data file.