Early development of synchrony in cortical activations in the human

Highlights • We study the early development of cortical activations synchrony index (ASI).• Cortical activations become increasingly synchronized during the last trimester.• Interhemispheric synchrony increases more than intrahemispheric synchrony.• Our EEG metric ASI can be directly translated to experimental animal studies.• ASI holds promise as an early functional biomarker of brain networks.


INTRODUCTION
Large-scale spatio-temporal correlations in neuronal activity are considered to provide the functional basis for a range of brain functions in distributed networks (Bressler and Menon, 2010;Uhlhaas et al., 2010;Palva and Palva, 2011). These correlations are readily observed in the neuronal activity, as well as in the fluctuation of cerebral blood flow (Biswal et al., 1995;Jerbi et al., 2010), and they also correlate with behavioral states (Raichle, 2010;Palva and Palva, 2011;Hutchison et al., 2013).
Little is known about the early ontogenesis of functional communication in the human neuronal networks. Recent anatomical studies have disclosed an account of the microscopic development of structural networks in the human fetus (Kostovic and Jovanov-Milosevic, 2006;Kostovic and Judas, 2010), which sets the physical frame to how the early neuronal dynamics may emerge during latter half of gestation. Some features of large-scale spatial coordination in the electrical activity of the brain have been reported in sleeping human newborns (Tokariev et al., 2012;Omidvarnia et al., 2014). It is known that two modes of brain activity (Vanhatalo and Kaila, 2006) alternate in sub-second time scales between a relative quiescence and its interruptions by spontaneous activity transients (SAT, a.k.a. burst; Vanhatalo et al., 2005). These SATs are thought to provide the endogenous driver needed for activitydependent wiring of the early brain networks, prior to onset of genuine sensory experience (Hanganu-Opatz, 2010;Kilb et al., 2011;Colonnese and Khazipov, 2012).
Early clinical studies on neonatal EEG established that the temporal co-incidence of these activity bursts between the hemispheres, commonly called ''interhemispheric synchrony" (IHS), is a good marker of normally developing EEG activity at term age. A developmental increase in IHS was found during the last trimester of pregnancy (Lombroso, 1979), however, the existing literature is based on subjective and largely qualitative EEG assessment that compromises the validity of E-mail address: ninah.koolen@esat.kuleuven.be (N. Koolen). Abbreviations: AC, Algebraic connectivity; ASI, activation synchrony index; cx-cx, cortico-cortical; EDTF, energy weighted temporal dependency function; GS, Global Synchrony; IHS, interhemispheric synchrony; iqr, interquartile range; IVH, intraventricular hemorrhage; MST, Minimum spanning tree; PMA, postmenstrual age; SAT, spontaneous activity transients.
Neuroscience 322 (2016) 298-307 detailed findings. Moreover, there are no reports on spatial differences in the development of synchrony between cortical areas, yet histological studies have clearly established distinct developmental trajectories in the growth of long-range cortico-cortical (cx-cx) pathways (Judas et al., 2005;Kostovic and Jovanov-Milosevic, 2006).
We have recently developed and validated a measure for IHS, called activation synchrony index (ASI), which statistically quantifies the temporal coincidence of SAT in the cortical activity (Ra¨sa¨nen et al., 2013;Koolen et al., 2014b). Since late 1970s, the clinical EEG review has involved visual assessment of IHS as one of the key parameters of EEG maturation and normality (Lombroso, 1979). Our benchmarking study with visual reading and other synchrony measures showed that ASI is most accurately emulating the clinically recognized phenomenon of cx-cx EEG synchrony. This has opened the possibility to study how the synchrony between cortical activations evolves during prematurity. In the present study, we aimed to characterize the developmental correlations of cx-cx synchrony during the last ten weeks of pregnancy, which is characterized by the rapid development of long-range cx-cx connections. In particular, we wanted to disclose potential spatial gradients, as well as assess whether the developmental changes are robust enough to even allow using the ASI-based cx-cx synchrony as a maturational measure.

Data acquisition
The main dataset consisted of 22 recordings in 20 infants, recorded at a postmenstrual age (PMA) of 30-44 weeks at the Neonatal Intensive Care Unit of the University Hospitals of Leuven, Belgium (Koolen, 2014a,b). In this pilot-study, we used broad inclusion criteria and have included 2 infants with intraventricular hemorrhage (IVH) grade III, however no infants had parenchymal lesions such as parenchymal infarction or periventricular leukomalacia. Most importantly, the infants were clinically stable by the time of EEG recording. More clinical details of the infant group are given in Table 1. Two infants had consecutive recordings performed for clinical reasons to assess their brain development. Their two EEG recordings were entered as independent samples, however we also computed the group results with only one recording from each infant, and we saw no meaningful differences in the results. The lower age limit was set to 30 weeks PMA, to assess the developmental window at an age when the majority of thalamocortical connections is already established, whereas an intensive growth of cortical-cortical connections exists (Kostovic and Jovanov-Milosevic, 2006;Jovanov-Milosˇevic´et al., 2009;Kostovic and Judas, 2010). One term infant was excluded because of missing trace´discontinue EEG patterns, the foundation of ASI analysis. The recording time was at minimum 4 h, and a clinical expert (A.D.) selected the most discontinuous EEG, resembling quiet sleep in older patients, for 2 Â 10 min in each recording. All EEG measurements were recorded at 250 Hz, with 8 electrodes (Fp1, Fp2, C3, C4, T3, T4, O1, O2) placed according to the 10-20 standard locations and reference electrode Cz (BRAIN RT, OSG equipment, Mechelen, Belgium). The protocol was approved by the Ethics Committee of the University Hospitals of Leuven, Belgium. Preprocessing the data involved applying a 50 and 100 Hz Notch filter and a 1-20-Hz band pass filter to capture the interesting burst information present in this frequency range.

Analysis of synchrony
We analyzed the synchrony between cortical areas by using the recently developed measure ASI, which estimates the temporal relationships between newborn cortical events (for further details, see Ra¨sa¨nen et al., 2013). Our primary aim was to study the development of ASI from prematurity to term age. In addition, we also studied spatial differences in ASI development, as well as the temporal fluctuations of ASI within each recording session in order to assess its methodological stability and potential use as a maturational measure of functional connectivity.
Technically, ASI (and mutual information in general) captures all non-linear correlations between the signals and is invariant to any coordinate transforms of the input data, which is not true for normal correlational analysis (see, e.g., Herzel and Große, 1995). We have directly compared ASI to the cross-correlation of the raw signal or the signal energy function, neither of which is suitable for capturing the temporal relations of interest in our context (Ra¨sa¨nen et al., 2013). Other somewhat comparable measures would include the correlation coefficient between amplitude envelopes (see e.g. Hipp et al., 2012;Tokariev et al., 2015b) and the computation of the slope between quantized amplitude envelopes (e.g. Omidvarnia et al., 2014Omidvarnia et al., , 2015, neither of which is designed to statistically test the presence of temporal delays, a key feature of interest in assessing cx-cx activation synchrony in the preterm EEG signal. Due to the computational robustness and the ability to emulate clinical EEG review, ASI is also a putative, practical biomarker to describe preterm maturation in the context of a realworld hospital environment.

Computation of ASI
Temporal relationships between pairs of EEG signals were computed using a previously described measure ASI (Ra¨sa¨nen et al., 2013), which was recently shown to perform well in distinguishing normality in EEG traces from term patients (Koolen et al., 2014b). ASI provides a statistical measure for the temporal coupling of two quantized EEG amplitudes with higher ASI values reflecting larger synchrony between the signals, whereas an ASI of zero means that the two signals are statistically independent. The algorithm consists of the following steps (Ra¨sa¨nen et al., 2013): 1. The signal is first down sampled to 50 Hz. Higher frequencies (which are related to the burst events) are pre-emphasized with a first order FIR high-pass filter (H(z) = 1-0.95 z À1 ).
2. Amplitude envelopes of the pre-processed signals are obtained by computing the fast Fourier transformation (FFT) using a sliding Hamming window of 2 s and a step size of 100 ms (following the optimization performed in Ra¨sa¨nen et al., 2013), and then summing up the amplitudes of the frequency bins in the relevant frequency range of 1.5-20 Hz for each 100-ms frame. 3. The signal envelopes are quantized into Q discrete amplitude levels (or ''states") by clustering a random subset of the samples with the standard k-means algorithm and then assigning each sample to the nearest resulting cluster. As a result, two discrete sequences corresponding to the two input signals A and B are obtained The discrete representation enables the estimation of the joint probabilities of the quantized amplitude values across the two signals. In the present work, Q = 8 quantization bins were used similarly to Ra¨sa¨nen et al. (2013) and Koolen et al. (2014b). 4. In the next step, a so-called energy weighted temporal dependency function (EDTF) between the two quantized envelopes is calculated. EDTF is a mutual information-based metric that measures the overall (non-logarithmic) statistical dependencies across all possible discrete signal state pairs at different relative temporal lags between the signals of interest, weighting the degree of statistical coupling of each state pair by the amount of signal energy associated with these states. More specifically, ETDF is calculated using the following formula: where Ps (a,b) denotes the probability of observing level a in the first channel and level b in the second channel when the second signal is delayed by s seconds relative to the first signal, and where s = [À5, 5] s was used in the present work. The resulting ETDF is offset-normalized to have a minimum-value of zero within this ±5-s range, assuming that the true coupling between the signals should be zero at larger temporal delays and thereby compensating for any estimation biases in the mutual information that could result from the use of finite data. 5. Finally, the ASI is obtained by calculating the ratio of the normalized EDTF value at zero time lag (s = 0) over the mean value of the normalized ETDF function over the 10-s range (s e [À5,5] s), effectively reflecting the magnitude of the zero-delay coupling between the two signals in comparison to delayed or temporally smeared concurrent activation in the signals.
For the present study, ASI was initially computed for EEG epochs (or ASI windows) of 1 min and 2.5 min. It was originally shown (Ra¨sa¨nen et al., 2013), that ASI is more stable with longer (>2 min) window lengths. However, we have later shown that averaging over shorter windows may be preferred (Koolen et al., 2014b), perhaps due to limited long-range temporal correlations in the newborn EEG. Using multiple shorter windows is also better suited for the analysis of older neonates with shorter quiet sleep epochs. As a compromise of the above considerations, we studied the ASI stability (see below) using

Assessment of ASI stability
The underlying general assumption in IHS is temporal stability, hence ASI would be ideally expected to yield relatively stable values regardless of the data length used for its computation. In order to test whether ASI estimates are really so stable, we first analyzed how the length of EEG epoch would affect the findings, and whether the choice of the ASI analysis window influences the observed correlations to the PMA. We found that a significant developmental increase in the ASI is seen for all tested amounts of EEG data, however, there was a slight increase in correlation coefficients when more EEG data was used (Fig. 1). In addition, increasing the ASI window length from 1 min to 2.5 min resulted in ASI vs PMA correlations with both stronger correlation coefficients and steeper slopes (1 min: r = 0.59, slope a = 0.12; 2.5 min: r = 0.79, slope a = 0.26) (Fig. 1). Based on these observations, all later analysis was performed using 2.5 min ASI windows, instead of 1 min windows.

Spatial analysis of ASI
We computed ASI estimates between all monopolar channel combinations. In addition, ASI estimates for bipolar derivations in both hemispheres were obtained for A scheme of these spatial subgroups is shown in Fig. 3. In addition, we assessed another possibility to reduce the number of pairwise ASI estimates by extracting the first component with the principal component analysis.

Network measures
The pairwise connectivity matrix reflects interactions in the network that can also be quantified using graph metrics. Within these metrics, each individual EEG signal can be considered as the node, the signal pair is the edge, and the connectivity measure, here ASI, can be considered the weight of this edge. The commonly used graph theoretical metrics (Bullmore and Sporns, 2009;Stam and van Straaten, 2012) may have limited utility in very sparse graphs, or graphs with highly varying edge levels (see for more details Stam et al., 2014), which characterize our present situation. Hence, we decided to quantify the global network properties with two alternative graph measures: the Minimum spanning tree (MST) and the Algebraic connectivity (AC). The MST is an acyclic sub-network between nodes (here, EEG signals) that allows a quantitative assessment and comparison of networks with Fig. 1. Comparison of ASI analysis settings with respect to developmental correlations. Left side graph shows correlation coefficient (r) between ASI and PMA in the same dataset when ASI is computed using different amount of data (x axis), different analysis windows (1 min vs 2.5 min), or different combinations of channels. A little increase in r-values is found when using longer EEG epochs. On the right side, PMA correlations are shown for 1 min and 2.5 min ASI windows. Both correlations are significant, however use of 2.5 min ASI windows gave clearly steeper developmental trends and higher correlation values.
presumably low bias as compared to more traditional graph metrics (Stam et al., 2014;Tewarie et al., 2015). We computed the metric MST mean using the freely available program BrainWave (http://home.kpn.nl/stam7883/ brainwave.html; Stam et al., 2014). This returns the mean value of all edges selected for the MST of the given individual recording. In contrast to the mean over all 28 pairwise channel connections (measure GS above), MST mean will only take a subset of the strongest connections, which may reduce sensitivity to random variability.
AC is a metric often used to quantify the connectivity; a low AC means that its cost to cut the connectivity graph in approximately two parts is low (cfr. smaller weights between nodes). For maturing babies, edges between the brain regions would be expected to become stronger, leading to more costly graph cuts and, consequently, higher AC values. AC is obtained as the second-smallest eigenvalue of the Laplacian matrix L of the original connectivity matrix, and it is strictly positive (under the assumption of a connected graph) and increasing if more links are added to L. (Bertrand and Moonen, 2013).

Intraindividual stability
We first examined the temporal and spatial variability of ASI by computing the interquartile range (iqr) of all ASI values within each connectivity matrix (spatial variability), and by computing the iqr of the GS in successive ASI windows (temporal variability, Fig. 2B).
Our findings show that spatial variability is large (Fig. 2C), in the order of 1 to 6 ASI units within an individual, but without a significant developmental trend. There were also no significant differences between specific channel combinations (Anova test: 0.96; Fig. 2D). Temporal variability of the GS was expectedly lower, in the order of one ASI unit, and it showed no significant developmental change (Fig. 2B). These findings suggest that the search of developmental correlations is likely more reliable after spatial averaging over regional groups.

Spatial ASI analysis
We observed a clear overall increase in GS values over the course of early development (Fig. 3A). To see whether this was systematically related to specific EEG signal combinations, we assessed the mean ASI of each EEG signal compared to the other seven EEG signals. As shown in Fig. 3B, each of the eight EEG signals showed a significant developmental increase in their mean connectivity, and there were relatively minor differences between individual channels in the slopes of ASI vs. PMA.
The above findings together suggested that the development of activation synchrony is global, and a spatial combination across wider channel groups is possible. We hence computed the median ASI across spatial groups (Fig. 3C), and we found significant developmental correlations for all groups. There were, however, clear differences between spatial groups: the interhemispheric ASI was highest (r = 0.81, slope a = 0.35), followed by the GS (r = 0.79, a = 0.26), while the intrahemispheric ASI was the lowest (r = 0.72, a = 0.26 and r = 0.69, a = 0.18 for left and right hemispheres, respectively). The ASI values of the two premature infants with early IVH grade III were in line with the others (Fig. 3C, patients labeled in gray for GS and interhemispheric ASI).  3. Spatial ASI analysis and its development. (A) In addition to the temporal variability seen as the interquartile range of GS values in successive epochs, there was also an overall increase in GS values with increasing PMA, (B) graphs depicting the developmental change in the mean ASI of each EEG channel compared to the other 7 channels, which have all significant correlations, (C) graphs representing the developmental change in the mean ASI over the given spatial subgroup as schematically shown in the topoplots. Premature infants with early IVH grade III are shown in gray (for GS and interhemispheric synchrony). The right most plot depicts developmental change of the first component of principal component analysis (PCA). Significance of the correlation is depicted with an asterisk after correlation coefficient 'r'. The value 'a' depicts the slope of linear regression computed for the given graph.

Global metrics of ASI-based connectivity
We first tested reduction of analytic dimensions and spatial variability by using the principal component analysis. Its first component was found to have a highly significant correlation (r = 0.77, p = 2.62 * 10 À5 ), as well as a steep relationship to PMA (slope a = 0.88), suggesting that global properties of the connectivity might also reveal robust developmental correlations. Indeed, both measures of global graphs, the mean MST and AC, were found to correlate significantly to PMA (mean MST: r = 0.74, p = 7.2 * 10 À5 ; AC: r = 0.82, p = 3.3 * 10 À6 ) (Fig. 4).

ASI in bipolar derivations
In order to provide an additional, clinically familiar benchmark, we computed ASI between bipolar signals, the preferred montage in the clinical EEG reading (Andre´et al., 2010). We have shown earlier that different bipolar derivations are not directly comparable in terms of their temporal stability (Koolen et al., 2014b) or sensitivity to therapeutic maneuvers (Vanhatalo, unpublished observations), hence they may also exhibit different developmental correlations. We found that frontal interhemispheric connections are stronger correlated to PMA compared to posterior brain regions (Fig. 5). Indeed, no significant correlations are found between PMA and bipolar derivations in the posterior region or in the intrahemispheric channel combinations (r < 0.3 in all cases).
Finally, we wanted to see if there is a systematic spatial asymmetry in the bipolar derivations. Intraindividual comparison of anterior and posterior interhemispheric ASI showed that 20 out of 22 infants had higher ASI values in their anterior areas, which was statistically highly significant (p < 6 * 10 À5 ; binomial statistics). In addition, intraindividual comparisons of left and right intrahemispheric derivations showed that 16 out of 22 infants had higher ASI values on their left hemisphere, which was statistically significant (p < 0.03; binomial statistics).

DISCUSSION
We show that the synchrony between cortical activations strongly correlates with development during the last two months before normal birth in the human. The present quantitative findings are fully compatible with the earlier, mainly qualitative visual EEG observations (Lombroso, 1979). Our work extends prior knowledge by presenting how the developmental change is global, and it reflects the recently shown histological maturation of the corresponding cx-cx networks. It is conceivable that the pace of early functional brain maturation may differ between infants that normally develop in utero vs. those that are born very prematurely and hence undergo longer ex utero development. The effects of ex utero experience on EEG maturation are, however, so small (Shany et al., 2014), that they likely fall within the normal inter-individual variability and the limits of analytical accuracy in our work. A prudent study on the relative contributions of these effects, including various clinical factors, will need larger prospective data collection with repeated EEG measures from each individual to minimize the effects of physiological interindividual variability.
While the age range in this study only extends to term age, the long-range network organization is known to develop many more years into adolescence (Uhlhaas et al., 2010). Studying these later developmental trajectories of brain connectivity will need different analytical methodology, because ASI is only applicable to measure temporal synchrony between intermittent cortical activities, the hallmark of neonatal EEG activity (Vanhatalo and Kaila, 2006;Andre´et al., 2010;Ra¨sa¨nen et al., 2013).
Comparison of brain regions showed that the ASIbased connectivity develops in a global manner, with relatively minor differences between cortical areas. This was somewhat surprising, given the spatial gradients in the development of structural long range cx-cx connections (Judas et al., 2005;Jovanov-Milosˇevicé t al., 2014). Our spatial findings cannot be directly compared to prior studies, because cortical activations comparable to ASI have been studied in central areas only (Meyerson, 1968;Lombroso, 1979;Marcano-Reik et al., Fig. 4. Developmental change of graph metrics, MST mean and algebraic connectivity, both of which showed a significant correlation with PMA. 2010). However, recent work on spatial amplitude relationships suggested that the high cortical activity mode shows developmental gradients and an emergence of frontal and posterior groups toward late gestation (Omidvarnia et al., 2014). Our present findings show an anterior-posterior gradient in ASI without apparent developmental trajectory. This is compatible with the idea that stronger anterior IHS cx-cx correlations may reflect more global functional connectivity in the precentral areas as compared to the more spatially segregated postcentral areas that consist of sensory and association cortices.
Our group level hemispheric comparison revealed a significant asymmetry with left hemisphere showing relatively higher ASI levels in a majority of infants. This ASI asymmetry was computed between two bipolar derivations within the given hemisphere, so it reflects hemisphere level temporal coordination of spontaneous activations. Prior studies have consistently shown that functional and structural hemispheric lateralization begins very early in development (reviewed by Behrmann and Plaut, 2015), however we are not aware of prior EEG studies showing hemispheric differences in connectivity measures in human preterm infants. We find it reasonable to speculate that the ASI asymmetry found in our work reflects the recently reported relative advance in the left side structural connectivity (Ratnarajah et al., 2013).
In addition to spatial averaging of ASI over channel combinations, we showed that graph metrics may disclose similar, significant developmental correlations. However, the limited number of (eight) electrodes available in our dataset did not allow a pertinent assessment of graph measures at the hemispheric or lobar level. Yet, the observations support prior studies (Omidvarnia et al., 2014(Omidvarnia et al., , 2015 in a way that these graph measures may provide a useful tool for the developmental indexing of functional connectivity. Future studies with many more EEG channels for pertinent developmental graph analysis are warranted. While our findings are readily explained in the context of other physiological and anatomical literature, there are some technical considerations that limit the quantitative accuracy of our results. First, the spatial comparisons may be influenced by the number of electrodes and montages. The present dataset consists of eight channel recordings, which is the standard clinical practice (Andre´et al., 2010), but limited with respect to spatial resolution available in the neonatal scalp EEG (Grieve et al., 2004;Odabaee et al., 2014;Tokariev et al., 2015a). We computed most spatial analysis with monopolar montage (Cz reference), which is not fully neutral but likely the best compromise. This may bias findings toward less differences between electrodes, which is likely given the findings from our further assessment with bipolar derivations that showed spatial gradients in both anterior-posterior and left-right direction. The interelectrode distance between monopolar reference and the 'recording electrode' should not have much influence, because of the high spatial specificity shown in the newborn scalp EEG (Odabaee et al., 2013(Odabaee et al., , 2014. Second, more EEG data available for the ASI analysis per patient could reduce the amount of random temporal variability. The amount of high quality EEG data available from human newborn infants is always limited, and our prior work (Koolen et al., 2014b) has shown 5-10 min of EEG yielding reasonably stable values. The possible bias from such technical variability would lead to the underestimation of correlations between ASI and development, which was already found to be highly significant with this data.
In addition to the physiological implications, our work does also suggest that ASI-based quantitation of functional cortical synchrony might offer a useful developmental measure in clinical studies. The recordings used in our present work are available in all medical centers doing routine neurophysiological service for neonates, and we demonstrate here the optimized analysis settings for such datasets. This opens a novel possibility to perform retrospective studies of brain functional connectivity in any developmental disorders that were recorded for clinical reasons during neonatal period. A particular extension of clinical research interest is the possibility to construct developmental indices from the ASI measures. Our work reports many ASI-based measures holding promise as a feature in developmental growth charts, which can be constructed using a well characterized, prospectively collected control population. Such growth charts would provide functional biomarkers that are very much needed in future studies on early development or in association with therapeutic interventions.