However, little information is provided on the developmental trajectories of the functional properties of the sensorimotor cortex, and on how sensorimotor processing is modulated by sensorimotor experience. Here we analyzed the MEG data collected in our previous study (Berchicci et al., 2011) using a measure of functional connectivity and graph theoretical concepts to contribute to the understanding of the functional organization of the developing sensorimotor cortex from early infancy to adulthood. In an attempt to verify whether sensorimotor processing contributes to the generation of perceptual-motor coupling, we analyzed the MEG data recorded in two experimental conditions: at rest and during the execution of a prehension task.

We typified the functional organization over the sensorimotor cortex by means of functional connectivity and efficiency measures. Since the developmental properties of the mu rhythm reflect sensorimotor processing (Pineda, 2005), all measures were calculated within an individual mu rhythm frequency band to ensure that all findings referred to the sensorimotor network. Synchronization likelihood (Stam and van Djik, 2002), used in brain studies to quantify the probability for the functional interdependencies between neural signals, was used to reconstruct functional connectivity maps over the sensorimotor cortex. The type of functional efficiency over the sensorimotor network was estimated by means of segregation and integration measures (Stam, 2000) calculated on the functional connectivity maps. Functional segregation measures (i.e., mean clustering coefficient and local efficiency) quantify the specialized information processing occurring within densely interconnected groups of brain regions (Rubinov and Sporns, 2010), whereas functional integration measures (i.e., characteristic path length, Watts and Strogatz, 1998 and global efficiency, Latora and Marchiori, 2001) reckon the ability of a functional network to rapidly combine specialized information processing from distributed brain regions.

We also tried to explore whether any differences could be observed in the functional properties of the sensorimotor cortex during the execution of a prehension task, in order to provide information useful to understand the contribution of sensorimotor processing to the development of perceptual-motor coupling in infants and children. To this aim, we compared the global measures of functional connectivity and efficiency over the sensorimotor cortex at rest with those related to the execution of a prehension task. Prehension is one of the most remarkable examples of perception and action coupling in infants, since reaching for an object is guided by perceptual information on the relation between the self and the environment, which continuously changes as the posture (e.g., the motor system) is adjusted to that information (Hatwell, 1987). Therefore, prehension is a motor task suitable to explore whether sensorimotor processing contributes to the generation of perceptual-motor coupling.

Our investigation adds to prior studies that have considered either the rest condition or the observation of actions (Gilmore et al., 2011; Virji-Babul et al., 2012; Damaraju et al., 2014; Rotem-Kohavi et al., 2014), since we studied the functional properties over the developing sensorimotor cortex not only at rest but also during prehension. Moreover, by including infants below 6 months of age in our study population, we provided a valuable adjunct to studies on the developing brain.

Based on the knowledge that children and young-adults' brains have "small-world" organization at the global level (Supekar et al., 2009), and on the notion that fronto-temporal connections develop more slowly than other regions (Lebel et al., 2008), we expect that the functional properties over the sensorimotor network evolve, from infants to children and adults, toward a more efficient "small-world" organization. As well, based on the hypothesis that sensorimotor processing plays a role in the generation of perceptual-motor coupling, we expect to observe different functional organizations over the sensorimotor cortex for the rest and prehension conditions.

# **MATERIALS AND METHODS**

### **PARTICIPANTS**

Data collection was performed in infants, children and adults following a cross-sectional design. Participants were selected from a larger database of subjects (see Berchicci et al., 2011). Subjects experiencing any serious illnesses or developmental problems since birth (i.e., traumatic brain injury, seizures, and congenital conditions), or receiving any long-term medication were excluded from the study. Out of the 43 healthy infants enrolled (*<*12 months of age), 25 infants met all inclusion and exclusion criteria. Chronological age at entry in the study ranged between 11 and 47 weeks (mean = 24.9, *SD* = 10*.*8). The functional development of all infants was examined with the Kent Inventory of Developmental Skills (KIDS; Reuter et al., 1996). All infants were found to function within the normal range for age. Parent report was used to assess Apgar score (Apgar, 1953) at birth, which were all within normal limits (i.e., 8–10). However, due to poor signal quality that prevented a reliable calculation of the functional connectivity and efficiency measures, 11 infants had to be excluded from the study (see details in Section Statistical Analysis). Eighteen healthy children aged between 24 and 60 months were enrolled and met all inclusion and exclusion criteria, but, due to motion artifacts, only 12 children were included in the study population (mean = 37.3, *SD* = 12). Six righthanded adults (mean = 28.3, *SD* = 7*.*8) participated in the study as control group.

Based on the results of our prior study (see Berchicci et al., 2011), we knew that a large variation of the mu rhythm peak frequency occurred during the first year of life, indicating a rapid functional development. For this reason, we decided to split the group of infants in two age groups, one including the infants *<*6 months of age, and another including the infants *<*12 months of age. To have a similar number of subjects in all study groups, we split also the group of children in two age groups, one including the children with an age range of 24–34 months, and another including the children with an age range of 36–60 months. Details on the age groups can be found in Section Statistical Analysis.

The protocol was reviewed and approved by the Human Research Review Committee at the University of New Mexico Health Sciences Center and written informed consent was obtained from participant's guardians or adult participants after the description of the study protocol. Infants and children were recruited at day-care centers and from the community using word of mouth, brochures posted on campus, and social networks.

### **PROCEDURE AND TASK**

The experimental setup is only briefly described herein; details can be found in our previous study (see Berchicci et al., 2011). Neuromagnetic activity was recorded using a multi-channel pediatric magnetoencephalography (MEG) system (Okada et al., 2006) for hemispheric recordings. The sensor array operated in a magnetically shielded room, consisted of 76 first-order axial gradiometers, and had a headrest with a smooth outer surface made of thick fiberglass. Inter-sensor distance was approximately 13 mm center-to-center. The headrest was based on a standard reference for the head size of babies. Since the thickness of the neonatal scalp and skull is about 3–4 mm, brain activity could be measured a few millimeters above the brain surface, providing excellent sensitivity.

Since participants performed the assigned tasks with their right hand, they were positioned on the MEG bed on their left side, and a pillow was used to support their back, if necessary. The left hemisphere of the head was positioned on the MEG headrest in order to cover the sensorimotor areas. To ensure child safety and to promote calmness during acquisition, a parent attended the child inside the magnetically shielded room. For each subject, two 5-min blocks of continuous MEG data were recorded. If the infant/child felt uncomfortable from lying still, or was getting drowsy and needed a break, data collection was stopped.

Two different experimental conditions were intermixed: *rest and prehension*. Under the *rest condition*, the participant remained motionless for about 10 s while the investigator stood in front of him/her at a distance of approximately 40 cm. Infants and children were visually engaged to prevent head motions. During the *prehension condition*, participants were invited (in case of adults and children) or directed (in case of infants) to squeeze a pipette placed at about 5 cm from their right hand. The pipette, small enough to be comfortably held and squeezed by infants and children, was connected to a pressure transducer to record the pressure exerted during prehension. Adults used a similar device with a button trigger. The pressure profile was synchronized with MEG data, hence permitting the identification of the time points at which squeezing started. This was necessary for MEG data post-processing purposes. The shielded room light was kept dim during all acquisition sessions to minimize distractions. The experimental sessions were recorded with a video camera synchronized with MEG recording to observe behavior and take note of any movement that might have occurred during acquisitions.

## **DATA ACQUISITION AND PRE-PROCESSING**

MEG data were recorded with a sampling rate of 500 Hz. Although data acquisition was halted in cases of significant displacement of the child's head, residual artifacts, including small movements due to chewing/sucking and arm displacements, were rejected during MEG data pre-processing.

First, MEG channels not working properly were excluded from further analysis. Second, MEG data were band-pass filtered between 0.5 and 40 Hz using a forward-reverse Butterworth filter of the third order. Filtered MEG data were then processed with PCA for dimension reduction, and Independent Component Analysis (ICA, BinICA algorithm with open source toolbox, EEGlab, www*.*sccn*.*ucsd*.*edu/eeglab) was used to separate the independent components (ICs) of interest and to reject those related to artefactual sources. We allowed for a total of 20 ICs, and retained for further processing only those ICs that satisfied two conditions: a topological distribution compatible with the activation of the contralateral sensorimotor cortex, and a clear mu rhythm peak at the individual frequency of the subject, as identified in our previous study (see Berchicci et al., 2011). In general, we retained 2–4 ICs per subject, since the majority of ICs were related to noise, artifacts, MEG channels not working properly, or were seemingly related to brain activity but with a too high content of noise. An example of ICs separated for an infant at 36 weeks is provided in **Figure 1**, where one can see how ICs with a topological distribution that could be related to the activation of the sensorimotor cortex (e.g., IC6 or IC9), have a frequency content that is mainly due to noise. This type of ICs was not retained to reconstruct the brain signals. The retained ICs (ICs 7, 8, and 14 in the example of **Figure 1**) were re-projected on the sensors' positions to reconstruct the MEG source signals related to the true brain activity.

### **DATA ANALYSIS**

#### *Functional connectivity representations*

For each subject, the pre-processed MEG signals were segmented according to the 4-s time intervals of rest and prehension identified in our previous study (see Berchicci et al., 2011) for both spontaneous mu rhythm (rest condition) and mu rhythm desynchronization (prehension condition). The rest time intervals were selected during the resting state far from prehension occurrences. The prehension time intervals were selected around prehension (as identified with the pressure transducer profile and the analysis of the video recordings) to include anticipatory cortical activity.

Hence, prehension time intervals started 1 s before movement onset.

The patterns of functional dependencies between the preprocessed MEG signals were estimated by means of the Synchronization Likelihood (SL) (data were analyzed using Brainwave v<IP_ADDRESS>, http://home*.*kpn*.*nl/stam7883/brainwave*.* html). This measure, introduced by Stam and van Djik (2002), gives a straightforward normalized estimate of the dynamical interdependencies between two or more simultaneously recorded time series. Differently from other connectivity measures such as coherence, this measure, closely related to the concept of generalized mutual information, overcomes the problems related to sub-systems dimensionality, and can be suitable for the analysis of non-stationary data, such as ours. In essence, SL describes how strongly the signal recorded by each channel is synchronized to the signals recorded by the other channels in the array at a given time instant or in a given time span.

For each subject and each experimental conditions (rest and prehension), we calculated one SL matrix for each 4-s time interval within a frequency band centered around the individual mu rhythm peak frequency (IMPF). The IMPF of each subject had been determined in our previous study (see Berchicci et al., 2011), and the frequency band used to calculate the SL matrices had high-pass and low-pass cut-off frequencies at (IMPF − 2 Hz) and (IMPF + 2 Hz) respectively. The SL matrices obtained for all time intervals in each condition were then averaged to obtain one average SL matrix. For each subject we then had two average SL matrices, one for the rest and one for the prehension condition.

To retain only significant functional connections across MEG signals, we thresholded each average SL matrix on the basis of its own SL values distribution, which is expected to be non-Gaussian. We calculated the Median and Median Absolute Deviation (MAD) of the distribution of SL values for each average SL matrix, and defined a new thresholded SL matrix (that we called SLMAD matrix) where only the SL values *>* (Median + 1 MAD) were retained. All other SL values in the SLMAD matrix were set equal to zero.

Since we analyzed the MEG signals in the sensor space, the SLMAD matrices could still have a bias due to the spatial proximity of adjacent sensors. For this reason, we reconstructed the matrix of Euclidean inter-sensor distances MEu for all retained MEG channels, and calculated the coefficient of determination *r*2 between MEu and the two SLMAD matrices (one for rest and one for prehension) for each subject. We then retained for further analysis only the subjects for whom both SLMAD matrices satisfied the condition *r*2 *<* 0.1, hence ensuring that the SL values in the SLMAD matrices do not depend on the Euclidean distances across neighboring MEG sensors.

To quantify the overall probability for functional connectivity over the sensorimotor areas at rest and during prehension, for each subject and each experimental condition we calculated the mean SL value in the SL*MAD* matrix as:

$$SL\_{MEAN} = \frac{1}{N} \sum\_{i,j} SL\_{MAD,i,j} \tag{1}$$

where N is the total number of SL values -= 0 in the SLMAD matrix.

## *Measures of functional organization*

We used graph theoretical concepts to study the topological features of the functional networks represented by the SLMAD matrices. Within this framework, patterns of functional connections are represented as graphs where nodes (in our case the MEG sensors) are linked with edges (if connected). For each subject, we estimated the level of functional organization over the sensorimotor cortex at rest and during prehension by means of segregation and integration measures described below. To this aim, we calculated the segregation and integration measures on the SLMAD matrices, where each SL value (when different from zero) represents the strength of the functional connection (edge) between two given MEG signals (nodes). As the comparison of graphs derived from brain networks requires a step of normalization, for instance by setting a fixed average degree (K) (van Wijk et al., 2010), prior to calculating segregation and integration measures we checked the degrees of the SLMAD matrices of each subject included in the study population, and considered for further processing only those subjects whose SLMAD matrices had the same degree in both rest and prehension (see also sub-Section Statistical Analysis).

*Measures of functional segregation.* Functional segregation in the brain refers to specialized information processing occurring within densely interconnected groups of brain regions (Rubinov and Sporns, 2010). In our case, functional connectivity was calculated by SLMAD matrices in the sensor space, hence segregated neural processing will be suggested by statistical dependencies between clustered channels. We calculated two weighted measures of functional segregation:

(1) *Mean clustering coefficient C*. For a given node *i* (in our case, for a given MEG channel *i*), *ci* is defined as the fraction of the node's neighbors (other MEG channels) that are directly connected with it (Watts and Strogatz, 1998). At the network level, the *mean clustering coefficient C* reflects the prevalence of clustered connectivity around individual nodes:

$$C = \frac{1}{n} \sum\_{i \in N} c\_i = \frac{1}{n} \sum\_{i \in N} \frac{2t\_i'''}{k\_i \ (k\_i - 1)} \tag{2}$$

where *N* is the set of all nodes in the network, and *n* is the number of nodes, *ki* is the degree of node *i* (i.e., the number of links connected to node *i*), and *t w i* is the weighted geometric mean of triangles around *i.* High values of C indicate that a high number of connections exists among neighboring nodes.

(2) *Local efficiency Eloc*. The local efficiency is defined as the average efficiency of the local subgraphs. For weighted connectivity matrices, such as SLMAD matrices, it is:

$$E\_{\text{loc}} = \frac{1}{2} \sum\_{i \in N} \frac{\sum\_{j, h \in N, j \neq i} \left(\boldsymbol{w}\_{i,j} \boldsymbol{w}\_{i,h} \left[\boldsymbol{d}\_{j,h}^{\boldsymbol{w}} \left(\boldsymbol{N}\_{i}\right)\right]^{-1}\right)^{1/3}}{k\_{i} \left(k\_{i} - 1\right)} \tag{3}$$

where, for a given triangle of vertices *i,j,h*, *wi,j*, and *wi,h* are the weights between nodes *i* and *j* and nodes *i* and *h* respectively, and *dw j,h* is the distance between nodes *j* and *h*, i.e., the minimum number of links connecting node *j* to node *h*. *Eloc* provides an information similar to the mean clustering coefficient C, but adds an indication on how much the system is *fault tolerant*, i.e., how efficient is the communication between the first neighbors of node *i* when it is removed. Local efficiency suggests a connection redundancy that protects the network from local errors and failures (Latora and Marchiori, 2001).

When high functional segregation is found in a functional connectivity matrix, the prevailing functional connections occur across neighboring brain areas (or, as in our case, across signals recorded by neighboring channels).

*Measures of functional integration.* In the brain, functional integration represents the ability to rapidly combine specialized information from distributed brain regions. Although our MEG recordings refer only to the hemisphere contra-lateral to the moving limb, we calculated two measures of functional integration to estimate how easily different brain regions communicate. The measures of functional integration are based on the concept of a path:

(1) *Characteristic path length L.* In a graph, the lengths of paths connecting two given nodes estimate the potential for functional integration in the network (Watts and Strogatz, 1998). The *characteristic path length L* is defined as the mean of the distances (i.e., the minimum path lengths) among all nodes pairs:

$$L = \frac{1}{n} \sum\_{i \in N} \frac{\sum\_{j \in N, j \neq i} d\_{i,j}^{\nu}}{n - 1} \tag{4}$$

where *dw i,j* is the weighted distance between nodes *i* and *j*, i.e., the minimum number of links connecting node *i* to node *j*. Short paths, corresponding to low values of L, indicate a strong potential for integration within the network.

(2) *Global efficiency Eglob*. In functional connectivity data, such as the SLMAD matrices, paths represent sequences of statistical associations between subsequent pairs of channels. The level of functional integration in the network can be represented by the *global efficiency Eglob*, which is the average inverse shortest path length (Latora and Marchiori, 2001; Achard and Bullmore, 2007). For weighted connectivity matrices:

$$E\_{glob} = \frac{1}{n} \sum\_{i \in N} \frac{\sum\_{j \in N, \ j \neq i} (d\_{i,j}^w)^{-1}}{n - 1} \tag{5}$$

When high functional integration is found in a functional connectivity matrix, the functional organization of the brain takes advantage of multiple specialized and densely connected areas that are linked with long distance functional connections for a more efficient information processing.

#### *Statistical analysis*

Linear regression analysis was performed on the individual measures to assess whether any linear correlation exists between age and/or mu rhythm peak frequency and the functional connectivity/efficiency measures.

Groups analysis was performed for each connectivity/efficiency measure on the age groups defined in Section Participants. In each age group, we retained only those subjects for whom the following two conditions were satisfied: (1) the determination coefficient r2 between the individual SLMAD matrices at rest and during prehension and the individual MEu matrix of the euclidean inter-sensor distances was smaller than 0.1, and (2) the degree of the individual SLMAD matrices for rest and prehension was the same. The first condition guarantees that the SLMAD matrices are independent on the inter-sensor distances, and the second condition satisfies the requirement for a reliable comparison of segregation and integration measures. Seven infants (age *<*6 months) and 4 infants (age *<*12 months) did not satisfy these conditions and were excluded from further analysis. The characteristics of the age groups are reported in **Table 1**.

Statistical analysis was performed using an ANOVA 5 (age groups) × 2 (conditions: rest, prehension) for each dependent variable, i.e., for the individual measure of functional connectivity (*SLMEAN*), and for functional segregation (*C* and *Eloc*), and integration (*L* and *Eglob*) measures. Post-hoc comparisons were performed using Bonferroni corrections. Statistical significance was set at *p <* 0*.*05.

#### **RESULTS**

Linear regression analysis on the individual functional connectivity/efficiency measures did not show any significant linear dependence on age or mu rhythm peak frequency. In particular, the coefficient of determination r2 between age and *SLMEAN*, *C*, *Eloc*, *L*, and *Eglob* ranged between 0.197 and 0.465, and r2 between mu rhythm peak frequency and *SLMEAN*, *C*, *Eloc*, *L*, and *Eglob* ranged between 0.008 and 0.269, indicating no significant linear correlation between any connectivity/efficiency measure and age or mu rhythm peak frequency.

#### **FUNCTIONAL CONNECTIVITY REPRESENTATIONS**

The changes of the group-averaged SLMEAN values with age in the two experimental conditions (rest and prehension) are summarized in **Figure 2**. ANOVA results on the individual measures of SLMEAN showed significant differences across age groups [*F*(4*,* 27) = 4*.*546 *p* = 0*.*006, η2 p = 0*.*402, power = 0.897], no significant differences between conditions, and no interaction between groups and conditions. During the first year of life (groups G1 and G2), the average SLMEAN values do not differ significantly, whereas around 2 years of age (group G3), SLMEAN values start to increase toward the adult values (group G5). *Post-hoc* analysis showed significant differences only between the

**Table 1 | Number of subjects included in each age group (Nstudy) out of those who met all inclusion and exclusion criteria (Nincl***/***excl).**

For each age group, the age range and the mu rhythm peak frequency (Mean ± SD) are also reported.

group of adults G5 (0.315 ± 0.093) and the groups of infants G1 (0.143 ± 0.056; *p* = 0*.*015) and G2 (0.129 ± 0.062; *p* = 0*.*007).

These results indicate that, during infancy, the likelihood for functional connections over the sensorimotor cortex is much lower than in adulthood, and that it increases with age. Although no statistically significant difference was observed between conditions, it is interesting to observe that, in our group of very young infants (G1), the average SLMEAN is slightly higher during prehension than at rest, whereas in the groups of children and in the group of adults the average SLMEAN slightly decreases during prehension. This fact could deserve further investigation, possibly in larger populations of infants and children. Examples of SLMAD matrices at rest and during prehension for each age group are provided in **Figure 3**. A typical layout of the MEG sensor array over the head is also shown. However, please consider that the same channels in different SLMAD matrices could be positioned over different brain regions, as no information on the relative positions of the MEG channels with respect to the subject's head could be collected.

#### **MEASURES OF FUNCTIONAL SEGREGATION**

The changes of the group-averaged *C* and *Eloc* values with age in the two experimental conditions (rest and prehension) are summarized in **Figures 4A,B**. ANOVA results on the individual measures of *C* showed significant differences across age groups [*F*(4*,* 27) = 4*.*975 *p* = 0*.*004, η2 p = 0*.*424, power = 0.924], no significant differences between conditions, and no interaction between groups and conditions. During the first year of life (groups G1 and G2), *C* remains almost unchanged, whereas around 2 years of age (group G3) it starts to increase toward the adult values (group G5). *Post-hoc* analysis showed significant differences between the group of adults G5 (0.231 ± 0.076) and the groups of infants G1 (0.094 ± 0.037; *p* = 0*.*011), G2 (0.083 ± 0.047; *p* = 0*.*005). A tendency toward significant difference was also observed between G5 and G3 (0.104 ± 0.041; *p* = 0*.*030).

Similarly, ANOVA results on the individual measures of *Eloc* showed significant differences across age groups [*F*(4*,* 27) = 5*.*223 *p* = 0*.*003, η2 p = 0*.*436, power = 0.937], no significant differences between conditions, and no interaction between groups and conditions. During the first year of life (groups G1 and G2), *Eloc* remains almost unchanged, and around 2 years of age (group G3) it starts to increase toward the adult values (group G5). *Posthoc* analysis showed significant differences between the group of adults G5 (0.266 ± 0.084) and the groups of infants G1 (0.112 ± 0.041; *p* = 0*.*008), G2 (0.100 ± 0.049; *p* = 0*.*004), whereas a tendency toward significant difference was observed between G5 and G3 (0.125 ± 0.045; *p* = 0*.*026).

These results indicate that, during infancy, global measures of functional segregation over the sensorimotor cortex are significantly lower than in adulthood, and that they increase with age. Although no statistically significant difference was observed between conditions, in the very young infants (G1), functional segregation seems to be higher during prehension than at rest, whereas in children and in adults it seems to decrease from rest to prehension. As for SLMEAN, this aspect deserves further investigation in larger populations of infants and children.

#### **MEASURES OF FUNCTIONAL INTEGRATION**

The changes of the group-averaged *L* and *Eglob* values with age in the two experimental conditions (rest and prehension) are shown in **Figures 4C,D**. ANOVA results on the individual measures of *L* showed significant differences across age groups [*F*(4*,* 27) = 4*.*245 *p* = 0*.*009, η2 p = 0*.*386, power = 0.873], no significant differences between conditions, and no interaction between groups and conditions. *L* tends to decrease with age, although discontinuously. *Post-hoc* analysis showed significant difference between the group of adults G5 (12.975 ± 3.034) and the group of older infants G2 (28.588 ± 7.575; *p* = 0*.*007), and a trend toward significant difference between G5 and the group of older children G4 (25.739 ± 10.494; *p* = 0*.*047).

ANOVA results on the individual measures of *Eglob* showed significant differences across age groups [*F*(4*,* 27) = 5*.*594 *p* = 0*.*002, η2 p = 0*.*453, power = 0.952], no significant differences between conditions, and no interaction between groups and conditions. We can see that *Eglob* remains almost unchanged in infancy (groups G1 and G2), and that it increases with age. *Post-hoc* analysis showed significant differences between the group of adults G5 (0.158 ± 0.037) and the groups of infants G1 (0.071 ± 0.022; *p* = 0*.*004) and G2 (0.068 ± 0.029; *p* = 0*.*003). A tendency toward significant difference was observed between G5 and G3 (0.085 ± 0.027; *p* = 0*.*030).

Both measures indicate a significant difference of functional integration over the sensorimotor cortex between infancy and adulthood, and that functional integration over the sensorimotor cortex increases with age. Although no statistically difference was observed between conditions, both *L* and *Eglob* seem to indicate a higher functional integration during prehension in the very young infants (G1), whereas in children and in adults functional integration seems to decrease from rest to prehension. As for the other connectivity and segregation measures, this aspect deserves further investigation in larger populations of infants and children.

# **DISCUSSION**

Our goals were to contribute to the identification of the developmental trajectories of the functional properties over the sensorimotor cortex, and to explore whether any differences could be observed between the rest condition and the execution of a prehension task. To ensure that all findings were related to the sensorimotor cortex, all analyses were performed within a frequency band centered on the individual mu rhythm peak frequency of each subject, as defined in our previous study (Berchicci et al., 2011), since the mu rhythm is known to be the idling rhythm of the sensorimotor cortex. Furthermore, the use of a prehension task for the active condition permitted us to investigate the contribution of sensorimotor processing to the generation of perception and action coupling in infants.

We first observed that the probability of functional connectivity across the sensorimotor areas, estimated by means of *SLMEAN*, remained almost unchanged during the first year of life, whereas it afterwards increased across groups to reach adult values, significantly different from the infant ones. This increase, occurring after 12 months of age, is in agreement with the suggested pattern of brain maturation in which areas with fronto-temporal connections develop more slowly than other regions (Lebel et al., 2008), and with the notion that motor experience increases the capacity of the whole network to work in an integrated way (Fair et al., 2009; Giedd et al., 2009; Boersma et al., 2011). Unfortunately, no information on the relative position of the MEG sensor array and the subject's head could be collected, therefore we could not precisely ascribe the observed functional connections across MEG signals to specific brain areas. Consequently, our observations on the probability of functional connectivity should be interpreted as global measures over the sensorimotor network.

To support the interpretation of these results, we used graph theoretical concepts to examine the changes occurring in the functional organization of the sensorimotor network across ages, and found that all segregation and integration measures showed a trend similar to *SLMEAN*: during the first year of life these properties remained almost unchanged, whereas they increased significantly across groups afterwards, as assessed by a marked increase of both segregation parameters (*C* and *Eloc*), by a tendency of *L* to decrease, and by a clear increase of *Eglob* across age groups.

The age-related increase of the functional segregation parameters (*C* and *Eloc*) is compatible with the observation that the widespread brain electrical activity typical of infants (8 months of age) becomes more localized during early childhood (3 years of age) (Bell and Wolfe, 2007). We hypothesize that the increase of functional connections among adjacent cortical areas, which contributes to the functional specialization of the sensorimotor network and to the onset of a connection redundancy that protects the network from local errors and failures (Latora and Marchiori, 2001), is due to brain maturation, a process that involves both physical growth and the intellectual and/or emotional process of

development. Also, the increase of local efficiency over the sensorimotor cortex is in line with the findings by Wu et al. (2013), who observed an age-related increase in the local efficiency of the whole functional networks, which may contribute to the development of modular information processing of functional systems.

Although the findings on the functional integration properties over the sensorimotor cortex should be taken with caution because the MEG sensor array covered only one hemisphere, it is interesting to observe that also these properties increased with age after the first year of life, and that the infant values were significantly different from the adult ones. The irregular trend of *L* is also worth noting, as it might be related to the occurrence of growth spurts that coincide with periods of discontinuous development in cognition (van Baal et al., 2001).

Overall, our findings on the functional segregation and integration properties over the sensorimotor cortex support the notion that maturation contributes to both the functional specialization of the sensorimotor network, and the wiring of more efficient long-range connections among different brain circuits. This interpretation is in agreement with the observations by Gao et al. (2011) on the increase of both local and global efficiency at rest from newborns to infants, and with the findings of Damaraju et al. (2014) that the connectivity strength between more distant networks increases with age in infants between 4 and 9 months. Conversely, our results are only partially compatible with the notion that strong local connectivity in the young brain gradually shifts toward stronger long-distance connectivity with maturation (van Baal et al., 2001; Lebel et al., 2008; Fair et al., 2009; Power et al., 2010; Yap et al., 2011).

Our results on the age-related increase of both local and global efficiency are consistent also with the outcome of other studies on the developing brain, which demonstrated a steady increase of global efficiency after 1 year of age (Fransson et al., 2007; Fan et al., 2011), a positive relationship between age and local efficiency throughout the life span (Supekar et al., 2009; Dennis et al., 2013; Wu et al., 2013), and the involvement of both segregation and integration in the development of the adult frontoparietal network for adaptive online task control (Fair et al., 2007; Sepulcre et al., 2010). Our findings seem to support this interpretation of the functional evolution of the fronto-parietal network, known to host sensorimotor processing (Pineda, 2005). From this perspective, our results seem to fit well also with the suggestion that perceptual and cognitive developments involve the simultaneous segregation and integration of informationprocessing streams (Bunge and Wright, 2007), which in turn support the learning processing in infants and children (de Klerk et al., 2014).

Some authors have also suggested that the child brain has a small-world organization that combines the high clustering properties of an ordered network and the short path length of a random network (Fair et al., 2009; Supekar et al., 2009; Power et al., 2010; Boersma et al., 2011). Other authors have observed that, in school children, the brain functional networks evolve from more random toward a more ordered configuration (Smit et al., 2010, 2012; Boersma et al., 2011). Although our hemispheric recordings do not allow for general conclusions on the functional organization of the whole brain, nonetheless our results complement the above mentioned findings by suggesting that, during infancy, the functional properties over the sensorimotor network have a more random organization that, after 1 year of age, shifts toward a more efficient small-world configuration characterized by higher local and global efficiency. It is interesting to observe that this change toward a more efficient functional organization over the sensorimotor network continued throughout childhood until adulthood. The shift from a more random to a small-world configuration could be ascribed to brain maturation, which derives from changes at the structural level, with the thickening of long-range connections among distant brain regions starting during the first year of life (Yap et al., 2011) and with the emergence of inter-hemispheric connectivity earlier than longer-range antero-posterior connections (Bell and Fox, 1992), but is also shaped by experience, which can support the capacity of the network to create long-range connections that are functional to goal-directed actions (Wu et al., 2013; Rotem-Kohavi et al., 2014).

Unfortunately, we could not draw any conclusions on the contribution of sensorimotor processing to the generation of perceptual-motor coupling, since we did not find any significant difference in the functional properties over the sensorimotor cortex between rest and prehension. However, it is worth noting that the pattern of functional organization over the sensorimotor cortex seems to change around 1 year of age. All measures seem to suggest that, during infancy, prehension could be characterized, with respect to rest, by an increased probability of functional connectivity and by increased segregation and integration. This pattern seems to reverse after 1 year of age. These observations, if confirmed in future studies on larger populations of infants and children supported by measurements of motor experience, could indicate that, in young infants, perceptual-motor coupling is accompanied by an increased connectivity within the sensorimotor network to compensate for the lack of functional specialization to accomplish the given task (prehension). The evolution of this functional pattern could be due to the developmental timeline of the infant grasping or to the tight correlation between grasping recognition and execution, which are poor until 6 months of age (Del Giudice et al., 2009), hence requiring the recruitment of more functional resources in early infancy with respect to later ages. On the other hand, the pattern reversal that seems to occur after 1 year of age might indicate that maturation had fostered a more efficient execution of goal-directed actions, in line with the experience-dependent position on perceptual-motor coupling development of the sensorimotor cortex (Heyes, 2001; van Elk et al., 2008; Del Giudice et al., 2009; Cook et al., 2014), and with the concept that sensorimotor experience contributes to the development of more efficient functional networks to support sensorimotor development (Fransson et al., 2011).

We are aware that the impact of our study is limited by a number of factors. First, the shape and size of the MEG sensor array allowed the recording of brain activity originating only from one hemisphere. Consequently, our observations on the smallworld topology of the child brain cannot be referred to the whole brain, and the interpretation of the integration measures needs to be cautious. Second, we employed a cross-sectional design for this study. Given that, during infancy and childhood, the brain undergoes a remarkable development at both structural and functional levels, the recorded brain activity could be influenced by the size of the scalp at different ages, which fits differently on the headrest covering the sensor array. Third, the measures of the topological features of the developing sensorimotor network were calculated in the sensor space rather than in the source space, because anatomical information on the baby's head and on its position with respect to the MEG sensor array could not be collected. This condition limits the interpretation of our results in terms of specific cortical areas. Fourth, we are aware of the small sample size of our age groups. However, study protocols like ours are difficult to perform in children and even more so in infants, and the analyses performed are extremely sensitive to noise and to other signal features that further reduced the number of retained subjects. Nonetheless, several research studies in adults and children have been considered reliable even when based on a small number of subjects (Nishitani and Hari, 2000; Simoes et al., 2004; Lepage and Theoret, 2006; van Schie et al., 2008).

In light of these limitations, further work is needed to validate our findings. However, we believe that our results add valuable information to the current knowledge on the functional development of the sensorimotor network during infancy and childhood.  Better understanding of the functional development of these networks during infancy may bring new insights on the pathophysiological mechanisms of neurological development, such as Down syndrome, autism spectrum disorders and cerebral palsy. Further studies employing whole head neuroimaging techniques during motor task execution and complemented with motor experience measurements are therefore needed to support the present findings, and to improve our understanding of the functional development of the sensorimotor network and its contribution to the generation of perceptual-motor coupling in infants.

# **AUTHOR CONTRIBUTIONS**

The study presented here has been conceived and designed by Marika Berchicci and Silvia Comani. The experiments were performed by Marika Berchicci. Data were analyzed by Gabriella Tamburro, Silvia Comani, and Marika Berchicci. Data have been interpreted by Marika Berchicci, Gabriella Tamburro, and Silvia Comani. The manuscript was written, revised and approved by Marika Berchicci, Gabriella Tamburro, and Silvia Comani.

# **REFERENCES**

late childhood to early adulthood. *Cereb. Cortex* 19, 640–657. doi: 10.1093/cercor/bhn117

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

*Received: 24 November 2014; accepted: 14 January 2015; published online: 17 February 2015.*

*Citation: Berchicci M, Tamburro G and Comani S (2015) The intrahemispheric functional properties of the developing sensorimotor cortex are influenced by maturation. Front. Hum. Neurosci. 9:39. doi: 10.3389/fnhum.2015.00039*

*This article was submitted to the journal Frontiers in Human Neuroscience.*

*Copyright © 2015 Berchicci, Tamburro and Comani. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Mothers say "baby" and their newborns do not choose to listen: a behavioral preference study to compare with ERP results

#### Christine Moon1 \*, Randall C. Zernzach2 and Patricia K. Kuhl 3

1 Department of Psychology, Pacific Lutheran University, Tacoma, WA, USA, 2 Wilford Hall Ambulatory Surgical Center, 59th MDW, JBSA-Lackland, TX, USA, 3 Institute for Learning and Brain Sciences, University of Washington, Seattle, WA, USA

Previously published results from neonatal brain evoked response potential (ERP) experiments revealed different brain responses to the single word "baby" depending on whether it was recorded by the mother or an unfamiliar female. These results are consistent with behavioral preference studies in which infants altered pacifier sucking to contingently activate recordings of the maternal vs. an unfamiliar female voice, but the speech samples were much longer and information-rich than in the ERP studies. Both types of neonatal voice recognition studies imply postnatal retention of prenatal learning. The preference studies require infant motor and motivation systems to mount a response in addition to voice recognition. The current contingent sucking preference study was designed to test neonatal motivation to alter behavior when the reward is the single word "baby" recorded by the mother or an unfamiliar speaker. Results showed an absent or weak contingent sucking response to the brief maternal voice sample, and they demonstrate the complementary value of electrophysiological and behavioral studies for very early development. Neonates can apparently recognize the maternal voice in brief recorded sample (previous ERP results) but they are not sufficiently motivated by it to alter sucking behavior.

#### Edited by:

Marika Berchicci, University of Rome "Foro Italico" - Rome, Italy

#### Reviewed by:

Christine Parsons, University of Oxford, UK Theofanis Panagiotaropoulos, Max Planck Institute for Biological Cybernetics, Germany

#### \*Correspondence:

Christine Moon, Department of Psychology, Pacific Lutheran University, 12180 Park Ave. S., Tacoma, WA 98447, USA<EMAIL_ADDRESS>
> Received: 24 November 2014 Accepted: 06 March 2015 Published: 25 March 2015

#### Citation:

Moon C, Zernzach RC and Kuhl PK (2015) Mothers say "baby" and their newborns do not choose to listen: a behavioral preference study to compare with ERP results. Front. Hum. Neurosci. 9:153. doi: 10.3389/fnhum.2015.00153 Keywords: auditory, behavior, fetal, learning, motivation, neonatal, non-nutritive sucking, preference

Recent electrophysiological and brain imaging studies have advanced our understanding of the development of very early perception and learning through measurement of brain activity in response to specific events (May et al., 2011; Partanen et al., 2013; Kuhl et al., 2014). When there are differential responses to familiar vs. novel stimuli, learning and memory can be inferred.  It is, of course, ideal when brain activity can be linked to behavior, and the perception-action circuit is illuminated. Documenting the link can be difficult in early development because the competence-performance distinction (Chomsky, 1965) is particularly relevant. If the brain shows that a stimulus is recognized, there may be no corresponding measurable behavior because the individual cannot yet execute a motor response. A second consideration is that the path to behavior includes not only motor competence but also motivation. The newborn brain may detect an event, recognize it, and may be capable of rendering a motor response, but if the event is not sufficiently motivating, the infant will not mount the response.

Results of several experiments show that the newborn brain responds differentially to even a brief snippet of the maternal vs. a stranger female voice (Deregnier et al., 2000; deRegnier et al., 2002; Siddappa et al., 2004; Therien et al., 2004). In several studies, scalp recordings of brain activity in sleeping newborns using evoked response potentials (ERPs) used the brain's mismatch detection response to show this discrimination. This procedure requires that the stimuli be brief so the brain's rapid response can be time-locked to stimulus presentation. In the newborn ERP studies, the mother and stranger voice recordings were each limited to a brief sample of the word ''baby'' that averaged 750 ms. The midline brain response to both the maternal and novel voices was a positive wave appearing at 290 ms and that has been interpreted as reflecting auditory stimulus detection of both types of stimuli. For the novel voices only, the positive wave was accompanied by a negative slow wave that was interpreted as a response to novelty (deRegnier, 2005). The lack of a similar response to the maternal voice implies recognition of a familiar sound. The recognition response has been detected in newborns from uncomplicated pregnancies and births, whereas in newborns of iron-deficient (Siddappa et al., 2004) mothers and in extremely premature newborns (Therien et al., 2004), the differences are attenuated or absent. These ERP differences were interpreted as showing compromised recognition memory development in the infants of complicated pregnancies or deliveries.

Prior to these ERP results, it was not known whether newborns could even recognize their mothers' voices without the prosodic information present in lengthy and acoustically rich samples of talking or reading. Doubt about the limits of newborn voice recognition arose from analysis of intrauterine voice recordings (Querleu et al., 1988) and from behavioral experiments using only lengthy and acoustically rich voice samples. For a review, see Moon and Fifer (2000). It is now apparent that brief and relatively uninformative voice samples are sufficient for a differential brain response.

There have been several behavior-based investigations of maternal voice recognition. DeCasper and Fifer showed that neonates alter their behavior to selectively activate the sound of the mother reading a nursery rhyme (DeCasper and Fifer, 1980). Because the results were based on infants who had had no more than 12 h postnatal experience with their mother, their selective response to the maternal voice implied postnatal retention of prenatal learning.  The authors compared the observed behavior to that of non-human infants who express ''perceptual preference'' and ''proximity-seeking behavior'' p.1176 (DeCasper and Fifer, 1980). Developmental psychobiologists and cognitive development researchers have described this kind of behavior as a ''listening bias'' (Vouloumanos and Werker, 2004), ''operant choice'' (Granier-Deferre et al., 2011), ''operant learning'' (Floccia et al., 1997) and ''operant-choice preference'' (Aldridge et al., 2001). The latter terms presuppose differential reinforcing value of the stimuli and motivation to respond. A contingent sound-sucking experiment using yoked-control methodology showed that the contingency is important in increasing sucking rates (Floccia et al., 1997). Thus, a measurable difference in neonatal preference or choice offers the opportunity to study motivation at the beginning of life.

Since 1980, experiments using contingent pacifier sucking have replicated and extended two main ideas: (i) humans learn before birth; and (ii) prenatal learning can affect postnatal behavior. The maternal voice preference has been robust in demonstrating learning in the womb. It has held whether the voice samples were highly melodic and rhythmic recordings of mothers reading nursery rhymes (DeCasper and Fifer, 1980) or relatively monotone adult conversation (Fifer and Moon, 1989). Although neonates respond to a simple maternal voice recording, they prefer a version that has been low pass filtered to simulate the voice in utero (Fifer and Moon, 1995; Spence and Freeman, 1996). A newborn experiment with fathers' and comparison male voices failed to show a preference, despite newborn discrimination of the two voices (DeCasper and Prescott, 1984). During the neonatal period, infants have expressed preferences for other sounds available in utero such as the mother's native vs. a foreign language (Mehler et al., 1988; Moon et al., 1993). The native over foreign language preference was replicated and extended in a study showing that, for infants of bilingual mothers, her two languages receive equivalent responses (Byers-Heinlein et al., 2010). Perhaps the most convincing example of prenatal learning's effect on postnatal preference are the results of an exposure study in which mothers read a nursery rhyme out loud during pregnancy. Newborns sucked more to activate the familiar rhyme vs. a novel rhyme, regardless of whether the voice was maternal (DeCasper and Spence, 1986). Thus, there is converging evidence from 35 years of laboratory research that shortly after birth infants are capable of employing prenatal experience, motor control, and the motivational system to mount a behavioral response when the consequence is hearing a familiar sound.

What has not been known up to now is whether the brief and recognizable sample of the maternal voice saying a single word is sufficiently motivating for newborns to mount a preference response. What follows is a description of an experiment using the preference procedure in which the consequence of infant sucking on a pacifier was the delivery through headphones of brief maternal and non-maternal voice samples of the word ''baby''. Thirty-six infants were in the experimental group of Mother-Stranger (MS) infants for whom the stimuli were the maternal vs. a stranger female voice. Twenty-four infants were in the Stranger-Stranger (SS) control group who heard two different non-maternal female voices. If the brief snippets of mother's voice are sufficiently motivating for newborns to alter sucking behavior to show a preference, then (Kuhl et al., 2014) sucking during opportunities to produce the maternal voice should exceed those for the stranger voice, and (May et al., 2011) sucking by the MS group should exceed that of the SS group.

# Materials and Methods

# Participants

Sixty term neonates (M age = 31 h, SD = 10.3) each completed an experimental session prior to discharge from the mother-baby postpartum unit of a military medical center. They had no documented antenatal or birth complications, no risk factors for hearing loss (American Academy of Pediatrics Joint Committee on Infant Hearing, 2007), and English was the primary language spoken in the home. Infants were assigned to one of two conditions. In the MS Condition (N = 36), the stimuli were mother's and a stranger female's voice. In the Stranger-Stranger Condition (N = 24), the stimuli were two voices of unfamiliar females. In order to be included in the analysis, infants were required to complete a 10 min session with no more than two consecutive or three non-consecutive minutes in which no sucking occurred. Study sessions were terminated immediately if the infant became excessively fussy or cried. The data from 29 infants were excluded from analysis due to drowsiness (n = 18), crying/fussiness (n = 3), inconsistent sucking not apparently due to arousal state (n = 7), or experimenter error/equipment problems (n = 1).

# Apparatus and Stimuli

Infants sucked on a Gerber Little Suzy Zoo pacifier fitted with a plastic tube connected to a Becton-Dickinson P23 × L pressure transducer that provided input to a Grass Telefactor CP122 Strain Gage Amplifier. The analogue output of the amplifier was converted to a digital signal by a Data Translation DT2814 data acquisition board that was connected to a Gateway 486 PC equipped with a ProAudio 16 sound card. Custom software recorded sucking pressure, controlled stimulus delivery, and created a summary data file. Analog pressure changes were converted into digital signals and served as input for a Gateway 486 PC. The computer delivered stimuli to the infants via Grado SR225 earphones suspended from a custom-made adjustable plexiglass frame that fit into the infant's bassinet. Voice stimuli were recorded using an Electro-Voice PL88H Microphone on a laptop computer. Analysis of the recorded speech samples was conducted using Signalyze Speech and Sound Analysis software on an Apple computer.

Each of 53 stimuli consisted of one naturalistic token of the word ''baby'' spoken in a woman's voice. Talkers had been instructed to say the declarative sentence, ''He's/she's a baby'' with a falling intonation contour. The statement was followed by four repetitions of the word ''baby'' using the same intonation contour. One voice token for each speaker was chosen on the basis of absence of background sound, clarity, loudness and presence of falling intonation contour. The mean duration of voice stimuli was 480.0 ms (s.d. 92.3). For the MS group, the Mother and Stranger stimuli did not significantly differ in duration, nor did the two voices in the SS condition. Pairs of stimuli were matched for loudness by two adult listeners and were presented at about 66 dBA as measured by a sound level meter (Bruel and Kjaer Model 2235) placed midway between the headphones.

# Design and Procedure

The study was conducted according to a protocol that was approved by the medical center and university (PLU) ethics boards. For the MS infants, the stimuli were available during five one-minute periods of stimulus voice alternations for a total of 10 min (Cowan et al., 1982; Sansavini et al., 1997; Vouloumanos and Werker, 2007). The MS stimuli were counterbalanced for order of presentation (Mother First, N = 17, Stranger First, N = 19) and were delivered on a partial reinforcement schedule of a minimum of two qualifying sucks in quick succession for reinforcement. The Stranger-Stranger infants (SS) heard two different unfamiliar female voices, the availability of which alternated in one-minute intervals for 10 min. Each voice was thus available for five one-minute periods.

Parental informed consent was obtained in the motherinfant hospital room according to university- and medical center-approved protocols. For infants in the MS group, a recording of maternal voice was obtained at the bedside for subsequent editing. Each mother received a walkie-talkie to call when her infant appeared to be in a quiet and alert state. Infants were transported to the study room in the motherbaby unit where the headphone frame was placed in the bassinette and headphones were fitted next to the infants' ears. Prior to the beginning of the session, an experimenter swaddled the infant either snugly or loosely, based on the best judgment about supporting a quiet and alert state for 10 min. The experimental pacifier was given to infants and was held in place by an experimenter who listened to masking music through headphones during the session. After subjects demonstrated consistent sucking on the pacifier, usually no more than two to three minutes, the study session began by recording 1 min of baseline sucking data at the conclusion of which a threshold was established for stimulus delivery. The threshold was set at the 30th percentile of sucking amplitudes produced during the baseline minute.  This is consistent with previous neonatal contingent sucking studies in which amplitude threshold has ranged from the 20th (Floccia et al., 2000; Vouloumanos and Werker, 2007) to the 50th (Floccia et al., 1997) percentiles of baseline.

# Results

The dependent measure was the number of qualifying (above threshold) sucks per minute. Analyses included the factor of time (sequential minutes of the session) because previous studies have shown the emergence of a preference response over time (Sansavini et al., 1997; Vouloumanos and Werker, 2004, 2007). A D'Agostino-Pearson test of normality was conducted on the mean number of sucks per minute for the sample of 60 newborns, and results were consistent with a normal distribution, K 2 = 2.21, p > 0.05.  There was a significant difference with greater mean sucks per minute of baseline for the MS group, M = 43.0, SD = 17.6 compared to the SS group, M = 31.9, SD = 19.2, t(58) = 2.3, p = 0.025. Baseline sucks were entered as a covariate in subsequent analyses.

voices alternated in five one-minute intervals of a 10 min presentation period.

TABLE 1 | Analysis of variance results for the Mother-Stranger group (N = 36).

Data from the MS group were used to test the hypothesis that sucks during minutes of opportunity to hear mother's voice would exceed those of stranger voice minutes. For the MS group (N = 36) a mixed 2 (**Voice**) × 2 (**Order** of voice presentation) × 5 (**Minutes**) ANOVA was conducted. There were no statistically significant main effects of the within-subjects variable Voice (Mother Voice sucks per minute: Mean = 28.7, SD = 16.3, Stranger Voice Mean = 27.0, SD = 14.4 or the betweensubjects variable of Order. There was no main effect of Minutes nor were there any significant interaction effects. See **Figure 1**; **Table 1**.

 A mixed two factor ANOVA included Stimulus Group (2) × Minutes (10). There was no main effect of Stimulus Group (F(1,57) = 1.9, p = 0.18, ηp 2 = 0.03), no main effect of Minutes (F(1,513) = 1.33, p = 0.22, ηp 2 = 0.02) and no interaction effect of Group X Minutes (F(9,513) = 0.89, p = 0.22, ηp 2 = 0.12). See **Figure 2**.

# Discussion

Neither the direct comparison of responses to the maternal vs. stranger voices nor the indirect comparison in present vs. absent maternal voice provides statistically significant evidence for a preference for the maternal voice. Visual inspection of the pattern of the two MS sub-groups in **Figure 2** provides some evidence for a preference response in the transitions mother-tostranger (Mother First group) and stranger-to-mother (Stranger First group) over 10 min.

One possible explanation for the absent or weak behavioral response to mother's voice that can be ruled out is that infants could not recognize the maternal voice in the brief sample. Previously published ERP studies demonstrate otherwise, at least in infants without complicated prenatal and birth histories.

Although the newborns were able to discriminate the maternal and the stranger voices, the brief, repetitive sample of the maternal voice was apparently not sufficiently motivating for them to suck to activate it significantly more frequently than the alternative. There are other published reports of newborn failure to show a voice preference in contingent sucking procedures.  Newborns did not show a maternal voice preference when mother's and stranger voices were whispered (Spence and Freeman, 1996) and they showed no preference when hearing recordings of their bilingual mothers' two languages (Byers-Heinlein et al., 2010).

Collecting behavioral data from neonates has limitations. There are many sources of variability that are difficult to control such as rapid changes in infant state, competing infant behaviors, and internal perceptual events that are often unspecifiable. In the current study, the attrition rate was 33 per cent. Although this rate is typical or even low for contingent sucking experiments (Floccia et al., 1997, 2000; Vouloumanos and Werker, 2007), it is an indicator of the inherent variability in behaviorbased data collection with neonates. Although the current study was conducted by experienced neonatal researchers using a standard contingent sucking protocol, it certainly merits replication, at best with extension to other forms of familiar and unfamiliar stimuli, especially those of caregivers over time as newborns adjust to postnatal life.

Taken all together, the current results add to the literature on neonatal behavioral preferences for sounds, and they complement previously published, electrophysiological results. They inform us about incipient capacities in very early development. In the absence of robust behavioral evidence for discrimination and therefore recognition of mother's voice, the presence of differential ERP responses confirms that the two auditory signals are, in fact, processed differently by the neonatal brain, at least for infants with uncomplicated pre- and early postnatal histories. Mother's voice is recognized, even in the impoverished form. The presence or, in this case, the absence of overt behavior informs about the relative salience of the signal and, perhaps, its hedonic valence, something that electrophysiological measures do not provide at present. Mother's familiar talking or reading voice with its unique and characteristic changes over time in pitch, rhythm and loudness is sufficiently salient and positively valenced that neonates are motivated to act to produce more of it (DeCasper and Fifer, 1980; Fifer and Moon, 1989). Without these characteristics in the sound of the maternal voice, they apparently are not.  It will be important to complement this understanding with increasing knowledge about the actions that result.

# Acknowledgments

We thank the families and staff at Madigan Army Medical Center for their support of the study. The views expressed in the article are those of the authors and do not reflect the official policy of the Department of the Army, the Department of Defense or the U.S. Government. This research was supported by a National Institutes of Health Grant HD 37954 to PKK. Software was generously provided by Anne Christophe of the Centre National de la Recherche Nationale, Paris, France. Student research assistants were: Jung Chong, R. Todd Davis, and Scott Nguy.

# References

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

Copyright © 2015 Moon, Zernzach and Kuhl. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution and reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Objective differentiation of neonatal EEG background grades using detrended fluctuation analysis

Vladimir Matic1, 2 \*, Perumpillichira Joseph Cherian3 , Ninah Koolen1, 2, Amir H. Ansari 1, 2 , Gunnar Naulaers 4 , Paul Govaert 5 , Sabine Van Huffel 1, 2, Maarten De Vos 6 and Sampsa Vanhatalo7

1 Department of Electrical Engineering (ESAT), STADIUS Centre for Dynamical Systems, Signal Processing and Data Analytics, KU Leuven, Leuven, Belgium, 2 iMinds Medical IT Department, Leuven, Belgium, 3 Section of Clinical Neurophysiology, Department of Neurology, Erasmus MC, University Medical Center, Rotterdam, The Netherlands, 4 Neonatal Intensive Care Unit, University Hospital Gasthuisberg, Leuven, Belgium, 5 Section of Neonatology, Department of Pediatrics, Erasmus MC-Sophia Children's Hospital, University Medical Center, Rotterdam, Netherlands, 6 Department of Engineering, Institute of Biomedical Engineering, University of Oxford, Oxford, UK, 7 Department of Children's Clinical Neurophysiology, HUS Medical Imaging Center and Children's Hospital, Helsinki University Central Hospital and University of Helsinki, Helsinki, Finland

#### Edited by:

Silvia Comani, Università degli Studi "G. d'Annunzio," Italy

#### Reviewed by:

Luc Berthouze, University of Sussex, UK Pål Gunnar Larsson, Oslo University Hospital, Norway

#### \*Correspondence:

Vladimir Matic, Division STADIUS-BIOMED, Department of Electrical Engineering, KU Leuven, Kasteelpark Arenberg 10 - bus 2446, 3000 Leuven, Belgium<EMAIL_ADDRESS><EMAIL_ADDRESS>
> Received: 25 November 2014 Accepted: 20 March 2015 Published: 23 April 2015

#### Citation:

Matic V, Cherian PJ, Koolen N, Ansari AH, Naulaers G, Govaert P, Van Huffel S, De Vos M and Vanhatalo S (2015) Objective differentiation of neonatal EEG background grades using detrended fluctuation analysis. Front. Hum.  We studied whether long range temporal correlations quantified by detrended fluctuation analysis (DFA) could be used in the neonatal EEG to distinguish different grades of abnormality in the background EEG activity. Long-term EEG records of 34 neonates were collected after perinatal asphyxia, and their background was scored in 1 h epochs (8 h in each neonate) as mild, moderate or severe. We applied DFA on 15 min long, non-overlapping EEG epochs (n = 1088) filtered from 3 to 8 Hz. Our formal feasibility study suggested that DFA exponent can be reliably assessed in only part of the EEG epochs, and in only relatively short time scales (10–60 s), while it becomes ambiguous if longer time scales are considered. This prompted further exploration whether paradigm used for quantifying multifractal DFA (MF-DFA) could be applied in a more efficient way, and whether metrics from MF-DFA paradigm could yield useful benchmark with existing clinical EEG gradings. Comparison of MF-DFA metrics showed a significant difference between three visually assessed background EEG grades. MF-DFA parameters were also significantly correlated to interburst intervals quantified with our previously developed automated detector. Finally, we piloted a monitoring application of MF-DFA metrics and showed their evolution during patient recovery from asphyxia. Our exploratory study showed that neonatal EEG can be quantified using multifractal metrics, which might offer a suitable parameter to quantify the grade of EEG background, or to monitor changes in brain state that take place during long-term brain monitoring.

 Scalp electroencephalography (EEG) is the most commonly used, non-invasive method in the continuous monitoring of brain function in the neonatal intensive care units (NICU). While neonatal EEG monitoring is becoming a standard of care in many NICUs, the interpretation of continuous EEG records in a 24/7 mode remains a global challenge, especially because of shortage of the special expertise required for an adequate visual reading (Boylan et al., 2010). One solution to this over the past decade has been the use of compressed displays, especially the amplitude integrated EEG (aEEG, a.k.a. CFM; De Vries and Hellström-Westas, 2005), which enables easier review of selected EEG features at bedside. However, it is now well recognized that aEEG trends are susceptible to artifacts that require special expertise in aEEG reading, and yet, the aEEG interpretation is qualitative and subjective.

The current main challenge in the EEG interpretation is to objectively and quantitatively characterize the spontaneous, ongoing brain activity, often called "background activity" in the EEG nomenclature. It has been shown that the EEG background is most informative when it comes to assessing acute state or predicting future outcome of the brain (Monod et al., 1972; Watanabe et al., 1999; Menache et al., 2002). Several background grading systems have been published over the years (Watanabe et al., 1999; Murray et al., 2009; Cherian et al., 2011; Walsh et al., 2011), and they combine visually (i.e., subjectively) observed EEG properties to yield one, holistic EEG grade. Perhaps due to the high interindividual variability and ambiguity in visual assessments, the clinically used EEG grading is remarkably rough. For instance, ranges of interburst intervals (IBI) between mild and moderate EEG grades go in 5 s increments from <5 to 5–10 s till >10 s, respectively (Murray et al., 2009; Cherian et al., 2011).

It is intuitively obvious, that so rough EEG grading cannot adequately reflect the time-varying physiological state of the brain, although it allows a technically straightforward translation of the criteria into automated classifiers. Intriguingly, classifiers based on these same criteria are not able to fully emulate the visual grading (Murray et al., 2009; Korotchikova et al., 2011; Stevenson et al., 2013), which indirectly calls for identification of novel EEG features with clinical relevance.

Recent work in basic neuroscience has provided ample evidence that many brain behaviors exhibit scale-free properties where dynamics of a given feature have no distinct spatial or temporal scale (Beggs and Plenz, 2003; Fransson et al., 2013; Iyer et al., 2014; Roberts et al., 2014). This was also recently shown to be the case for the EEG activity of full-term neonates that recover from perinatal asphyxia (Iyer et al., 2014; Roberts et al., 2014), one of the most common reasons for continuous EEG monitoring in the NICUs.

Scale-free dynamics in a complex system can give rise to selfsimilarity over temporal scales, i.e., long-range temporal correlations (LRTC), which may be assessed from the EEG using detrended fluctuation analysis (DFA; Peng et al., 1994; Hardstone et al., 2012). Recent work has shown that LRTC in brain function is significantly affected by various clinical conditions (Linkenkaer-Hansen et al., 2001; Parish et al., 2004; Stam et al., 2005; Monto et al., 2007). Compared to the conventional, visually defined EEG features, LRTC reaches beyond the visually perceived time scales, while it also integrates multiple time scales (from seconds to minutes) into one mathematically and conceptually transparent estimate. Moreover, DFA as a paradigm is very strongly supported by theoretical, experimental and clinical studies (for review, see Hardstone et al., 2012). In addition, a small number of neonates have been studied and long-range temporal correlations have been explored (Berthouze et al., 2010). Recent development of signal analysis methods has extended DFA to multifractal DFA (MF-DFA), which characterizes time series with multiple co-existent dynamic processes that may give rise to temporally local fluctuations including both extreme small and large magnitudes. Several studies have reported MF-DFA to offer additional insight in experimental and clinical neuroscience contexts (Kantelhardt et al., 2002; Zorick and Mandelkern, 2013). The neonatal EEG is, indeed, an example of signal with extreme, apparently erratic fluctuations in signal amplitudes, making it potentially suitable for MF-DFA paradigm as well. Notably, a recent study provided evidence of fractal behavior in time series of interburst intervals in preterm babies (Hartley et al., 2012).

The present study was set out to examine the possibility that neonatal EEG exhibits LRTC, which can be used as a feature to assess brain condition in a clinically relevant context. This entails answering the following questions: First, is it possible to assess reliable DFA exponents from the neonatal EEG at clinically relevant time scales from seconds to tens of minutes? Second, if DFA is suboptimal, can metrics from an existing MF-DFA paradigm reflect neonatal EEG in a meaningful way? Third, is it possible to use these measures to distinguish background grades of the neonatal EEG? The last question offers a benchmark to the existing analysis paradigms, while it is also directly relevant for development of novel features for automated background EEG classifiers.

# Materials and Methods

The present study was designed to proceed sequentially. We first examined the mathematical and physiological feasibility of the DFA approach in neonatal EEG analysis, and we studied the DFA fluctuation plots for differences in our clinical populations (different EEG grades). After observing mathematical limitations with the DFA paradigm, the work was continued to test MF-DFA paradigm, and finally to benchmark MF-DFA metrics with clinical EEG grading.

# Dataset

We used altogether 1088 15 min epochs (total 272 h) from 34 neonates that were recorded for clinical reasons at Sophia Children's Hospital, Erasmus MC (Rotterdam, the Netherlands). These neonates were part of a larger cohort of term asphyxiated neonates who were monitored from March 2003 till August 2007. Other results from this cohort have been published previously (Deburchgraeve et al., 2008; Cherian et al., 2011; De Vos et al., 2011; Matic et al., 2014). The clinical inclusion criteria for EEG monitoring were: gestational age of 37–43 weeks with clinical features of encephalopathy, and having at least one of the following features of birth asphyxia: (a) arterial pH of umbilical cord blood ≤ 7.1, (b) Apgar score ≤5 at 5 min, and (c) high clinical suspicion (like fetal distress, umbilical cord prolapse, difficult labor, or a history of convulsions). Exclusion criteria were: congenital cardiac abnormalities or other multiple anomalies, inborn errors of metabolism, and constantly isoelectric (<5µV) EEG activity. The study had the approval of the Erasmus MC Medical Ethical Review Board. Informed consent was obtained from the parents/guardians prior to the onset of the registration. The EEG recordings started 2–48 (median 19) h post-partum. The median duration of EEG was 31 (range 20–96) h. Seizures were recorded in 18 patients, however this clinical information was not considered in the later analyses further than confirming that distribution of DFA/MF-DFA metrics was not statistically significantly different between babies with vs. without seizures (calculated separately for each background category).

# EEG Sampling

The EEG was recorded with NicOne EEG recorder (Cardinal Healthcare, Madison, WI, USA) at sampling frequency of 256 Hz and our current data analysis was performed using bipolar montage with 12 derivations. Four epochs of 2 h continuous EEG recordings were selected per neonate with a minimal distance between them longer than 4 h.

# Feasibility of DFA in the Neonatal EEG

Detrended Fluctuation Analysis (DFA), originally introduced by Peng et al. (1994), is a widely applied method to quantify LRTC in various biological time series that exhibit self-similarity over multiple time scales. In the context of EEG analysis, DFA is typically computed by measuring amplitude dynamics of selected frequency bands (cf. Hardstone et al., 2012). In the first step, after band-specific amplitude extraction, the DFA procedure splits the signal into time windows (time scale—s), detrends the EEG segment, and computes the variance of amplitude fluctuations. In the second step, a DFA plot is generated by plotting the average signal variance per each window size on log-log coordinates (time scale—s, Fluctuation function—F in **Figures 1C,E**). Finally, a linear fit of this fluctuation function gives the DFA output, the scaling exponent α (cf. Appendix A; Hardstone et al., 2012).

Use of DFA paradigm requires defining several parameters so that they comply with the a priori known characteristics of the given biological signal. These include definitions of (1) the frequency band which is used for the signal envelope extraction (2) the epoch lengths, time scales, that are used to calculate the fluctuation function as well as (3) the range of time windows over which the slope of the linear fit (exponent α) is calculated (see **Figure 1C**).

# Epoch Length

We analyzed data epochs of 300, 600, 900, and 3600 s duration, but the final analysis was then carried out using 900 s (15 min) epochs. The choice was to strike a practical compromise between obtaining enough data for a reliable DFA plot, yet having likely a consistent background state in the given window. Our clinical background grading was performed initially in 15 min epochs, hence 900 s was practical for benchmarking purposes.

# Frequency Band

Neonatal EEG, especially in the acute state after asphyxic injury as in our data series, is dominated by intermittent bursting activity that consists of oscillations mostly at about 0.1–10 Hz. The lowest frequencies may be unpredictably contaminated by movement and other artifacts, so they were filtered out using a FIR highpass filter at 3 Hz. We chose to use a lowpass FIR filter at 8 Hz to maximally limit our analysis to the previously described 3–8 Hz band (Palmu et al., 2010; Tokariev et al., 2012; Omidvarnia et al., 2014). However, we also examined other frequency bands (1–3 and 1–5 Hz), and we got qualitatively comparable findings (data not shown).

# Time Scale (DFA Window)

Feasibility of the DFA method requires that the signal does exhibit genuine LRTC, which in an ideal case could be seen over all time scales up to the limits of the recording time. The shortest time scale is defined by the need to include multiple fluctuation cycles of interest to really assess its temporal correlations (cf. Hardstone et al., 2012). In the present context, the amplitude fluctuations in the early recovering neonatal EEG consist of intermittent bursting at multi-second scale. The selected frequency band (3–8 Hz) gives slowest fluctuations with a duration of around 1 s (Tokariev et al., 2012), which suggested to use the lowest time scale of 10 s. The upper range can be physiologically limited by the spontaneous changes in the brain state over time, especially vigilance state cycling that occurs 1–3 times per hour in the neonatal brain (Stevenson et al., 2014). However, the upper limit of feasible scales can also be reasoned mathematically, since estimation of the DFA exponent α is based on linear fitting in the DFA plot over the given time scale (see below; **Figure 1C**).

Next, we continued quantitatively by analysing how well the linear fitting suits to different ranges of windows (5–60, 10–60, 10–100, and 10–300 s) in the fluctuation functions. Two independent techniques were used: The first method was developed by Botcharova et al. (2013) based on the maximum likelihood approach that compares a set of alternative models: linear, polynomial, root, exponential, logarithmic, and spline functions. This paradigm identifies the best fit model when the involved number of parameters is minimal and over-fitting is penalized. The second method was less strict, and adopted from prior studies (e.g., Linkenkaer-Hansen et al., 2001), based on computing the standard linear regression fit, R 2. , and using a predefined threshold of 0.97 to indicate acceptance of the epochs for DFA analysis. More details about both results are shown in the Supplemental Data and they confirm that results obtained with the maximum likelihood based approach are valid.

# Feasibility Of MF-DFA in the Neonatal EEG

Various complex systems and physiological time series may show LRTC although they cannot be characterized with a single parameter (monofractal; for a review, see Ihlen, 2012). For instance, there may be clear monofractal dynamics at distinct, narrow time scales bounded by deflections in the log-log fluctuation plot (so called "crossover points"; Kantelhardt et al., 2002). Then, instead of a single linear approximation, the fluctuation function can be better described as a set of linear approximations, yielding different scaling exponents for different time scales. They are called multifractal time series, where the set of scaling exponents,

or parameters, are called q-order (generalized) Hurst exponents (Ihlen, 2012). For monofractal time series, the q-order Hurst exponents will have equal values, whereas in multifractal time

displayed in the X axis. The red line presents an example of linear fluctuation function that is here fitted for the time scale 10–60 s to obtain its slope, the DFA scaling exponent α. (D) DFA fluctuation functions are plotted for 15 min EEG

> series, the q-order Hurst exponents will be significantly diversified. A common reason for multifractality is that the signal may exhibit different LRTC for temporally local small and large

background EEG grades shows a significant difference between grades 2 and 3. (F) An example of the DFA fluctuation function (moderate EEG grade) where

clear existence of the crossover points are illustrated.

fluctuations. In contrast, monofractal time series are mainly noise like (e.g., white noise; cf. Ihlen, 2012) and they do not show local periods of either extremely small or large fluctuations. Hence, the known non-stationarity of the fluctuations in the neonatal EEG does, indeed, suggest that MF-DFA may offer better handling of the temporal EEG dynamics (e.g., local suppression and burst periods). Mathematically, MF-DFA differs from the DFA only in the calculation of the q-order fluctuation functions, Fq (**Figures 2B1,B2**), which is followed with the extraction of the qorder Hurst exponents (**Figures 2C1,C2**; see also Appendix A). Next, MF-DFA can be plotted as a so called multifractal spectrum, which can be quantified with mathematically robust metrics. The conventional MF-DFA metrics included mean hq, width hq, mean Dq, and height Dq (see **Figure 3D** for their graphical presentation). The width\_hq is calculated as the difference between maximal and minimal hq values, whereas height\_Dqis calculated as the difference between maximal and minimal Dqvalues (as suggested in Zorick and Mandelkern (2013). Mean hqand mean Dq are calculated by averaging all data points across spectra and across 12 channels. The statistical results are presented in **Figures 3E–H**.

In addition, the difference between DFA and MF-DFA can be observed only within the relatively short time scales (e.g., less than 100 s) and within MF-DFA concept the linear fitting from 1 s till 128 s was applied in this work.

# Benchmarking with EEG Grades

The ultimate goal in the present study was to establish the applicability of DFA or MF-DFA metrics in quantification of background EEG abnormalities. To this end, we had two benchmarking methods, visual and automated background classification.

Both DFA exponent and four MF-DFA metrics were calculated exactly the same way for each 15 min EEG epochs that also had a known visually scored background EEG grade inherited jointly for all four 15 min epochs within the given hour of EEG.

The visual background grading was performed by an experienced clinical neurophysiologist (PJC), blinded to the clinical details of the neonates. A three level EEG grading was adapted from Murray et al. (2009) as follows: (a) normal/mild (IBIs < 5 s; poorly organized sleep-wake cycle; recovered, continuous background EEG), (b) moderate (5 s ≤ IBIs < 10 s), and (c) severe (10 s ≤ IBI s< 60 s; dominated by prolonged and marked suppressions).

The automated background grading was taken from the automatically computed inter-burst intervals (IBI). This was generated using the automated IBI detection method described in Matic et al. (2012). The main steps of this algorithm are summarized in the Appendix B. Subsequently, we used mean IBI values and compared them with MF-DFA based parameters. All data points were extracted from 15 min cEEG epochs, thereby providing for every neonate 32 data values = 4 (epochs of 2 h) × 2 (h) × 4 (15 min segments in 1h cEEG).

# Statistical Analysis

Group comparisons between background grades were performed with the Kruskal-Wallis test, and we applied Bonferroni correction for multiple comparisons for the post-hoc tests. There were always three comparisons, hence the corrected significance level were p = 0.05/3 = 0.01 and p = 0.005/3 = 0.001. Correlations between MF-DFA metrics and IBI values were computed using Spearman correlation coefficient, and here we used more stringent significance level set to 0.01. Binomial statistics were used to evaluate whether significant findings in individual subjects can be considered statistically significant on the group level (see also Vanhatalo et al., 2004; Palva et al., 2005).

# Results

# DFA

# Visual Observations

DFA fluctuation function plots (an example shown in **Figure 1C**) were computed and examined first visually from hundreds of EEG segments. As a general observation, the ranges of window lengths with possible linear relationships were variable and limited. Our systematic qualitative comparison of window lengths ranging from 5–60 to 10–300 s (data not shown) suggested that linear fit should preferably be attempted at time scales of up to about 60 s. Then, we also reasoned that the lowest feasible limit should be around 10 s to avoid measuring trivial autocorrelations in the amplitude fluctuations (for further reasoning regarding the choice of minimum length, please see Hardstone et al., 2012). A direct comparison of hundreds of DFA fluctuation plots from EEG epochs with different EEG backgrounds showed that they all follow the same overall course, however, they also appear to differ in the segment corresponding to the shorter time windows (**Figure 1E**).

# Testing the Feasibility of DFA

In above, we had defined the epoch length to 900 s, which was a partly theoretical, partly practical choice. It is known that the DFA fluctuation plots provide more stable estimates of variance when longer epochs are analyzed. The practical reasoning was that the epoch length could be chosen as a multiple of our manual EEG grade classification that was set to 1 h. Obviously, a 900 s long, fixed epoch lengths cannot follow the more rapid dynamic changes known to take place in brain state of these babies (see also Berthouze and Farmer, 2012). Hence, other epoch lengths may be worth trying in future neonatal applications.

The actual feasibility of conventional DFA was tested by examining whether linear fit can adequately characterize the fluctuation plot within the tested range of time scales.

First, our analysis of 1088 epochs using the "maximum likelihood"-based technique suggested that only minority of them could be modeled with a linear fit. This finding was consistent even when we systematically varied frequency band (1–3, 1–5, and 3–8 Hz) or the range of fitting scale (5–60, 10–60, 10–100 and 10–300 s). The best set of parameters, which we fixed and consistently used in the study, was 3–8 Hz frequency band and 10–60 s window fitting size, and even then, only 200 out of 1088 15 min cEEG epochs (18%) were successfully modeled as linear. We also found that continuous background EEG (grade 1)

C1,C2 q-order Hurst exponent, Hq, is obtained as the linear fitting coefficient of the Fq fluctuation lines from B1,B2. Wider range of Hq values is

Longer tails and wider width between the tails in D2 compared to D1 are

associated with a higher degree of multifractality.

\*\*p < 0.01, \*\*\*p < 0.001.

all 1088 multifractal spectra lines is displayed for three background

is most likely to fit linear model (data not shown). The detailed results are presented in the Supplemental Figure S1.

Second, we also studied the more commonly applied measure, the linear fit R 2 using the previously adopted threshold of R 2 = 0.97. We found that about two thirds (68%; 731 out of 1088) of EEG epochs passed this test. For detailed comparison of R 2 values and maximum likelihood—based models, please see Supplemental Figure S4.

### Benchmarking with Clinical EEG Grades

For benchmarking with clinical EEG grades, we chose only those EEG epochs that had passed our linear fitting test (R 2 > 0.97), and we computed their DFA exponents (α). Group-wise comparison of α shows a wide variation within all EEG grades, while the mean α increases toward more severe EEG background. Notably, the normal/mild EEG grade (continuous EEG activity) was close to α = 0.5, **Figure 1D**, (an "uncorrelated" process; cf. Hardstone et al., 2012), while the most severe EEG background approached 1.0. Statistical group comparison of α shows that group 3 is significantly different from the other two groups (p < 0.01 after post-hoc correction; Kruskal-Wallis test).

# Crossover Points

Finally, we visually inspected a large number of DFA plots to observe possible systematic differences between EEG grades. **Figure 1E** shows an example where the plots have different overall forms: Those computed from the severe EEG grade epochs (black) are at the top, and the plots from the mild EEG grades (blue) are on the bottom of the group. It was also apparent, that the range of linear-appearing trends varied between EEG grades, and that they were often clearly multiphasic: short time scales with steep rise, middle window lengths with flatter rise, and again steep rise at the longest time scales. An example of an EEG epoch with very clear "crossover points," i.e., change of slope in the fluctuation function (Hardstone et al., 2012) is shown in **Figure 1F**. Finding such crossover points, suggests that the signal might exhibit multifractal behavior (Kantelhardt et al., 2002; Ihlen, 2012; however, see Hu et al., 2001). Describing such signal with one exponent, as in the conventional monofractal DFA, will overlook key dynamic properties, and the results will be particularly sensitive to chance factors in the conventional DFA analysis.

In our dataset, the crossover point in the burst-suppression background (red and black DFA fluctuation functions in **Figure 1E**) was usually encountered within 30–60 s time scale [1.5–2 in log10(time) s]. Continuous EEG background (blue DFA fluctuation functions) were successfully modeled as linear for longer window sizes [up to 100 s or 2 in log10(time) s], and the crossover points mainly appeared in the interval of 100–300 s [2– 2.5 in log10(time) s]. In comparing different EEG grades, the different crossover points hence preclude use of single, monofractal DFA exponent, and supports exploration of MF-DFA metrics.

# Multifractal DFA (MF-DFA)

We next studied the feasibility of MF-DFA in neonatal EEG by applying a previously published method (Kantelhardt et al., 2002; Ihlen, 2012). The special clinical attraction in MF-DFA comes from its clear metrics that have been recently shown to aid discrimination of brain states (Zorick and Mandelkern, 2013) and pathologies (Zheng et al., 2005).

# Feasibility of MF-DFA

The first step here was to show that neonatal EEG signal is characterized by multiple rather than single exponents. The presence of multifractal behavior in at least some neonatal EEG data becomes apparent when comparing the temporal behavior of two common neonatal EEG patterns in asphyxiated infants, the normal continuous EEG signal and the severely abnormal EEG with marked intermittency called burst-suppression (see **Figures 2A1,A2**). For those signals, we calculated the q-order fluctuation functions Fq (**Figures 2B1,B2**) with q values from −5 till 5, equidistantly sampled with a 0.5 unit-step as suggested by Ihlen (2012), notably, different range could be examined as well. The q factor is able to detect mixed occurrence of large fluctuations with periods of negligible fluctuations, which in our case characterizes the burst-suppression background. In contrast, the continuous EEG consists of neither very large nor very small fluctuations (**Figure 2A1**).

Comparison of q-order fluctuation functions (Fq) taken from these two different EEG patterns shows their fundamental difference. In the normal, continuous EEG activity, the slopes of the fluctuation functions appear comparable to each other. Hence, they could also be reasonably described with only one DFA fluctuation function (**Figure 2B1**; q = 2), and be considered as monofractal process. In contrast to this, EEG periods with burstsuppression are much more variable, and different time epochs exhibit a notably wide spread of Fq (**Figure 2B2**), suggestive of an underlying multifractal process (see also Ihlen, 2012). Taken together, these observations suggest that measuring LRTC in neonatal EEG datasets like ours should require quantification of multifractal dynamics.

Subsequently, we calculated q-order Hurst exponents, Hq (Ihlen, 2012), which returns linear approximation of the Fq fluctuation functions (**Figures 2C1,C2**). This again discloses a clear difference between background EEG activity grades. In the continuous EEG, the relation between Hq and q factor shows a monotonic decrease (characteristic of a monofractal process; Ihlen, 2012), whereas the Hq plot of burst-suppression EEG (**Figure 2C2**) shows a relatively sharp and nonlinear decline (characteristic of a multifractal process). These q-order Hurst exponent functions can be readily transformed into a multifractal spectrum (**Figures 2D1,D2**), which outlines the relation of local temporal changes of the Hq and the distribution of its probability density function (see Ihlen, 2012). Multifractal spectrum represents two measurable dimensions, the Dq (q-order singularity/fractal dimension) and the hq (q-order singularity exponent). In the visual assessment of multifractal spectra, differences can be noticed in a horizontal and vertical positioning (hq, Dq values), a width (width\_hq), as well as in the general shape of the multifractal spectra reflecting temporal variations of local Hurst exponents.

The MF-DFA spectra were then computed for each (n = 1088) of the 15 min long EEG epochs. Example of these spectra are plotted from each EEG background grade (**Figures 3A–C**). Visual inspection of the spectra shows that worsening of EEG background (from mild to severe) is associated with a shift to the right (higher hq value). Grand averages of the spectra are shown in **Figure 3D** (1088 epochs from 12 channels: total 13056 spectra).

Other differences among EEGs were also observed visually within the tails of the spectra: the left tails, corresponding to positive q values, are longer in the EEG with local large fluctuations, i.e., the severe EEG grade with IBIs > 10 s (**Figure 3C**). Moderate and milder EEG background traces do not have such pronounced structural differences and the corresponding tails are respectively shorter (**Figures 3A,B**). On the other hand, the right tails, reflect the amplitude fluctuations within local periods of low activity (e.g., suppressed EEG periods).

# Benchmarking with Clinical EEG Grades

Comparison of MF-DFA metrics between EEG grades showed highly significant differences between the three EEG grades. The mean height of spectra, mean hq, was significantly (p < 0.01; Kruskall-Wallis test) increased with EEG severity (mild 0.63 ± 0.10; moderate 0.67 ± 0.09; severe 0.70 ± 0.09). The width of the spectra, width hq, was also significantly (p < 0.001) increased with increasing severity of EEG background (mild 0.56 ± 0.13; moderate 0.67 ± 0.16; severe 0.73 ± 0.19). At the same time, the mean Dq decreased significantly (p < 0.01) with increasing severity of EEG background (mild 0.75 ± 0.05; moderate 0.72 ± 0.05; severe 0.70 ± 0.06). Height of the MF-DFA spectra, height\_Dq, was significantly (p < 0.01) different between moderate and severe EEG grades (moderate 0.70 ± 0.09; severe 0.76 ± 0.11). These findings together suggest that MF-DFA is superior to the conventional DFA in distinguishing between neonatal EEG grades.

A further inspection of the shapes of MF-DFA spectra shows additional qualitative spectral features that are different between background grades. Namely, the mild/normal EEG grade does typically correspond to asymmetric spectra with a longer right side tail, whereas the spectra of severe EEG grade are more asymmetric with longer left side tail. This asymmetry is not measured by the currently used, established MF-DFA metrics. It may, however, readily explain the lesser differences between groups seen in mean\_Dq and height\_Dq, which are insensitive to asymmetry. To examine this further in a quantitative manner, we also computed the metrics from the right side tail only (see below).

# Benchmarking and Application into Intraindividual Context

While the group level findings above are by themselves strong, clinical applicability of any novel study paradigm requires performance at individual level as well (see discussion in Ahtola et al., 2014). In order to preliminarily probe this issue, we studied those 34 neonates in our dataset that had four EEG epochs of 2 h duration randomly selected at time points with at least 4 h intervals. We then computed MF-DFA metrics from 15 min segments within each 2 h epochs (total 32 epochs in each infant). The long intervals between epochs yielded us datasets where EEG grade, or at least the mean IBI, was changing due to variation in the clinical state of the patient. As a benchmark to MF-DFA metrics, we calculated the mean IBI for each 15 min epoch using our previously described automated detector (Matic et al., 2012). Nine patients were excluded because their range of IBIs was too narrow to reasonably compute correlation coefficients with another parameter, i.e., IBI range was only within 3 to 6 s interval. This left 25 babies for our analysis. Finally, non-parametric, ranking correlation (Spearman) was computed between mean IBI levels and MF-DFA spectra metrics in each patient.

# **Measures of the full MF-DFA spectra**

Individual results are shown in Supplementary Table 1, Supplemental Data. We observed that in 10 out of 25 babies there was a significant (p < 0.01) correlation at the individual level, which finding is statistically very significant at the group level (p < 1e–10; Binomial statistics). As shown in **Figure 4**, for an example infant (n = 21), a visually apparent temporal relationship was observed between decreasing IBI trend and changes in at least three of the MF-DFA parameters (mean\_hq, width\_hq, and mean\_Dq).

# **Measures of the right side tail of MF-DFA spectra**

Closer inspection showed that the signs of correlation coefficients (ρ) were inconsistent even for infants with statistically significant correlations (Supplementary Table 1, Supplemental Data). To search for reasons we inspected further the shapes of MF-DFA spectra, and observed their asymmetry, especially in the severe EEG grades. For instance, visual comparison of individual MF-DFA spectra shown in **Figure 4B** from different time epochs during recovery shows that the shape of spectra changes over time, most notably as a shift from longer left tail to longer right tail as well as the shift of the right tail to the left when going from severe to normal EEG background, respectively. Since the current MF-DFA metrics are insensitive to asymmetry of the multifractal spectra, it is possible that change in symmetry with EEG grades would confound the correlations. Hence, we investigated the parameterization of the right tail of the MF-DFA spectra, and we found these metrics to exhibit much stronger relationship to the mean IBI values. In particular, for 13 out of 25 patients, we found that at least 3 metrics namely width hq (ρ = −0.70±0.15), mean Dq (ρ = 0.68 ± 0.13) and height Dq (ρ = −0.67 ± 0.11), were significantly correlated with the mean IBI values, showing consistent range of correlation coefficient values (ρ) as detailed at subject level in Supplementary Table 2, Supplemental Data.

Additionally, we illustrate how the MF-DFA parameters (extracted from the right tail of the spectra) are correlated with IBI values for five infants that showed higher Spearman correlation p-values (selected from the Supplementary Table 2, Supplemental Data) (**Figure 5**).

# Discussion

Our study shows that changes in long range temporal dynamics in the early neonatal EEG may be characterized with the DFA paradigm, however the conventional DFA is compromised and the recently introduced multifractal DFA can better disclose

epochs].

neonatal EEG dynamics in a way that correlates with the current clinical EEG grading systems. Our clinical benchmarking of DFA feasibility extends prior literature to include full-term neonatal EEG as a potential application of DFA/MF-DFA paradigms. The present findings corroborate earlier reports that DFA or MF-DFA metrics can quantify brain dynamics associated with prematurity (Hartley et al., 2012) and neurological conditions (Linkenkaer-Hansen et al., 2001; Parish et al., 2004; Monto et al., 2007). We also suggest a novel MF-DFA metric for EEG analysis based on observing asymmetric change in MF-DFA spectra between clinically relevant brain states.

three states quantifying 15 min of burst-suppression (black-recorded at 6 h

The physiological underpinnings of multifractal behavior in the neonatal brain remain elusive, and our dataset is suboptimal for a tour de force exploration of the phenomenon of multifractality in the neonatal brain activity. However, our observations provide some tentative ideas to how MF-DFA may be suitable in such data: The well-known hallmark of both the early preterm EEG and the severely abnormal EEG of a term newborn is a salient intermittency (Vanhatalo and Kaila, 2006; André et al., 2010): seemingly erratic ruptures of transient cortical activity alternate with periods of relative quiescence. Considering the structural immaturity of underlying brain networks (Kostovic´ and Judaš, 2010) and/or the pathological mechanisms related to resource constraints (Roberts et al., 2014), it is conceivable that they would readily compromise the extreme long range temporal correlations previously shown in the adult brain (Linkenkaer-Hansen et al., 2001). There is ample evidence from the neonatal EEG studies that cortical activity may continue even after extensive structural lesions that lead to physical disconnection of cortical areas. The multifractality could hence reflect a spatial and/or temporal fragmentation of large scale brain activity, resulting in multiple partly independent mechanisms in brain activity. Nevertheless, more rigorous confirmation with more extensive datasets and stringent mathematical testing of alternative models will be necessary to confirm fundamental behaviors of this kind (for instance, see Clauset et al., 2009; Roberts et al., 2014).

In addition to showing the significant discrimination of background grades with the established MF-DFA metrics, our closer visual inspection of MF-DFA spectra suggested further differences. The highest difference between grades was achieved for the relatively small window size values [<30 s or 1.5 in log10(time)s]

(**Figure 2B2**), which is most sensitive to temporally local variations in signal fluctuations. We also observed that the different background grades yield saliently different shapes, especially the asymmetry, in their corresponding MF-DFA spectra. This is not quantified by the established MF-DFA metrics (cf. Zorick and Mandelkern, 2013), hence we proposed quantification of the right side tail of the spectra for an improved sensitivity. These measures were, indeed, significantly correlating with IBIs in a large proportion of infants. It is possible that the prominent intermittency of the burst suppression signal has an important effect on the MF-DFA metrics. Conceivably, burst suppression EEG consists of two overlaid, but at least partly independent physiological processes: intermittent bursting behavior and the varying levels of spontaneous activity between the bursts, both with unique scaling properties (Roberts et al., 2014). A similar situation is also found during the early intermittent EEG activity of preterm babies where two mechanistically different brain activities are superimposed on each other (Vanhatalo and Kaila, 2006), giving rise to unique long range temporal correlations in the inter-event time series (Hartley et al., 2012). The physiological meaning of asymmetry remains to be established in simulations and experimental studies, however our data suggests that it relates to continuity and perhaps the relative shapes of bursts that interrupt the interburst silence.

In search of clinical relevance, we benchmarked the MF-DFA metrics with both visual and automated background EEG grading, and show that all our four metrics (mean hq, width hq, mean Dq, height Dq) do significantly co-vary with the visually scored EEG background grades. Three (mean hq, width hq, mean Dq) of our metrics were significantly different between all EEG groups. In addition, the application within patients showed that MF-DFA metrics might even reflect temporal change of background EEG states at individual level as suggested by the significant correlations with IBI levels in 10 out of 25 infants. These correlations suggest that MF-DFA metrics and IBI would at least partly reflect similar pathological brain dynamics, or at least the amplitude fluctuations of burst suppression are reflected in the MF-DFA spectra. However, despite the general correspondence with IBI levels and its evolution, the match between IBI and MF-DFA metrics is neither fully consistent nor complete. Such partial mismatch are consistent with the idea that MF-DFA metrics reflect some as yet latent EEG signal dynamics.

Beyond the physiological interest, our present findings have potential implications in future construction of computational classifiers for neonatal EEG. A number of studies have reported development of various automated EEG background classifiers. Some of them are designed to directly replicate the visual criteria. Using this approach, the recent study by Korotchikova et al., 2011, combined a larger set of quantified EEG features and found them indiscriminative with respect to background EEG classes, with special challenges in distinguishing mild vs. moderate grades. A more recent attempt along this line was able to discriminate normal vs. abnormal EEG backgrounds as a dichotomous choice (Stevenson et al., 2013), which does not yet support its use in EEG monitoring or for more refined needs. Another approach has been to optimize a learning algorithm with a smaller dataset to detect extreme classes of EEG background, such as severe burst-suppression patterns (Löfhede et al., 2010), which is further used to output IBI detections. The design of such classifier will ideally allow accurate quantitation of IBIs (Flisberg et al., 2011), however, it will not allow constructing methods that can monitor evolution of EEG activity over background grade boundaries, which obviously happens among a majority of clinical patients during their brain monitoring.

Prior extensive engineering works have mostly focused on identifying discrete EEG states (e.g., burst suppression), or on discrimination of normal from abnormal backgrounds as benchmarked with visual criteria. One of the shortcomings in benchmarking with visual EEG grading is its temporal coarseness. Visual grades are commonly assigned hourly for slow changes such as sleep-wake cycling (Murray et al., 2009; Matic et al., 2014; Stevenson et al., 2014). Other EEG grading systems have been developed to include even longer term evolution, up to over 24 h, to characterize recovery from hypoxic insults (Cherian et al., 2011), which as a principle reminds of the conventional approach used in the aEEG assessment (Hellström-Westas and Rosén, 2006). We are not aware of objective computational, graded measures of background EEG activity although the clinical experience shows that the key part of brain monitoring is the follow-up of patients through different stages over hours. Our findings suggest that measures based on MF-DFA metrics may offer an advantage here: for instance, these metrics may distinguish moderately abnormal EEG from the severe and normal backgrounds, which is clinically very important. Brain states change more rapidly in the real life, hence we wanted to compare MF-DFA-based metrics also with the automated algorithm that detects IBI levels in 15 min epochs. This showed that MF-DFA metrics do co-vary with IBI, and they do evolve over time. A systematic characterization and benchmarking with existing trend paradigms, especially aEEG (Hellström-Westas and Rosén, 2006) will help to identify the clinical added value of MF-DFA measures.

We consider our study as exploratory, providing only proof of concept. Further prospective and larger studies are needed to determine the clinical robustness and added value obtained from the proposed method, as well as to validate its use in the prospective clinical context. The first practical consideration is the age range of infants. Our dataset was limited to full-term neonates, but we do not see any a priori physiological or technical reason why the same temporal dynamics could not be clinically meaningful in the preterm babies as well. The second practical consideration is the unknown relative sensitivity of our metrics to artifacts that are unavoidably present in real life recordings from intensive care units. This can be assessed formally by simulated addition of real EEG artifacts (cf. Matic et al., 2009; De Vos et al., 2011; Räsänen et al., 2013). The future artifact avoidance strategies may include rejection by amplitude criteria (cf. Palmu et al., 2013), although missing data segments may also introduce confounders. The third practical consideration is the added clinical value relative to existing measures, especially the aEEG paradigm. It is conceivable that our metrics provide a theoretically and functionally complementary view on neonatal brain dynamics, yet the actual added bedside value needs to be determined in future clinical studies. There are also technical considerations for future studies. First, the known alternative ways in the MF-DFA analysis, such as wavelet transform modus maxima (Muzy et al., 1993), might be useful to further improve the robustness and consistency (Zorick and Mandelkern, 2013). Second, our study demonstrate as yet uncharacterized changes in the spectral shape, which

References

suggests that further theoretical and practical studies are needed to define spectral measures that optimally reflect changes in the neurophysiological and/or clinical state of the brain.

# Acknowledgments

We thank our patients and their parents for their cooperation and the neurotechnologists who made the high quality EEG recordings.VM, AA, NK, and SV are supported by: Research Council KUL: CoE PFV/10/002 (OPTEC); PhD/Postdoc grants Flemish Government: FWO: projects: G.0427.10N (Integrated EEG-fMRI), G.0108.11 (Compressed Sensing) G.0869.12N (Tumor imaging) G.0A5513N (Deep brain stimulation); PhD/Postdoc grants IWT: projects: TBM 080658- MRI (EEG-fMRI), TBM 110697-NeoGuard; PhD/Postdoc grants. iMinds Medical Information Technologies SBO 2015, ICON: NXT\_Sleep. Belgian Federal Science Policy Office: IUAP P7/19/(DYSCO, "Dynamical systems, control and optimization," 2012–2017). Belgian Foreign Affairs-Development Cooperation: VLIR UOS programs. EU: The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007–2013)/ERC Advanced Grant: BIOTENSORS (n◦ 339804). This paper reflects only the authors' views and the Union is not liable for any use that may be made of the contained information. Other EU funding: RECAP 209G within INTERREG IVB NWE program, EU MC ITN TRANS-ACT 2012 (n◦ 316679), ERASMUS EQR: Community service engineer (n◦ 539642-LLP-1-2013) SV was supported by the Academy of Finland grant 253130, Finnish Cultural Foundation, Päivikki ja Sakari Sohlberg Foundation, and Juselius Foundation.

# Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fnhum. 2015.00189/abstract

Figure B1 | Steps of the automated IBI detection method are

illustrated. First, EEG signal is segmented using an adaptive segmentation algorithm (A1) Next, the segments are classified into the low amplitude class (approximately <10 µV) and depicted in blue. (A2) The temporal profile signal counts the prevalence of low amplitude segments and it is used to detect IBI patterns (A3). The time that the temporal profile signal exceeds the threshold determines the IBI duration and the IBI is detected.

glossary. Neurophysiol. Clin. 40, 59–124. doi: 10.1016/j.neucli.2010. 02.002

Beta bands across a wide age range. Clin. Neurophysiol. 121, 1187–1197 doi: 10.1016/j.clinph.2010.02.163

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

Copyright © 2015 Matic, Cherian, Koolen, Ansari, Naulaers, Govaert, Van Huffel, De Vos and Vanhatalo. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Interhemispheric synchrony in the neonatal EEG revisited: activation synchrony index as a promising classifier

# *Ninah Koolen1,2\*†, Anneleen Dereymaeker 3†, Okko Räsänen4, Katrien Jansen3, Jan Vervisch3, Vladimir Matic 1,2, Maarten De Vos 5,6, Sabine Van Huffel 1,2, Gunnar Naulaers 3 and Sampsa Vanhatalo7*

#### *Edited by:*

Silvia Comani, Università degli Studi "G. d'Annunzio," Italy

#### *Reviewed by:*

Simon Farmer, National Hospital for Neurology and Neurosurgery and University College London, UK Andrea Guzzetta, University of Pisa, Italy

#### *\*Correspondence:*

Ninah Koolen, Division STADIUS-BIOMED, Department of Electrical Engineering, University of Leuven, Kasteelpark Arenberg 10 - Bus 2446, 3000 Leuven, Belgium e-mail: ninah.koolen@ esat.kuleuven.be

†These authors are joint first authors.

A key feature of normal neonatal EEG at term age is interhemispheric synchrony (IHS), which refers to the temporal co-incidence of bursting across hemispheres during trace alternant EEG activity. The assessment of IHS in both clinical and scientific work relies on visual, qualitative EEG assessment without clearly quantifiable definitions. A quantitative measure, activation synchrony index (ASI), was recently shown to perform well as compared to visual assessments. The present study was set out to test whether IHS is stable enough for clinical use, and whether it could be an objective feature of EEG normality. We analyzed 31 neonatal EEG recordings that had been clinically classified as normal (n = 14) or abnormal (n = 17) using holistic, conventional visual criteria including amplitude, focal differences, qualitative synchrony, and focal abnormalities. We selected 20-min epochs of discontinuous background pattern. ASI values were computed separately for different channel pair combinations and window lengths to define them for the optimal ASI intraindividual stability. Finally, ROC curves were computed to find trade-offs related to compromised data lengths, a common challenge in neonatal EEG studies. Using the average of four consecutive 2.5-min epochs in the centro-occipital bipolar derivations gave ASI estimates that very accurately distinguished babies clinically classified as normal vs. abnormal. It was even possible to draw a cut-off limit (ASI∼3.6) which correctly classified the EEGs in 97% of all cases. Finally, we showed that compromising the length of EEG segments from 20 to 5 min leads to increased variability in ASI-based classification. Our findings support the prior literature that IHS is an important feature of normal neonatal brain function. We show that ASI may provide diagnostic value even at individual level, which strongly supports its use in prospective clinical studies on neonatal EEG as well as in the feature set of upcoming EEG classifiers.

**Keywords: interhemispheric synchrony, biomarker, preterm infant, brain monitoring, neonatal EEG**

# **INTRODUCTION**

Recent progress in basic neuroscience, neuroimaging, and neonatal care has raised interest in the understanding of physiological and pathological processes in the preterm and neonatal brain. Electroencephalography (EEG) is a non-invasive and sensitive tool for evaluating brain function in the neonatal period. A key component in early brain functional development is the emergence of functional networks, the long range connectivity between and within brain hemispheres, which lay the basis for the development of neurocognitive capabilities (Uhlhaas et al., 2010; Lubsen et al., 2011; Omidvarnia et al., 2014). In this context, it is intriguing that synchrony between hemispheres has been considered as an important feature of normal neonatal brain function since the seminal studies in the late 1970s (Lombroso, 1979).

Clinical classification of spontaneous neonatal EEG is traditionally based on visual assessment of multiple features: continuity, quality of sleep wake organization, interhemispheric synchrony (IHS), symmetry and amplitude. However, visual EEG interpretation requires expertise, and there are no objective standards for classification schemes, nor any other established method to yield appropriate diagnostic accuracy to support modern neuroscience or evidence based medicine. This has greatly compromised comparisons of results between different studies. A quantitative analysis of EEG activity, including automated analysis of selected features of cortical function, could create

1 Division STADIUS, Department of Electrical Engineering (ESAT), University of Leuven, Leuven, Belgium

2 iMinds-KU Leuven Medical IT Department, Leuven, Belgium

3 Department of Development and Regeneration, University of Leuven, Leuven, Belgium

an objective and appropriate classification scheme for neonatal EEG. Moreover, it could also support and accelerate visual EEG analyses in medical centers without access to highly skilled EEG interpretation.

Visual estimation of IHS is based on observation of coincidence of bursts in spontaneous background EEG activity during tracé alternant/discontinue (Holmes and Lombroso, 1993). IHS is interpreted as a sign of connectivity or functional interaction between hemispheres and is hence considered an important feature of normal brain function, resulting from the development of callosal connections (Kostovic and Jovanov-Milosevic, 2006; Dudink et al., 2008).

Despite the long and widespread clinical use of IHS, there is no physiologically plausible definition of IHS for visual EEG reading (for more details, please see Räsänen et al., 2013). A quantitative measure for IHS, the activation synchrony index (ASI), was recently developed by Räsänen et al. (2013). The ASI is based on statistically measuring the temporal delay between two signal energies, and was shown to correlate with visually rated IHS grades. It was also shown to clearly outperform other methods proposed for the same purpose in the recent literature (Räsänen et al., 2013).

Our ultimate aim is to use ASI more widely as a clinical biomarker or as a feature in EEG classifiers, which poses special requirements for individual stability. In addition, the analysis settings such as the epoch lengths need to be practical with respect to clinical reality, where obtaining longer high quality epochs is often challenging due to intermittent, trivial artifacts. The aim of this study was to test the ability of the ASI to distinguish normal and abnormal neonatal EEG recordings. Therefore, we first defined the optimal ASI parameters to maximize its intraindividual and technical stability, and then tested ASI as an input for a classifier to identify abnormal EEG activity.

# **MATERIALS AND METHODS**

The present study consisted of two discrete phases (see **Figure 1**). First, we optimized the intraindividual stability of ASI. The earlier development of ASI had aimed to optimize parameter settings for maximal distinction between asynchrony grades as well as robustness against artifacts. Our present study continued ASI optimization by searching for EEG epoch lengths and channel derivations that would strike a balance between technical reliability and clinical practice. For clinical practice, an important consideration is that the algorithm should work equally well on EEG samples either with longer continuous epochs or with multiple shorter epochs. In addition, we tested our optimized ASI protocol on a new EEG dataset that had been classified using conventional clinical criteria, with the aim of seeing how well ASI alone is able to distinguish normal and abnormal EEG recordings. This latter aim is also a rigorous test of the conventional, albeit quantitatively untested assumption, that interhemispheric asynchrony would be an essential property in the abnormal EEG near term age.

### **SUBJECTS AND EEG RECORDINGS**

A total of 31 neonatal EEG traces [postmenstrual age (PMA) 36–42 weeks], were retrospectively selected by an expert (A.D.) to assess the spontaneous background EEG activity. These EEG recordings were collected from a larger EEG dataset, recorded for clinical purposes. All EEG data were recorded at 256 Hz with a video-synchronized EEG device (BRAIN RT, OSG equipment, Mechelen, Belgium). After skin preparation (Nuprep Gel), 10–17 Ag/AgCl cup electrodes were placed according to the international 10–20 standard locations. The reference electrode was Cz. Electrode impedance was below 10 k at the start, and the signal quality was monitored visually throughout the recording. The minimum recording time was 4 h to record multiple vigilance states. The protocol was reviewed and approved by the relevant Ethics Committee of the University Hospitals of Leuven, Belgium.

The EEG records were reviewed for clinical purposes by two independent raters (A.D. and K.J.) who were not aware of the later use of the EEG for this study. Hence, they were fully blinded to both the overall design as well as the numerical ASI results. The clinical assessment is based on holistic, conventional visual criteria (see Table 1 in Supplementary Material). If there was any difference in the interpretation of an EEG, a consensus was reached after re-evaluation (2 patients, Table 1 in Supplementary Material). The following parameters were assessed on the whole measurement, according to standard developmental features in preterm and term EEG (Scher et al., 2005; André et al., 2010; Shellhaas et al., 2011; Hayashi-Kurahashi et al., 2012): (1) brain activity cycling as being normal or abnormal for the age; (2) grade of continuity; (3) age-specific landmarks: amplitude, disorganization, and dysmaturation patterns; and (4) qualitative synchrony in quiet sleep. The overall clinical classification of the EEG recording as normal (*n* = 14) vs. abnormal (*n* = 17) was based on observing at least two features that were not appropriate for age.

# **ACTIVATION SYNCHRONY INDEX**