Greater Repertoire and Temporal Variability of Cross-Frequency Coupling (CFC) Modes in Resting-State Neuromagnetic Recordings among Children with Reading Difficulties

Introduction

The neurophysiological basis and functional significance of spontaneous brain activity remains disputable among neuroscientists, ever since the initial demonstration of focally reduced activity associated with increased cognitive demands (Pinneo, 1966). More recently, clinical applications of spontaneous connectivity patterns have been explored in search of indices of abnormal brain function using a variety of imaging methods, including structural MRI, diffusion MRI, and functional MRI, electroencephalography (EEG) and magnetoencephalography (MEG; for reviews see Bassett and Bullmore, 2009; Sporns, 2014; Braun et al., 2015). Resting-state functional connectivity has been studied in relation to disorders of consciousness (Heine et al., 2012), Alzheimer's disease (Buckner et al., 2008; Stam et al., 2009; Khazaee et al., 2015), schizophrenia (Bassett et al., 2012; Damaraju et al., 2014), depression (Jin et al., 2011; for reviews see Rosazza and Minati, 2011; Dutta et al., 2014), in mTBI (Dimitriadis et al., 2015d; Antonakakis et al., 2016), in dyslexia (Koyama et al., 2010, 2011; Dimitriadis et al., 2013b), and in normal populations (Raichle et al., 2001; Greicius et al., 2003; Fair et al., 2008; De Pasquale et al., 2010; Brookes et al., 2011).

In this line of research, the brain is viewed as a complex network comprised of sets of smaller subsystems, interacting in a dynamic manner, both locally and at longer ranges (Tognoli and Kelso, 2014a,b). Local neuronal assemblies specialized for a particular neuronal process are tuned to generate oscillations in specific frequency bands, and that the interaction of a given assembly with more distant neuronal assemblies is “carried” by signals in the prominent oscillating frequency of the source population (Shew et al., 2009). This cross-frequency coupling (CFC) can support the encoding of nested hierarchical relations, which is crucial for the representation of composite objects, likely following specific syntactic rules known to both sender and receiver, and enables interpretation of long uninterrupted signals by the latter (Buzsáki, 2010). The same principle may apply to the communication between neuronal assemblies located within the same cortical region. This mechanism promotes accurate timing between different oscillatory rhythms and dynamic control of distributed functional networks (Varela et al., 2001; Buzsáki, 2006; Canolty and Knight, 2010). There is rapidly accumulating experimental evidence supporting this role of CFC in cognition (Jensen and Colgin, 2007; Canolty and Knight, 2010; Palva and Palva, 2011; Buzsáki and Watson, 2012; Jirsa and Müller, 2013; Dimitriadis et al., 2015a,c; Dimitriadis et al., 2016).

Findings pooled from several studies indicate that temporal parsing of neuronal activity in different frequency ranges is very well-conserved in mammalian evolution (Northcutt, 1995). This preservation suggests that temporal coordination of distributed brain processes, as reflected by oscillatory patterning, synchronization, phase locking, and cross-frequency coupling, might have important functions and not be epiphenomenal (Buzsáki et al., 2013). Such fine-tuned processes can be disrupted by a number of developmental variations including: (i) Changes in the subunit composition of ligand- and voltage-gated membrane channels altering the time constants and resonance properties of neurons and microcircuits; (ii) alterations in white matter integrity such as myelination and axonal diameter affecting conduction velocity and overall efficiency of neuronal signaling. In all of these variations one expects to find alterations in variables reflecting brain dynamics, such as the power of particular oscillations and the extent and precision of synchronization in the various frequency ranges and their cross-frequency relationships (Buzsáki and Watson, 2012).

In this context, appropriate descriptors of brain dynamics derived from EEG and MEG recordings may help characterize the pathophysiological substrate of psychiatric and developmental disorders in the form of “oscillopathies” or “dysrhythmias.” With respect to reading and associated developmental difficulties, there is accumulating evidence that the degree of myelination in left hemisphere cortico-cortical tracts correlates positively with reading skill (Niogi and McCandliss, 2006; Hoeft et al., 2011). In a recent fMRI study (Richards et al., 2015) reported strong associations between measures of white matter integrity assessed through Diffusion Tensor Imaging and functional connectivity using fMRI in several brain regions in children varying in reading and spelling ability. Interestingly, degree of connectivity was both down- and up-regulated among struggling readers in complex ways.

Distinct forms of cross-frequency interactions have been described: (i) power to power, (ii) phase to phase, (iii) phase to frequency and, (iv) phase to power (Jensen and Colgin, 2007). The latter type, also referred to as phase-amplitude modulation, has been repeatedly described in both animal and human data (Osipova et al., 2008; Cohen et al., 2009a,b; Colgin et al., 2009; Axmacher et al., 2010a,b; Tort et al., 2010) and is perhaps the most suitable CFC type to support high-precision timing of signaling between neuronal assemblies located either within the same or in remote cortical regions. This form of CFC can be estimated through the phase-to-amplitude coupling index (PAC) which quantifies the modulation of the amplitude of higher frequencies by the phase of slower rhythms. PAC captures the phenomenon whereby a specific number of cycles of (higher) frequency rhythms in the target population can be oscillated within the phase of slower rhythms of the source population. This number is presumed to reflect the amount of information that can be exchanged between the two neuronal oscillators operating at different characteristic frequencies. The PAC index quantifies both the strength of the interaction and the prominent frequencies displayed by each interacting pair of oscillators at each point in time (i.e., time window of recorded activity). In the present study PAC serves as an index of the neuronal communication mechanism whereby slower rhythms can reset ongoing neurophysiological processes taking place within a different neuronal population (Buzsáki, 2010; Canolty and Knight, 2010; Buzsáki et al., 2013). Further analyses on PAC data can be used to establish the repertoire of CFC modes that characterize resting activity in a particular group of participants (i.e., the different pairs of frequencies characterizing the modulating and the modulated neural oscillators across sensors). To date there is very little research on the typical patterns of CFC interactions in the normal mature brain, while even less is known regarding the repertoire of prominent CFC modes in typically developing children.

The amount of information being exchanged between two neuronal oscillators may increase as a function of the spectral separation between the two prominent oscillating frequencies (Jensen and Colgin, 2007; Canolty and Knight, 2010; Florin and Baillet, 2015). It is not known, however, whether the typical functional substrate for cortico-cortical signaling features a small number of CFC modes or if a larger repertoire of CFC modes is related to more efficient communication between neuronal assemblies. In this newly emerging field, establishing distinct profiles of CFC modes that vary systematically as a function of reading ability in children would provide important clues regarding the functional significance of PAC metrics.

In the case of neurophysiological activity recorded at rest, CFC patterns should be relatively stable in time if they are to serve as the substrate of functional interactions between specific neuronal assemblies. Although by definition neuronal communication is a dynamic process, only recently appropriate analytical tools were made available to model connectivity patterns as they evolve over time during the relatively restricted time window of a resting activity recording session. Analytic tools developed in the context of the theory of Coordination Dynamics (Tognoli and Kelso, 2014a,b) and used in the present study permit estimation of the probability of change from one prominent CFC mode to another over successive time windows, separately at each sensor and between sensor pairs, through the transition rate index (TR). Next, TR values are aggregated over time windows and sensors to provide global indices of the degree of CFC variability that characterizes resting state recordings of each participant. The derived indices can then be used to classify participants based on reading-group membership and further validated through associations with reading achievement measures.

This study had three aims: First, to identify prominent patterns of CFC interactions in typically developing school-aged children and to determine if similar modes characterize the resting state neuromagnetic recordings of an age-matched group of students who experience severe reading difficulties (RD); Second, to test the prediction that struggling readers would show more variable patterns of CFC over time than typical readers, suggesting a less stable substrate available to support the integration of signaling between neuronal assemblies. Finally, the potential functional significance of putative aberrant dynamic patterns of CFC for reading achievement was assessed through brain-behavior correlational analyses.

Materials and Methods

Participants

The target group included 25 children (12 boys) with reading difficulties (RD group) as indicated by scores below the 16th percentile level (standard score of 85) on the Basic Reading composite index (average of Word Attack and Letter–Word Identification subtest scores of the Woodcock–Johnson Tests of Achievement-III; Woodcock et al., 2001; WJ-III). These scores would be consistent with the presence of dyslexia and is lower than in previous studies (Rezaie et al., 2011; Simos et al., 2011) of this cohort because we wanted to focus on severely impaired children. They were recruited from a larger Grade 6–8 educational study (Vaughn et al., 2010) at-risk for further academic failure because they did not pass the school- administered Texas Assessment of Knowledge and Skills (TAKS). The TAKS is a criterion-referenced reading comprehension assessment that represents the state accountability test. Additional standardized decoding assessments were administered to select this sample of participants because failure on the TAKS is not always because of decoding difficulties. Because entire middle schools were screened, this should ensure better representation of socioeconomic (SES), ethnic backgrounds, achievement, and general cognitive ability scores than recruiting through a clinic. While all of the struggling readers were invited and actively recruited to participate through letters to parents and follow-up phone calls, response rates for imaging studies in school—identified samples are typically very small (e.g., Simos et al., 2007; Meyler et al., 2008).

A second group of 27 children (9 boys) who had never experienced difficulties in reading (NI group) served as comparisons having standard scores >90 on the Reading Composite index (corresponding to the 36th percentile). These students were recruited from the same schools as the RD group in an attempt to control for educational history, ethnicity, and SES factors. All participants had Full Scale IQ scores >80 on the Wechsler Abbreviated Scale of Intelligence (Wechsler, 1999).

Detailed demographic and psychoeducational data are presented in Table 1 and graphically in Section 10 in Supplementary Material. The two groups were comparable on age, ethnicity, handedness, and Performance IQ (p > 0.1 in all cases). As expected the RD group scored significantly lower than the NI group on measures of reading, spelling, and Verbal IQ. Additionally, participants were selected for inclusion in either group only if they had T scores below 65 on the Attention and Hyperactivity CBCL scales (Child Behavior Checklist, Parental Form; Achenbach, 1991) or a mean score lower than 1.67 points on the parental ADHD-C SNAP-IV scale (Swanson et al., 2005), indicating low risk for ADHD (Chen et al., 1994). Informed assent and consent forms were signed by all participating children and their parents or guardians, respectively. The study had been approved by the University of Texas, Institutional Review Board.

TABLE 1

Table 1. Demographic and psychoeducational data for typical (NI) and students with severe reading difficulties (RD).

MEG Recordings

Three minutes of continuous brain activity was acquired with a whole-head neuromagnetometer array (4-D Neuroimaging Magnes WH3600), consisting of 248 first-order axial gradiometer coils, housed in a magnetically shielded chamber. Participants were placed in a supine position and instructed to keep their eyes closed during the recording. Data were originally collected at a sampling rate of 1017.25 Hz and bandpass filtered within 0.1–200 Hz. Axial gradiometer recordings were then transformed to planar gradiometer field approximations using the sincos method available in Fieldtrip (Oostenveld et al., 2011).

Data Preprocessing

Ongoing activity was corrected for artifacts through a two-step procedure implemented in Matlab (The MathWorks, Inc., Natick, MA, USA) and Fieldtrip (Oostenveld et al., 2011). Line noise was first removed using a notch filter at 60 Hz and the single-subject data was whitened and reduced in dimensionality by means of Principal Component Analysis (PCA) with a threshold corresponding to 95% of total variance (Delorme and Makeig, 2004; Escudero et al., 2011; Antonakakis et al., 2013). The resulting signals were submitted to independent components analysis (ICA) using the extended Infomax algorithm as implemented in EEGLAB (Delorme and Makeig, 2004). A given independent component was considered to reflect ocular or cardiac artifacts if more than 30% of its z-score kurtosis or skewness values, respectively, were outside ± 2 of the distribution mean (Escudero et al., 2011; Antonakakis et al., 2013). The remaining ICs were used to reconstruct a relatively artifact-free signal. The average number of artifactual ICs was 5.5 for NI and 5.1 for RD students. To minimize the contribution of volume conduction to the estimates of intrinsic coupling modes (Engel et al., 2013) and to the identification of the dominant cross-frequency pairs, MEG recordings were orthogonalized prior to dynamic cross-frequency coupling analysis (Hipp et al., 2012). Next, the data were bandpass-filtered in the following frequency ranges using a Third-order Butterworth filter (in zero-phase mode): 0.5–4, 4–8, 8–10, 10–13, 13–15, 15–19, 20–29, and 30–45 Hz corresponding to δ, θ, α₁, α_2, β₁, β₂, β₃, and γ bands.

CFC Metric Computation

CFC estimates the strength of pairwise interactions and identifies the prominent interacting pair of frequencies, both between and within sensors (Buzsáki, 2010; Canolty and Knight, 2010; Buzsáki et al., 2013). Among available CFC descriptors, phase-amplitude coupling (PAC), which relies on phase coherence, is the one most commonly encountered in research (Cohen et al., 2009a,b; Voytek et al., 2010). The PAC algorithm, as adapted to continuous MEG multichannel recordings, is described below.

The more general case of CFC (i.e., between-sensor coupling) is described here—within-sensor CFC is derived by collapsing the two sensor indices to a common index. Let x(i_sensor, t) be the MEG activity recorder at the i_sensor-th site, and t = 1, 2,.…T the successive time points. Given two frequency-limited signals x(i_sensor, t) and x(j_sensor, t), cross-frequency coupling is estimated by allowing the phase of the lower frequency (LF) oscillation to modulate the amplitude of the higher frequency (HF) oscillation. The complex analytic representations of each signal z_LF(t) and z_HF(t) are derived via the Hilbert transform (HT[.]).

$\begin{array}{l}{\text{z}}_{\text{LF}}(\text{t})=\text{HT}[{\text{x}}_{\text{LF}}(\text{t})]=\left|{\text{z}}_{\text{LF}}(\text{t})\right|{\text{e}}^{\text{i\hspace{0.17em}}{\varphi }_{\text{LF}}(t)}=A_{\text{LF}}(\text{t}){\text{e}}^{\text{i\hspace{0.17em}}{\varphi }_{\text{LF}}(t)},\\ {\text{z}}_{\text{HF}}(\text{t})=\text{HT}[{\text{x}}_{\text{HF}}(\text{t})]=\left|{\text{z}}_{\text{HF}}(\text{t})\right|{\text{e}}^{\text{i\hspace{0.17em}}{\varphi }_{\text{HF}}(t)}=A_{\text{HF}}(\text{t}){\text{e}}^{\text{i\hspace{0.17em}}{\varphi }_{\text{HF}}(t)}\end{array}$

Next, the envelope of the higher-frequency oscillation A_HF(t) is bandpass-filtered within the range of the LF oscillation and the resulting signal is submitted to an additional Hilbert transform to derive its phase dynamics component ϕ'(t):

$\begin{array}{l}\text{z}\prime (\text{t})=\text{HT}[A_{\text{HF},\text{LF}}(\text{t})]=\left|\text{z}\prime (\text{t})\right|{\text{e}}^{\text{i\hspace{0.17em}}\varphi {\prime }_{\text{HF}}(t)}=\left|\text{z}\prime (\text{t})\right|{\text{e}}^{\text{i\hspace{0.17em}}{\varphi }_{\text{LF}\to \text{HF}}(t)}\end{array}$

which expresses the modulation of the amplitude of the HF by the phase of the LF oscillation. Phase consistency between the two time series is measured by means of both the original definition (Lachaux et al., 1999) and the imaginary portion of PLV (iPLV), as synchronization indexes to quantify the strength of PAC.

The original PLV is defined as follows:

$\begin{array}{ll}PLV=\frac{1}{T}\ast \left|{\displaystyle \sum _{t=1}^T{e}^{i({\varphi }_{LF}(t)-{\varphi }_{HF}(t))}}\right|& (1)\end{array}$

and the imaginary part of PLV as follows:

$\begin{array}{ll}\text{Im}PLV=\frac{1}{T}\ast \left|Im\left({\displaystyle \sum _{t=1}^T{e}^{i({\varphi }_{LF}(t)-{\varphi }_{HF}(t))}}\right)\right|& (2)\end{array}$

The imaginary portion of PLV is considered to be less susceptible to volume conduction effects in assessing CFC interactions and was used in all subsequent analyses. While iPLV is not affected by volume conduction, it could be sensitive to changes in the angle between two signals, which do not necessarily imply a PLV change. In general, iPLV is only sensitive to non-zero-phase lags and is thus resistant to instantaneous self-interactions associated with volume conductance (Nolte et al., 2004). For comparison purposes results based on PLV are shown as Supplementary Material (Section 7).

Figure 1 demonstrates the aforementioned algorithmic steps for a pair of sensors located over left and right frontal areas in a specific 2-s sliding window from the recording in a NI participant. In this example we show PAC interactions between oscillations in the δ band (LF) and oscillations in the β₁ band (HF). The original recordings are shown in Figure 1A. The HF component of signal x1(t) with its envelope is depicted in Figures 1B,D shows the low-pass filtered (within the δ band) version of the previous envelope [signal denoted A_{β1, θ}(t)]. The saw-like trace corresponds to its instantaneous phases ${\varphi }_{\beta }^{\prime }$₁(t). The LF version of the signal x2(t) is depicted in Figure 1C, along with the trace of the corresponding instantaneous phases ϕ_θ(t). The ϕ_θ(t) and ${\varphi }_{\beta }^{\prime }$₁(t) traces are superimposed in Figure 1E, to reveal the instantaneous phase differences which are plotted in Figure 1F. It is this sequence of phase-differences Δϕ(t) that enters in Equation (1) and will be “integrated” across time via averaging the corresponding directional vectors e^iΔϕ(t) in the complex domain. It is important that the length T of this sequence is sufficiently long so that the derived iPLV index provides a reliable estimate of PAC.

FIGURE 1

Figure 1. An example of the algorithmic steps involved in the estimation of iPLV between the phase of δ oscillations recorded from a left frontal sensor (LF signal) and the amplitude of β1 oscillations recorded from a right frontal sensor (HF signal). The two continuous MEG time series are shown in (A). The HF signal filtered in the β₁-band is shown in (B; red trace) with its amplitude modulation envelope superimposed (dark trace). Panel (C) displays the LF signal filtered in the δ-band (blue trace) along with its corresponding instantaneous phases φ'δ(t) (dark trace). The δ oscillations in the HF signal are shown in (D; red trace) with the corresponding instantaneous phases superimposed [φ'δ(t)]. The evolution of instantaneous phases of the LF rhythm overlaid on the corresponding phases of the HF oscillation envelope is shown in (E). The evolution of phase differences between the two signals (Δφ) is used as an index of frequency-specific phase-amplitude coupling (PAC) between the two sensors (F).

For each pair among the 248 MEG sensors, at each sliding window and for each subject, the maximum PAC score was computed for each combination of fφ (frequency band of the phase modulators from δ to γ) and fa (frequency modulated oscillator ranging from 0.5 to 45 Hz in 1 Hz steps. The maximum PAC value for each pair of sensors and at each sliding window in both the amplitude and phase domain is given by the following equation:

$\begin{array}{ll}{(}^{\wedge }{\text{f}}_{\varphi }{,}^{\wedge }{f}_{\alpha })=\underset{{\text{f}}_{\varphi }\text{-}{f}_{\alpha }}{\arg \max }(PAC({\text{f}}_{\varphi },f_{\alpha }))& (3)\end{array}$

In this manner we were able to assess if the prominent CFC interaction (indexed by the maximum PAC value), and specifically the frequency of the low-frequency phase, was associated with the frequency of maximum power. If the two frequencies were identical this would imply that the observed prominent CFC interaction was driven by the power of the dominant frequency and not its phase. In order to ensure that we did not include in further analyses CFC interactions of this type, we required a minimum frequency difference of 1 Hz distance between the two frequencies (frequency of the low-frequency phase, frequency of the higher power), the one identified by Equation 3 (maximum PAC value) and the one associated with the maximum spectral power, within the range of high amplitude. Section 1 in Supplementary Material describes the method employed to assess whether PAC estimates were not biased by the maximum power in more detail.

Dynamic PAC Estimates: The Time-Varying CFC Graph (^TVPAC-Graph)

The goal of the analytic procedures described in this section was to capture the repertoire of CFC interactions and their temporal evolution, while taking into account the quasi-instantaneous spatio-temporal distribution of PAC estimates. This was achieved by computing a set of PAC estimates within each series of sliding, partially overlapping windows spanning the entire 3-min continuous MEG recording. The width of the temporal window was set equal to the duration of two cycles of δ activity (i.e., 2 s) ensuring that modulations of activity by the phase of the lowest frequency band (δ) would be preserved. The center of the sliding window moved forward at 0.5-s steps and the CFC interactions between every possible pair of frequencies were reestimated leading to a temporally integrated PAC graph. In this manner a series of 308 sets of PAC estimates were computed per sensor pair for each participant.

This procedure, the implementation details of which can be found elsewhere (Dimitriadis et al., 2010a), resulted in 217 time-varying PAC graphs per participant (^TVPAC), each serving as an instantaneous snapshot of the surface network. There were 41, 37, 35, 32, 30, 26, and 16 (= 217) possible pairs in the δ, θ, α1, α2, β1, β2, and β3 bands. For instance, using 1 Hz frequency bins there were 45–4 = 41 interacting pairs for δ phase. ^TVPAC tabulates PAC estimates between and within sensors.

Surrogate Data Analysis of PAC Estimates

To identify significant PAC-interactions which were estimated for every pair of frequencies, within and between all 248 sensors, and at each successive sliding window, we employed surrogate data (Theiler et al., 1992). Surrogate data analyses determined: (a) if a given PAC value differed from what would be expected by chance alone, and (b) if a given non-zero PAC indicated coupling that was, at least statistically, non-spurious.

For every time window, sensor-pair, and pair of frequencies, we tested the null hypothesis H₀ that the observed PAC value came from the same distribution as the distribution of surrogate PAC-values. One thousand surrogate time-series ${\varphi }_{LF}^s(t)$ were generated by cutting at single point at a random location and exchanging the two resulting time courses (Canolty et al., 2006; Aru et al., 2015). Repeating this procedure produced a set of surrogates with minimal distortion of the original phase dynamics and impact on the non-stationarity of brain activity as compared to either merely shuffling the time series or cutting and rebuilding the time series in more than one time points. This procedure ensures that the observed and surrogate indices shared the same statistical properties. For each data set the surrogate PAC (^sPAC) was computed. We then determined a one-sided p-value expressing the likelihood that the observed PAC value could belong to the surrogate distribution, and corresponded to the proportion of “surrogate”' PAC^s which was higher than the observed PAC value (Theiler et al., 1992). PAC values associated with statistically significant p-values were considered unlikely to reflect signals not entailing PAC coupling.

The FDR method (Benjamini and Hochberg, 1995) was employed to control for multiple comparisons (across all possible pairs of frequencies) with the expected proportion of false positives set to q ≤ 0.01. Finally, the PAC mode that characterized a specific pair of frequencies was determined based on the highest, statistically significant PAC value from surrogates. To be included in further analyses a given PAC value should meet an additional criterion, namely that the frequency of the modulated HF should be distinct from the frequency associated with maximum power (see Section 2.4. and Section 1 in Supplementary Material). Here, we adopted a criterion of at least 1 Hz difference between the two frequencies.

Finally, the significant, dominant PAC values for each pair of sensors and across sliding windows were integrated within each frequency band yielding 28 possible pairwise PAC estimates among the eight frequency bands. In the example of Figure 2A only one frequency pair met this criterion. In the case of two frequency pairs both exceeding the statistical threshold (Figure 2B), the one with the highest PAC value was identified as the characteristic PAC mode for this pair of sensors at that particular time window. If none of the cross-frequency pairs exceed the statistical threshold (Figure 2C), a value of zero was assigned to this pair of sensors with not identified characteristic coupling mode. For each participant the resulting ^TVPAC profiles constituted a 4D array of size [28 (pairs of frequencies) × 308 (time windows) × 248 (sensors) × 248 (sensors)]. The identity of prominent frequency pairs for every pair of sensors (including within-sensor interactions) at each time window was finally stored in a second 4D array of size [28 × 308 × 248 × 248]. In the latter array significant PAC interactions were indicated by a value of 1, with zeros indicating non-significant PAC interactions.

FIGURE 2

Figure 2. Three different scenarios in the process of identifying prominent PAC modes associated with a particular pair of sensors at three different time windows (A–C). The modulating phase (LF) is plotted on the horizontal axis. iPLV values for each frequency pair are shown on the diagonal and are marked by (*) if they exceed the statistical threshold. A significant θ–α₁ PAC value is shown in (A), two significant frequency pairs (θ–α₁ and α₁–α₂) were found in (B) with the former identified as the prominent pair, whereas PAC values failed to reach significance during the third window (C). The sensor-frame location where data originated from is marked in (D). The modulation of iPLV values across the entire MEG time series (successive time windows t1 …tN) is also shown in (D). At t1 the phase of the LF (δ) of the signal recorded at a left frontal sensor appears to modulate θ amplitude at a right frontal sensor. At t2, there is modulation of right frontal γ1 amplitude by the phase of left frontal θ oscillations. At t3 this PAC mode persists, although its strength is higher. Subsequently (t4) the prominent PAC mode involves modulation of γ2 amplitude by β1 phase. For the two sensors shown here, the last time window (tN) is characterized by modulation of γ2 amplitude by the phase of α1 oscillations.

Group Comparisons on CFC Metrics

Dynamic Strength of ^TVPAC

Participants in the RD and NI groups were initially compared on the time-series of grand total PAC strength aggregated over all pairs of sensors. The dynamic strength of ^TVPAC was estimated by aggregating the dominant cross-frequency pair and its strength estimated via iPLV over all pairs of sensors for each sliding window. We adopted two complementary indices, namely Dynamic Information Exchange Rate (dIER) and Dynamic Weighted Information Exchange Rate (wdIER), to aggregate both the dominant CFC mode and the strength over all pairs of sensors. Each approach resulted in a single time series per participant describing the temporal evolution of information exchange rate.

Repertoire of Characteristic Cross-Frequency Interactions

The two groups were further compared on the repertoire of prominent PAC modes which were derived with the procedure described in section 2.4.2. In this step, in order to reduce the number of comparisons, a single 7 × 7 matrix was obtained per participant by aggregating over all successive time windows and sensor pairs. The resulting empirical probability distributions (PDs) of all interacting pairs of frequencies were then averaged across participants to identify the prominent cross-frequency phase-amplitude coupling pairs for each group. Additionally, we constructed separate comodulograms (7 × 7 matrices) aggregated over time or sensor space (see Figure S2 in Supplementary Materials).

Transition Dynamics Based on Dynamic Prominent CFC Profiles

Next, ^TVPAC measurements were used in order to assess group differences in the stability of PAC modes. An index of the dynamic characteristics of CFC was computed in the form of the transition rate (TR) defined as the number of times that a given PAC profile at latency t (e.g., θ → β1) changed during the next latency into a different characteristic interaction (e.g., θ → γ), divided by the total number of possible transitions. TR can take values between 0 and 1 and is defined as follows:

$\begin{array}{rc}TR=\frac{numberoftransitionsbetweenCFCpairs}{numberoftemporalsegments-1}& (4)\end{array}$

Using this formula a single TR value was computed for each participant, sensor (within-sensor CFC) and sensor pairs (between-sensor CFC). An example of cross-frequency pairs fluctuating over time is presented in Figure 2D, where the TR for the first four time windows would be TR = 2/(4–1) = 2/3 = 0.66. Transitions were counted only between two latencies where a significant CFC was detected.

Dynamic Information Exchange Rate (dIER)

Under the notion that PAC modes reflect processes through which “packets of neural information” are exchanged between populations of neurons operating at different characteristic frequencies, we developed two novel measures to summarize the rate of information transfer between neural assemblies throughout the brain. These information theory estimates were based on the dominant PAC mode for each pair of sensors over time. The “instantaneous” (dynamic) Information Exchange Rate (dIER) is the sum of the number of cycles of the higher frequency that can be nested within the phase of the slower frequency across prominent frequency pairs and sensor pairs at each time window. The dIER index is defined as follows:

$\begin{array}{rc}dIER(t)={\displaystyle \sum _{i_{sensor}=1}^{Noofsensors}}{\displaystyle \sum _{j_{sensor}=1}^{Noofsensors}}\frac{NoofcyclesofHF}{NoofcyclesofLF}& (5)\end{array}$

A weighted version of dIER (wdIER) was also computed that takes into account the strength of dominant PAC modes as indexed by the corresponding iPLV, with higher values implying stronger modulation of higher frequency rhythms by lower frequency oscillations throughout the recording epoch. Thus, in addition to providing information regarding how many cycles of the higher frequency can be “travelled” within the cycle of the lower frequency, wdIER also indexes how active is a “channel” between or within sensors for information exchange at quasi-stable time instances, and is defined as follows:

$\begin{array}{rc}wdIER(t)={\displaystyle \sum _{i_{sensor}=1}^{Noofsensors}}{\displaystyle \sum _{j_{sensor}=1}^{Noofsensors}}\frac{NoofcyclesofHF}{NoofcyclesofLF}*iPLV& (6)\end{array}$

In the example of Figure 1, there were 2*14 cycles of β₁ oscillations (HF) and 2 cycles for δ (LF) within the window of 2 s with a corresponding IER = 28/2 = 14. With a iPLV value of 0.21, wdIER = IER*iPLV = 14*0.21 = 2.94. The frequency bin (1 Hz resolution) of the faster oscillation where the maximum PAC was detected was used to compute this non-integer ratio of the two frequencies. A cycle of recorded activity in a particular frequency range corresponds to the period of the corresponding sinusoidal oscillation (see section 11 in supp. material).

Classification of RD and NI Readers Based on dIER/wdIER, TR Profiles and PD^PAC of Prominent CFC

Four complementary features were validated for their capacity to discriminate between NI and RD groups using a k-Nearest Neighbor (k-NN) classifier: the dynamic information exchange rate (dIER and wdIER indices) for each participant, the TR profiles, and the probability distribution (PD^PAC) of prominent CFC (see Figure 4). The corresponding feature-vectors (FVs) input to the k-NN classifier consisted of: the dIER and wdIER time series (each consisting of 308 samples corresponding to the number of time windows), of a 248 × 248 matrix of TR values, and of a vector consisting of 28 PD^PAC values.

Classification performance was assessed in four adjunct validation schemes. The first scheme entailed leave-one-out cross-validation where each participant was considered of unknown group membership. On each run, the k-NN classifier was trained on data from a subset of 51 participants (27 NI + 25 RD –1) and tested on its capacity to correctly classify the remaining participant based on each FV separately (dIER, wdIER, TR, or PD^PAC).

The second scheme was a 2-fold cross-validation, where the k-NN classifier was first trained on data from a randomly created training data set from 26 participants with known group membership and tested for classification accuracy on the data from the remaining 26 participants. Again classification accuracy was tested separately for each of the four FVs.

The third and fourth validation schemes were conducted in order to evaluate the within-subject reliability of the four sets of features (dIER or wdIERs, TR and PD^PAC) for discriminating RD from NI students. Initially, the dIER and wdIER time series were split into two segments of equal length and TR and PD^PAC values were recalculated for each segment. FVs input to the k-NN classifier consisted of either the 154D split-segment dIER and wdIER time series, the 248D TR value matrices, or the 28 PD^PAC values from each time segment and participant. The third validation scheme employed a leave-one-out procedure where indices derived from 51 × 2 segments served as the training set and the two segments from the remaining participant as the testing set. In the 2-fold cross-validation, training and testing sets included data from 52 randomly mixed A and B segments from each participant with the remaining 52 segments serving as the testing set.

To properly account for the dynamics of information exchange rate represented by the dIER and wdIER time series, each 308D set (or 154D set in the split-temporal domain schemes) was first transformed into a dynamic trajectory using a time-delay embedding procedure and used to compute the Wald-Wolfowitz (WW) Dissimilarity Index (w; described in more detail in Section 8 of the Supplementary Material).

Results for each classification scheme were pooled over 200 repetitions employing different training and testing sets each time. The Euclidean distance served as the optimization criterion in classification runs based on TR and PD^PAC profiles, whereas the Dissimilarity Index was employed for classification runs conducted on the dIER and wdIER data sets, as a more appropriate distance estimator between time series (Friedman and Rafsky, 1981; Rigas et al., 2015; see section 8 in Supplementary Material). Classification results were robust for k parameter values ranging between 7 and 11.

Predictability of Dynamic Information Flow

To access the predictability (or randomness) of dynamic information flow in both groups, we also estimated sample entropy (Richman and Moorman, 2000). Sample entropy was estimated over both dIER and wdIER indices in both groups while the pair of embedding dimensions and time delays were optimized for each time series accordingly (see Sections 8 and 9 in Supplementary Material).

Results

Dynamic Strength of ^TVPAC

To ensure that PAC values were not driven by the power of the low-frequency signal, we estimated the difference between: (a) the low frequency involved in computing each PAC value and (b) the frequency associated with the highest power in the signal, separately for each time window, sensor, and sensor pair (see Section 1, Supplementary Material). In more than 90% of all instances across subjects, the two frequencies differed by at least 1 Hz and in more than 70% of instances by at least 2 Hz. This implies that the frequency of the slower oscillator rarely corresponded to the dominant frequency in the power spectrum, supporting the claim that PAC reflects a phase-modulated phenomenon and is not driven by the amplitude of the dominant frequency.

Figure 3 presents the group-averaged traces of dynamic ^TVPAC strength showing higher PAC levels for RD children. It is important to mention here that, since we are dealing with resting state data, the underlying brain activity is not aligned across subjects. Although group averaging can only be used for demonstration purposes, the distributions of individual ^TVPAC values were clearly different between groups (Wilxocon test, p < 10⁻¹²).

FIGURE 3

Figure 3. Group-averaged strength of dynamic CFC graphs showing significantly greater overall strength in the RD compared to the NI group (Wilcoxon score, p < 10⁻¹³). The sums of PAC values over the entire network (between every possible pair of sensors) are plotted as solid lines surrounded by the corresponding ± 1 SD ranges.

The Repertoire of Characteristic CFC Interactions

Prominent frequency pairs (i.e., PAC modes) across all sensors, sensor pairs, and time windows were summarized for each subject in the form of an 8 x 8 matrix. Group-averaged for demonstration purposes, these empirical probability distributions presented in Figure 4 reveal a richer repertoire of CFC-interactions among RD children. In contrast, CFC interactions typical of NI children occur between brain rhythms in spectrally adjacent frequency bands (e.g., δ → θ, θ → α1 etc).

FIGURE 4

Figure 4. Group-averaged empirical Probability Distribution values of cross-frequency phase-amplitude coupling pairs, revealing a broader repertoire of CFC pairs for RD (B) compared to NI (A) participants who showed a more restricted range of prominent CFC pairs. Lower PD values indicate more significant iPLV values aggregated over time windows and sensor pairs.

Transition Dynamics of Prominent CFC Interactions

The two reading groups differed not only on the cross-frequency distributions of prominent PAC modes but also on their stability over time. TR values (between and within sensors) were near zero for the majority of NI children exceeding the threshold of ~0.003 (corresponding to 2 transitions across steps of the sliding windows) in only three participants. In contrast, TRs frequently exceeded the empirically-derived threshold in RD children. The topographic distribution of between-sensor TR values was visualized by applying a threshold of (mean + 2 SDs) to the set of 248 × 248 TRs in the group-averaged TR data. Figure 5 shows sensor pairs associated with “significant” TR values topographically organized by “lobe” and hemisphere. The figure indicates significant temporal variability in the CFC interactions driving right hemisphere brain oscillations by (lower frequency) left hemisphere brain rhythms. Temporally variable CFC interactions driven by right hemisphere oscillations were restricted within the same hemisphere. Significant temporal variability of within-sensor CFC interactions was also more prominent at right hemisphere sensors as shown in Figure 6.

FIGURE 5

Figure 5. Group-averaged transition rate (TR) for pairs of sensors located roughly over the four lobes (A–D and E–H) in each hemisphere for RD children. For better visualization results were thresholded based on mean +2SD from the set of 248 × 248 TR values in order to reveal the most significant TRs over the entire sensor space. Corresponding TR values failed to reach significance in the NI group (not shown).

FIGURE 6

Figure 6. Group-averaged distribution of transition rates (TR) for within-sensor PAC-interactions in the group of RD participants. Thresholding (based on mean +2 SDs of the 248 TR values) has been applied revealing sensors located over right hemisphere temporo-parietal and occipito-temporal regions and the midline, that show maximum instability of local PAC interactions.

Dynamic Information Exchange Rate

As expected on the basis of the aforementioned group differences on TR values, the RD group showed significantly higher values on both global measures of dynamic information exchange rate (dIER and wdIER; Wilcoxon test, p < 10⁻¹², Figure 7).

FIGURE 7

Figure 7. Group-averaged dynamic (dIER; A) and weighted dynamic (wdIER; B) Information Exchange Rate showing significant differences between the RD and NI groups (Wilcoxon score, p < 10⁻¹¹). Value sums over the entire network (between every possible pair of sensors) are shown on the vertical axis. Thinner lines indicate ± 1 SD intervals.

Classification of RD and NI Readers Based on Dynamic Information Exchange Rates

Results of each of the adopted cross-validation schemes indicated perfect discrimination between the two groups based on both dIERs and wdIERs profiles (see Section 8 in Supplementary Material). Similar results (i.e., 100% classification accuracy) were obtained when using the 28D vector from the individual probability distribution (PD) of prominent PAC modes (PD^PAC) and also the TR profiles (see Figure 4). For comparison, classification accuracy based on relative power in the same frequency bands used to estimate CFC did not exceed 70% (see Section 2 in Supplementary Material).

Associations between TR Values and Cognitive Measures

To assess the functional significance of indices of dynamic information exchange developed in the present work, we computed partial correlations between reading achievement scores and TRs (adjusting for age) computed for pairs of sensors where significant group differences were detected (sensors located over right temporo-parietal areas). Significant correlations, restricted to the RD group for WJ3-Reading Composite scores, were found at five sensor pairs, the topography of which is shown in Figure 8. Significant partial correlations between reading achievement/IQ measures and TR values from the five sensors located over right temporo- parietal sensors are tabulated in Table 2.

FIGURE 8

Figure 8. Topography layout of the five pairs of sensors, located over right temporo-parietal regions, where significant positive correlations were found between TR values and reading measures among RD participants.

TABLE 2

Table 2. Significant partial correlations between reading achievement/IQ measures and TR or PD values at right temporo- parietal sensors (RD group).

Associations between PD Values and Cognitive Measures

Partial correlation coefficients (adjusting for age) were computed between PDs integrated over time and sensors located above right temporo-parietal areas and reading achievement/IQ scores. Significant positive associations (p < 0.01) were found between WJ3-Reading Composite scores and PD values for three frequency pairs: δ-θ, θ-γ, and α2-γ which were restricted to the RD group (Figures 9, 10). Among NI participants there was a moderate-size correlation (r = 0.61) which, however, did not reach significance (Figure 9A). Significant partial correlations between reading achievement/IQ measures and PD values for the three frequency pairs are tabulated in Table 2.

FIGURE 9

Figure 9. Regression scatter plots of Woodcock–Johnson III (WJ3) Reading Composite scores over Probability Distribution values for three frequency pairs [(A) δ-θ, (B) θ-γ, (C) α2-γ] aggregated over right temporo-parietal sensors among RD participants (blue circles; B,C). A significant association noted for NI participants (green circles) is included in (A).

FIGURE 10

Figure 10. Regression scatter plots of Woodcock–Johnson III (WJ3) Reading Composite scores over Probability Distribution values aggregated over right temporo-parietal sensors among RD participants.

Predictability of Dynamic Information Flow Based on Dominant Cross-Frequency Pairs

Sample entropy values estimated over dIER and wdIER were significantly higher for NI as compared to the RD group implying lower overall predictability of CFC interactions in the latter group (p < 0.001, Wilcoxon rank-sum test; see Figure 11). Interestingly, there was a significant increase in entropy values in this group, when estimated for sensors located over the right temporal brain region (where the highest transition rates were originally observed).

FIGURE 11

Figure 11. Average entropy estimates derived from dIER and wdIER values aggregated across all sensors for the RD (blue bars) and NI groups (red bars). Notice the significant increase in entropy in the RD group after excluding sensors located over right temporo-parietal regions (RD-RT; green bars). *p < 0.001, Wilcoxon rank-sum test. Standard error in bars.

Discussion

Here we report results from a novel analytic approach to study resting-state functional connectivity derived from neuromagnetic data in sensor space. The first key finding is that the prominent mode of phase-amplitude coupling between neuronal assemblies in RD children at rest involves rhythmic oscillations that vary widely in frequency. In contrast, the typical mode of CFC interactions among NI readers involved driving the amplitude of a particular brain oscillation by the phase of oscillatory activity that differed very little in frequency from the former. The second notable finding is that the spectral mode of these interactions was highly stable among typically-achieving readers during the 3-min recording session. In contrast, cross-frequency interactions in RD students showed high levels of temporal instability and lower levels of entropy (predictability). Correlational analyses revealed positive associations between indices of dynamic CFC interactions in recordings over right hemisphere temporo-parietal brain areas and reading achievement scores among struggling readers. The proposed dynamic connectomic features (dIER, wdIER, TR, and PD^PAC) appeared sufficient to discriminate NI from RD students with a high degree of accuracy.

Repertoire of CFC Modes

Phase-amplitude phenomena, also referred to as “nested oscillations,” occur when the amplitude of an oscillation at a particular frequency is modulated by the phase of a lower-frequency oscillation. This phenomenon has been demonstrated in various species and brain areas and is postulated to be an important attribute of the coordinated activity taking place in distinct neuronal populations, ensuring efficient transfer of information (i.e., signaling that is accurate in timing and high in gain; Varela et al., 2001; Buzsáki, 2006; Buzsáki and Watson, 2012; Buzsáki et al., 2013). Such functional connectivity patterns that characterize resting activity may serve as the dynamic substrate of complex cognitive functions, such as reading, which involving multisensory integration and coordination of several distinct component operations (visual/graphemic processing and recognition, phonological processing, and articulation).

Different oscillation frequencies are characterized by different spatial scales for the interacting neuronal populations, with lower frequencies having a higher capacity to modulate activity in remote brain regions than higher frequencies (von Stein and Sarnthein, 2000). The spectral separation between the two oscillating frequencies may also be important for efficient information exchange (Buzsáki and Watson, 2012) given that, at last in principle, the number of cycles of the higher frequency encapsulated within the phase of the slower frequency is related to the amount of information being exchanged between the two oscillators (Jensen and Colgin, 2007; Canolty and Knight, 2010; Florin and Baillet, 2015). Cross-frequency phase-to-amplitude coupling may have evolved as a compensatory mechanism that enables efficient information transmission without the need to slow-down the processing rate of the receiving neuronal population to match that of the signal source (Buzsáki and Watson, 2012; Buzsáki et al., 2013). Presently we can only surmise that the substrate of relatively stable CFC-related dynamic information exchange observed among children with normal developmental histories reflects a typical and, presumably, optimal working level ensuring efficient neuronal transmission (Deco and Corbetta, 2011; Deco et al., 2013).

Temporal Variability of CFC Profiles

An example of the increased temporal variability of CFC profiles characterizing RD participants, is that for a given pair of sensors A, B the amplitude of the neuromagnetic signal in the β1 frequency band in sensor B may be initially driven by the phase of δ oscillations in sensor A. During subsequent time intervals, activity in sensor B may be characterized by the driving of γ amplitude by α₁ phase in sensor A whereas, only seconds later, α₂ amplitude in sensor B may be driven by the phase of θ oscillations in sensor A. The same was true for cross-frequency interactions at each sensor separately. This type of temporal variability was summarized by global indices of temporal variability in CFC interactions (dIER and dwIER profiles) which clearly differentiated NI from RD students with very high levels of accuracy in the context of different, complementary cross-validated classification schemes. A complementary, global index of predictability in CFC interactions, namely entropy, corroborated these findings by demonstrating higher randomness in the recordings obtained from RD participants.

Importantly, the patterns of temporal variability in cross-frequency interactions among struggling readers showed distinct topographies. PAC modes derived from activity recorded by each sensor separately were indicative of highly unstable cross-frequency interactions between neuronal assemblies located in right hemisphere temporo-parietal and occipito-temporal regions. Between-sensor analyses complemented this finding, suggesting the presence of highly variable CFC interactions of neuronal assemblies distributed within the right hemisphere (e.g., between frontal-temporal, and between temporal-occipital regions). Significant temporal variability was also found in the CFC interactions between left hemisphere occipital, frontal, and temporal sites “driving” higher frequency oscillations in corresponding right hemisphere sensors. Entropy findings further point to right temporo-parietal sensors in RD students as showing the highest levels of randomness in CFC interactions in the entire sensor space.

Although firm conclusions regarding the functional significance of these findings cannot be drawn at present, results raise a number of intriguing possibilities regarding the relevance of ongoing patterns of neuronal coupling at rest for task-related engagement of cortical areas during reading in RD. Overexitability of right hemisphere areas at rest may account for previous reports concerning task-related overactivation of right hemisphere “reading” areas in this population (Pugh et al., 2001; Milne et al., 2002; Hoeft et al., 2007, 2011; Simos et al., 2011). Both the extended repertoire of PAC interactions and their pronounced temporal variability may suggest an inefficient substrate for neuronal communication within and between right hemisphere regions. Correlational analyses support this notion by revealing positive associations between both the degree of temporal variability and the strength of temporally aggregated CFC interactions over right temporo-parietal sensors and reading achievement scores in the group of struggling readers. The fact that positive associations were found implies a compensatory mechanism, which is nevertheless not sufficient to support efficient reading ability. Moreover, the relative specificity to reading measures suggests a more prominent role of this mechanism for reading (as opposed to more general language or cognitive ability). The frequency range of CFC interactions associated with these significant correlations (δ-θ, θ-γ, and α2-γ) suggests a pattern of spectrally widespread alterations in the resting-state substrate of functional interactions in the RD group. These results complement our previous finding suggesting increased local efficiency of cortico-cortical interactions over left temporo-parietal regions in typical as compared to struggling readers revealed by resting-state synchronization measures (Dimitriadis et al., 2013b). Interestingly, the positive correlations between left hemisphere synchronization efficiency indices and reading achievement were restricted to the group of typical readers.

In more general terms, these results are consistent with an emerging view in cognitive neuroscience that the brain is a self-organizing complex system which, under normal circumstances, displays a limited range of variability of spatiotemporal configurations (Sporns et al., 2004; Shanahan, 2010; Deco and Corbetta, 2011; Tagliazucchi et al., 2012; Carhart-Harris et al., 2014; Hansen et al., 2015; Dimitriadis et al., 2015a). Higher levels of variability and temporal instability of signaling modes, such as those characterizing the neuromagnetic activity of RD students at rest, may signify a critical level of function. Under certain conditions, for instance when struggling readers are required to perform more demanding, yet age-appropriate reading tasks, the extensive brain network that is normally responsible for decoding and word recognition is more likely to operate in a transition zone between order and chaos, resulting in suboptimal performance (Tononi et al., 1994; Shanahan, 2010; Deco and Jirsa, 2012; Tagliazucchi et al., 2012).

Limitations of the Study

The present results should be considered in light of a number of methodological limitations. First and foremost functional connectivity estimates were obtained at the sensor level therefore limiting the type of conclusions that can be drawn regarding the efficiency of interactions between specific, underlying neuronal populations. Moreover, volume conduction effects can produce artificial results even in the case of CFC, which is considered to be less susceptible to volume conduction effects (Engel et al., 2013). However, several steps were taken in order to minimize the effect of volume conduction in the calculation of CFC estimates, namely transformation of axial gradiometer data to planar, orthogonalization of the signals and use of the imaginary portion of PLV. The success of these measures is indicated by negligible correlations between PAC time-series computed between adjacent sensors (see Figure S3 in Supplementary Material). Failure to find substantial correlations between relative power and PAC strength in the temporal domain (see Figure S4 in Supplementary Material) lends further support to the claim that our results were not contaminated by volume conduction.

A second issue to be considered when interpreting the current results concerns the sensitivity to low-frequency phase components with asymmetric cycles, a phenomenon noted in intracortical recordings (Steriade et al., 1996; Buzsáki and Draguhn, 2004; Haider et al., 2006). Moreover, edges in the data may cause spurious PAC phenomena (Kramer et al., 2008). It should be noted however that PAC has been found in numerous studies using ECoG, LFP, and MEG data (Canolty et al., 2006; Tort et al., 2010; Özkurt and Schnitzler, 2011) making it unlikely that this phenomenon is artifactual. The application of very steep bandpass filters (Third order forward and backward Infinite Impulse Response filter) served as an additional step in the present study against the likelihood of spectral cross-contamination between adjacent frequency bands. Moreover, the spectral separation between the low (modulating) frequency and the frequency associated with the highest power in the signal was probably sufficient to rule out the possibility that PAC phenomena where induced by high-power low frequency oscillations dominating the resting state recordings.

Strengths of the Study and Future Directions

Despite these limitations, the present study has a number of novel features. The main strength of the study lies in the implementation of a dynamic functional connectivity methodology (Dimitriadis et al., 2009, 2010a,b, 2012a,b, 2013a,b, 2015a,b,c; Deco et al., 2011, 2013; Damaraju et al., 2014; Kopell et al., 2014) to explore how resting-state is shaped by predominant cross-frequency coupling (Jirsa and Müller, 2013; Florin and Baillet, 2015). This methodology was adapted to MEG data that are uniquely suited to provide estimates of instantaneous information exchange between neuronal oscillators affording the requisite subsecond temporal resolution.

Although analyses were conducted in sensor space, MEG recordings were first converted to planar gradiometer field approximations, ensuring that the maximum signal lies directly above active brain areas. Moreover, our method represents a step forward from traditional investigations of functional connectivity, by taking into account temporal variability in cross-frequency coupling (Dimitriadis et al., 2009, 2010a,b, 2012a,b,c, 2013a,b, 2015b, 2016; Deco et al., 2011, 2013; Damaraju et al., 2014). Importantly for future applications, PAC indices appeared to be both sensitive and selective to the presence of developmental reading difficulties, and were associated with very effective discrimination of RD from age-matched typical readers. In contrast, relative power measures were substantially less effective in discriminating participants from the two groups.

It would be interesting in the future to study both NI and RD participants with the same methodology during experimental paradigms tailored to dyslexia and to compare the topologies, the extracted PAC-oriented markers and the information transmission between rest and cognitive tasks. Moreover, a source localization strategy would improve interpretability of the findings (Ioannides et al., 2012; Florin and Baillet, 2015) in relation to indices of anatomical connectivity DTI (e.g., Richards et al., 2015). At a more basic level, it is important to establish developmental trends in the CFC (and especially PAC) resting-state patterns and examine their functional significance through brain-behavior associations.

To conclude, resting-state phase-amplitude coupling metrics developed in the present study are some of several potentially useful indices of atypical functional brain organization in developmental reading disability. Their functional significance, specificity to component reading processes, which are known to be underdeveloped in struggling readers, and their sensitivity to educational interventions requires further systematic investigation, before their value as diagnostic or prognostic tools can be fully evaluated.

Author Contributions

Conception of the research: SD, PS, JF, AP; Methods design and data analysis: SD; Drafting the manuscript: SD, NL; Critical revision of the manuscript: PS, JF, AP; Every author read and approved the final version of the manuscript.

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.

Acknowledgments

This research was supported in part by grant P50 HD052117 from the Eunice Kennedy Shriver National Institute of Child Health and Human Development (NICHD). The content is solely the responsibility of the authors and does not necessarily represent the official views of the NICHD or the National Institutes of Health. This research was also supported in part by Greek State Scholarships Foundation (IKY) as part of post-doc research fellowship of SD.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fnhum.2016.00163

References

Achenbach, T. M. (1991). Manual for the Child Behavior Checklist/4–18 and 1991 Profile. Burlington, VT: University of Vermont Department of Psychiatry.

Google Scholar

Antonakakis, M., Dimitriadis, S. I., Zervakis, M., Rezaie, R., Babajani-Feremi, A., Micheloyannis, S., et al. (2016). Detecting mild traumatic brain injury from resting-state MEG recordings based on Cross-frequency Interactions. Int. J. Psychophysiol. 102, 1–11. doi: 10.1016/j.ijpsycho.2016.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Antonakakis, M., Giannakakis, G., Tsiknakis, M., Micheloyannis, S., and Zervakis, M. (2013). “Synchronization coupling investigation using ICA cluster analysis in resting MEG signals in reading difficulties,” in Bioinfo and Bioeng (BIBE), IEEE 13th International Conference (Chania), 1–5.

Google Scholar

Aru, J., Aru, J., Priesemann, V., Wibral, M., Lana, L., Pipa, G., et al. (2015). Untangling cross-frequency coupling in neuroscience. Curr. Opin. Neurobiol. 31, 51–61. doi: 10.1016/j.conb.2014.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Axmacher, N., Cohen, M. X., Fell, J., Haupt, S., Dümpelmann, M., and Elger, C. E. (2010a). Intracranial EEG correlates of expectancy and memory formation in the human hippocampus and nucleus accumbens. Neuron 65, 541–549. doi: 10.1016/j.neuron.2010.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Axmacher, N., Henseler, M. M., Jensen, O., Weinreich, I., Elger, C. E., and Fell, J. (2010b). Cross-frequency coupling supports multi item working memory in the human hippocampus. Proc. Natl. Acad. Sci. U.S.A. 107, 3228–3233. doi: 10.1073/pnas.0911531107

PubMed Abstract | CrossRef Full Text | Google Scholar

Bassett, D. S., and Bullmore, E. T. (2009). Human brain networks in health and disease. Curr. Opin. Neurol. 22, 340–347. doi: 10.1097/WCO.0b013e32832d93dd

PubMed Abstract | CrossRef Full Text | Google Scholar

Bassett, D. S., Nelson, B. G., Mueller, B. A., Camchong, J., and Lim, K. O. (2012). Altered resting state complexity in schizophrenia. Neuroimage 59, 2196–2207. doi: 10.1016/j.neuroimage.2011.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery ratea practical and powerful approach to multiple testing. J. R. Stat. Soc. B Stat. Methodol. 57, 289–300. doi: 10.1093/acprof:oso/9780195301069.001.0001

CrossRef Full Text | Google Scholar

Braun, U., Muldoon, S. F., and Bassett, D. S. (2015). On human brain networks in health and disease. Encyclopedia Life Sci. (eLS-Wiley). 1–9. doi: 10.1002/9780470015902.a0025783

PubMed Abstract | CrossRef Full Text | Google Scholar

Brookes, M. J., Woolrich, M., Luckhoo, H., Price, D., Hale, J. R., Stephenson, M. C., et al. (2011). Investigating the electrophysiological basis of resting state networks using magnetoencephalography. Proc. Natl. Acad. Sci. U.S.A. 108, 16783–16788. doi: 10.1073/pnas.1112685108

PubMed Abstract | CrossRef Full Text | Google Scholar

Buckner, R. L., Andrews-Hanna, J. R., and Schacter, D. L. (2008). Anatomy, function, and relevance to disease. Ann. N. Y. Acad. Sci. 1124, 1–38. doi: 10.1196/annals.1440.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G. (2006). Rhythms of the Brain. New York, NY: Oxford University Press.

Google Scholar

Buzsáki, G. (2010). Neural syntax: cell assemblies, synapsembles, and readers. Neuron 68, 362–385. doi: 10.1016/j.neuron.2010.09.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G., and Draguhn, A. (2004). Neuronal oscillations in cortical networks. Science 304, 1926–1929. doi: 10.1126/science.1099745

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G., Logothetis, N., and Singer, W. (2013). Scaling brain size, keeping timing: evolutionary preservation of brain rhythms. Neuron 80, 751–764. doi: 10.1016/j.neuron.2013.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G., and Watson, B. O. (2012). Brain rhythms and neural syntax: implications for efficient coding of cognitive content and neuropsychiatric disease. Dialogues Clin. Neurosci. 14, 345–367.

Google Scholar

Canolty, R. T., Edwards, E., Dalal, S. S., Soltani, M., Nagarajan, S. S., Kirsch, H. E., et al. (2006). High gamma power is phase-locked to theta oscillations in human neocortex. Science 313, 1626–1628. doi: 10.1126/science.1128115

PubMed Abstract | CrossRef Full Text | Google Scholar

Canolty, R. T., and Knight, R. T. (2010). The functional role of cross-frequency coupling. Trends Cogn. Sci. 14, 506–515. doi: 10.1016/j.tics.2010.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Carhart-Harris, R. L., Leech, R., Hellyer, P. J., Shanahan, M., Feilding, A., Tagliazucchi, E., et al. (2014). The entropic brain: a theory of conscious states informed by neuroimaging research with psychedelic drugs. Front. Hum. Neurosci. 8:20. doi: 10.3389/fnhum.2014.00020

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, W. J., Faraone, S. V., Biederman, J., and Tsuang, M. T. (1994). Diagnostic accuracy of the Child Behavior Checklist scales for attention-deficit hyperactivity disorder: a receiver-operating characteristic analysis. J. Consult. Clin. Psychol. 62, 1017–1025. doi: 10.1037/0022-006X.62.5.1017

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, M. X., Axmacher, N., Lenartz, D., Elger, C. E., Sturm, V., and Schlaepfer, T. E. (2009a). Good vibrations: cross-frequency coupling in the human nucleus accumbens during reward processing. J. Cogn. Neurosci. 21, 875–889. doi: 10.1162/jocn.2009.21062

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, M. X., Elger, C. E., and Fell, J. (2009b). Oscillatory activity and phase-amplitude coupling in the human medial frontal cortex during decision making. J. Cogn. Neurosci. 21, 390–402. doi: 10.1162/jocn.2008.21020

PubMed Abstract | CrossRef Full Text | Google Scholar

Colgin, L. L., Denninger, T., Fyhn, M., Hafting, T., Bonnevie, T., and Jensen, O. (2009). Frequency of gamma oscillations routes flow of information in the hippocampus. Nature 462, 353–357. doi: 10.1038/nature08573

PubMed Abstract | CrossRef Full Text | Google Scholar

Damaraju, E., Allen, E. A., Belger, A., Ford, J. M., Mathalon, D. H., McEwen, S., et al. (2014). Dynamic functional connectivity analysis reveals transient states of dysconnectivity in schizophrenia. Neuroimage Clin. 5, 298–308. doi: 10.1016/j.nicl.2014.07.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Deco, G., and Corbetta, M. (2011). The dynamical balance of the brain at rest. Neuroscientist 17, 107–123. doi: 10.1177/1073858409354384

PubMed Abstract | CrossRef Full Text | Google Scholar

Deco, G., and Jirsa, V. K. (2012). Ongoing cortical activity at rest: criticality, multistability, and ghost attractors. J. Neurosci. 32, 3366–3375. doi: 10.1523/JNEUROSCI.2523-11.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Deco, G., Jirsa, V. K., and McIntosh, A. R. (2011). Emerging concepts for the dynamical organization of resting-state activity in the brain. Nat. Rev. Neurosci. 12, 43–56. doi: 10.1038/nrn2961

PubMed Abstract | CrossRef Full Text | Google Scholar

Deco, G., Jirsa, V. K., and McIntosh, A. R. (2013). Resting brains never rest: computational insights into potential cognitive architectures. Trends Neurosci. 36, 268–274. doi: 10.1016/j.tins.2013.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Delorme, A., and Makeig, S. (2004). EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. J. Neurosci. Methods 134, 9–21. doi: 10.1016/j.jneumeth.2003.10.009

PubMed Abstract | CrossRef Full Text | Google Scholar

De Pasquale, F., Penna, S. D., Snyder, A. Z., Lewis, C., Mantini, D., Marzetti, L., et al. (2010). Temporal dynamics of spontaneous MEG activity in brain networks. Proc. Natl. Acad. Sci. U.S.A. 107, 6040–6045. doi: 10.1073/pnas.0913863107

PubMed Abstract | CrossRef Full Text | Google Scholar

Dimitriadis, S. I., Kanatsouli, K., Laskaris, N. A., Tsirka, V., Vourkas, M., and Micheloyannis, S. (2012b). Surface EEG shows that Functional Segregation via phase coupling contributes to the neural substrate of mental calculations. Brain Cogn. 80, 45–52. doi: 10.1016/j.bandc.2012.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Dimitriadis, S. I., Laskaris, N. A., Bitzidou, M. P., Tarnanas, I., and Tsolaki, M. (2015c). A novel biomarker of amnestic MCI based on dynamic Cross-Frequency Coupling patterns during cognitive brain responses. Front. Neurosci. 9:350. doi: 10.3389/fnins.2015.00350

PubMed Abstract | CrossRef Full Text | Google Scholar

Dimitriadis, S. I., Laskaris, N. A., Del Rio-Portilla, Y., and Koudounis, G. C. (2009). Characterizing dynamic functional connectivity across sleep stages from EEG. Brain Topogr. 22, 119–133. doi: 10.1007/s10548-008-0071-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Dimitriadis, S. I., Laskaris, N. A., and Micheloyannis, S. (2015b). Transition dynamics of EEG-based Network Microstates unmask developmental and task differences during mental arithmetic and resting wakefulness. Cogn. Neurodyn. 9, 371–387. doi: 10.1007/s11571-015-9330-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Dimitriadis, S. I., Laskaris, N. A., Simos, P. G., Micheloyannis, S., Fletcher, J. M., Rezaie, R., et al. (2013b). Altered temporal correlations in resting-state connectivity fluctuations in children with reading difficulties detected via MEG. Neuroimage 83, 307–317. doi: 10.1016/j.neuroimage.2013.06.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Dimitriadis, S. I., Laskaris, N. A., Tsirka, V., Vourkas, M., and Micheloyannis, S. (2010b). What does delta band tell us about cognitive Processes: a mental calculation study? Neurosci. Lett. 483, 11–15. doi: 10.1016/j.neulet.2010.07.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Dimitriadis, S. I., Laskaris, N. A., Tsirka, V., Vourkas, M., and Micheloyannis, S. (2012a). An EEG study of brain connectivity dynamics at the resting state. Nonlinear Dyn. Psychol. Life Sci. 16, 5–22.

PubMed Abstract | Google Scholar

Dimitriadis, S. I., Laskaris, N. A., Tsirka, V., Vourkas, M., Micheloyannis, S., and Fotopoulos, S. (2010a). Tracking brain dynamics via time-dependent network analysis. J. Neurosci. Methods 193, 145–155. doi: 10.1016/j.jneumeth.2010.08.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Dimitriadis, S. I., Laskaris, N. A., and Tzelepi, A. (2013a). On the quantization of time-varying phase synchrony patterns into distinct Functional Connectivity Microstates (FCμstates) in a multi-trial visual ERP paradigm. Brain Topogr. 3, 397–409. doi: 10.1007/s10548-013-0276-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Dimitriadis, S. I., Laskaris, N. A., Tzelepi, A., and Economou, G. (2012c). Analyzing functional brain connectivity by means of commute times: a new approach and its application to track event-related dynamics. IEEE (TBE) Trans. Biomed. Eng. 59, 1302–1309. doi: 10.1109/TBME.2012.2186568

PubMed Abstract | CrossRef Full Text | Google Scholar

Dimitriadis, S. I., Sun, Y. U., Kwok, K., Laskaris, N. A., Thakor, N., and Bezerianos, A. (2015a). Cognitive workload assessment based on the tensorial treatment of EEG estimates of cross-frequency phase interactions. Ann. Biomed. Eng. 43, 977–989. doi: 10.1007/s10439-014-1143-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Dimitriadis, S. I., Zouridakis, G., Rezaie, R., Babajani-Feremi, A., and Papanicolaou, A. C. (2015d). Functional connectivity changes detected with magnetoencephalography after mild traumatic brain injury. Neuroimage Clin. 9, 519–531. doi: 10.1016/j.nicl.2015.09.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Dimitriadis, S., Sun, Y., Laskaris, N., Thakor, N., and Bezerianos, A. (2016). Revealing cross-frequency causal interactions during a mental arithmetic task through symbolic transfer entropy: a novel vector-quantization approach. IEEE Trans. Neural Syst. Rehabil. Eng. doi: 10.1109/TNSRE.2016.2516107. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

Dutta, A., McKie, S., and Deakin, J. F. W. (2014). Resting state networks in major depressive disorder. Psychiatry Res. 224, 139–151. doi: 10.1016/j.pscychresns.2014.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Engel, A. K., Gerloff, C., Hilgetag, C. C., and Nolte, G. (2013). Intrinsic coupling modes: multiscale interactions in ongoing brain activity. Neuron 80, 867–886. doi: 10.1016/j.neuron.2013.09.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Escudero, J., Hornero, R., Abásolo, D., and Fernández, A. (2011). Quantitative evaluation of artifact removal in real magnetoencephalogram signals with blind source separation. Ann. Biomed. Eng. 39, 2274–2286. doi: 10.1007/s10439-011-0312-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Fair, D. A., Cohen, A. L., Dosenbach, N. U., Church, J. A., Miezin, F. M., Barch, D. M., et al. (2008). The maturing architecture of the brain's default network. Proc. Natl. Acad. Sci. U.S.A. 105, 4028–4032. doi: 10.1073/pnas.0800376105

PubMed Abstract | CrossRef Full Text | Google Scholar

Florin, E., and Baillet, S. (2015). The brain's resting-state activity is shaped by synchronized cross-frequency coupling of neural oscillations. Neuroimage 111, 26–35. doi: 10.1016/j.neuroimage.2015.01.054

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedman, J., and Rafsky, L. (1981). Graphics for the multivariate two-sample problem. J. Am. Stat. Assoc. 76, 277–287. doi: 10.1080/01621459.1981.10477643

CrossRef Full Text | Google Scholar

Greicius, M. D., Krasnow, B., Reiss, A. L., and Menon, V. (2003). Functional connectivity in the resting brain: a network analysis of the default mode hypothesis. Proc. Natl Acad. Sci. U.S.A. 100, 253–258. doi: 10.1073/pnas.0135058100

PubMed Abstract | CrossRef Full Text | Google Scholar

Haider, B., Duque, A., Hasenstaub, A. R., and McCormick, D. A. (2006). Neocortical network activity in vivo is generated through a dynamic balance of excitation and inhibition. J. Neurosci. 26, 4535–4545. doi: 10.1523/JNEUROSCI.5297-05.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansen, E. C. A., Battaglia, D., Spiegler, A., Deco, G., and Jirsa, V. K. (2015). Functional connectivity dynamics: modeling the switching behavior of the resting state. Neuroimage 105, 525–535. doi: 10.1016/j.neuroimage.2014.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Heine, L., Soddu, A., Gómez, F., Vanhaudenhuyse, A., Tshibanda, L., Thonnard, M., et al. (2012). Resting state networks and consciousness: alterations of multiple resting state network connectivity in physiological, pharmacological, and pathological consciousness states. Front. Psychol. 3:295. doi: 10.3389/fpsyg.2012.00295.

PubMed Abstract | CrossRef Full Text | Google Scholar

Hipp, J. F., Hawellek, D. J., Corbetta, M., Siegel, M., and Engel, A. K. (2012). Large-scale cortical correlation structure of spontaneous oscillatory activity. Nat. Neurosci. 15, 884–889. doi: 10.1038/nn.3101

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoeft, F., McCandliss, B. D., Black, J. M., Gantman, A., Zakerani, N., Hulme, C., et al. (2011). Neural systems predicting long-term outcome in dyslexia. Proc. Natl. Acad. Sci. U.S.A. 108, 361–368. doi: 10.1073/pnas.1008950108

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoeft, F., Meyler, A., Hernandez, A., Juel, C., Taylor Hill, H., Martindale, J., et al. (2007). Functional and morphometric brain dissociation between dyslexia and reading ability. Proc. Natl. Acad. Sci. U.S.A. 104, 4234–4239. doi: 10.1073/pnas.0609399104

PubMed Abstract | CrossRef Full Text | Google Scholar

Ioannides, A. A., Dimitriadis, S. I., Saridis, G. A., Voultsidou, M., Poghosyan, V., Liu, L., et al. (2012). Source space analysis of event-related dynamic reorganization of brain networks. Comput. Math. Methods Med. 2012:452503 doi: 10.1155/2012/452503

PubMed Abstract | CrossRef Full Text | Google Scholar

Jensen, O., and Colgin, L. L. (2007). Cross-frequency coupling between neuronal oscillations. Trends Cogn. Sci. 11, 267–269. doi: 10.1016/j.tics.2007.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, C., Gao, C., Chen, C., Ma, S., Netra, R., Wang, Y., et al. (2011). A preliminary study of the dysregulation of the resting networks in first-episode medication-naïve adolescent depression. Neurosci. Lett. 503, 105–109. doi: 10.1016/j.neulet.2011.08.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Jirsa, V., and Müller, V. (2013). Cross-frequency coupling in real and virtual brain networks. Front. Comput. Neurosci. 7:78. doi: 10.3389/fncom.2013.00078

PubMed Abstract | CrossRef Full Text | Google Scholar

Khazaee, A., Ebrahimzadeh, A., and Babajani-Feremi, A. (2015). Identifying patients with Alzheimer's disease using resting-state fMRI and graph theory. Clin. Neurophysiol. 126, 2132–2141. doi: 10.1016/j.clinph.2015.02.060

PubMed Abstract | CrossRef Full Text | Google Scholar

Kopell, N. J., Gritton, H. J., Whittington, M. A., and Kramer, M. A. (2014). Beyond the connectome: the dynome. Neuron 83, 1319–1328. doi: 10.1016/j.neuron.2014.08.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Koyama, M., Martino, A., Zuo, X.-N., Kelly, C., Mennes, M., Jutagir, D., et al. (2011). Resting-state functional connectivity indexes reading competence in children and adults. J. Neurosci. 31, 8617–8624. doi: 10.1523/JNEUROSCI.4865-10.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Koyama, M. S., Kelly, C., Shehzad, Z., Penesetti, D., Castellanos, F. X., and Milham, M. P. (2010). Reading networks at rest. Cereb. Cortex 10, 2549–2559. doi: 10.1093/cercor/bhq005

PubMed Abstract | CrossRef Full Text | Google Scholar

Kramer, M. A., Tort, A. B., and Kopell, N. J. (2008). Sharp edge artifacts and spurious coupling in EEG frequency comodulation measures. J. Neurosci. Methods 170, 352–357. doi: 10.1016/j.jneumeth.2008.01.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Lachaux, J. P., Rodriguez, E., Martinerie, J., and Varela, F. J. (1999). Measuring phase synchrony in brain signals. Hum. Brain Mapp. 8, 194–208.

PubMed Abstract | Google Scholar

Meyler, A., Keller, T. A., Cherkassky, V. L., Gabrieli, J. D., and Just, M. A. (2008). Modifying the brain activation of poor readers during sentence comprehension with extended remedial instruction: a longitudinal study of neuroplasticity. Neuropsychologia 46, 2580–2592. doi: 10.1016/j.neuropsychologia.2008.03.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Milne, R. D., Syngeniotis, A., Jackson, G., and Corballis, M. C. (2002). Mixed lateralization of phonological assembly in developmental dyslexia. Neurocase 8, 205–209. doi: 10.1093/neucas/8.3.205

PubMed Abstract | CrossRef Full Text | Google Scholar

Niogi, S. N., and McCandliss, B. D. (2006). Left lateralized white matter microstructure accounts for individual differences in reading ability and disability. Neuropsychologia 44, 2178–2188. doi: 10.1016/j.neuropsychologia.2006.01.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Nolte, G., Bai, O., Wheaton, L., Mari, Z., Vorbach, S., and Hallett, M. (2004). Identifying true brain interaction from EEG data using the imaginary part of coherency. Clin. Neurophysiol. 115, 2292–2307. doi: 10.1016/j.clinph.2004.04.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Northcutt, G., and Kaas, J. H. (1995). The emergence and evolution of mammalian neocortex. Trends Neurosci. 18, 373–379. doi: 10.1016/0166-2236(95)93932-N

PubMed Abstract | CrossRef Full Text | Google Scholar

Oostenveld, R., Fries, P., Maris, E., and Schoelen, J. M. (2011). Fieldtrip: open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data. Comput. Intell. Neurosci. 2011:156869. doi: 10.1155/2011/156869

PubMed Abstract | CrossRef Full Text | Google Scholar

Osipova, D., Hermes, D., and Jensen, O. (2008). Gamma power is phase-locked to posterior alpha activity. PLoS ONE 3:e3990. doi: 10.1371/journal.pone.0003990

PubMed Abstract | CrossRef Full Text | Google Scholar

Özkurt, T.E., and Schnitzler, A. (2011). A critical note on the definition of phase–amplitude cross-frequency coupling. J. Neurosci. Methods 201, 438–443. doi: 10.1016/j.jneumeth.2011.08.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Palva, J. M., and Palva, S. (2011). Roles of multiscale brain activity fluctuations in shaping the variability and dynamics of psychophysical performance. Prog. Brain Res. 193, 335–350. doi: 10.1016/B978-0-444-53839-0.00022-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinneo, L. R. (1966). On noise in the nervous system. Psychol. Rev. 73, 242–247. doi: 10.1037/h0023240

PubMed Abstract | CrossRef Full Text | Google Scholar

Pugh, K. R., Mencl, W. E., Jenner, A. R., Katz, L., Frost, S. J., Lee, J. R., et al. (2001). Neurobiological studies of reading and reading disability. J. Commun. Dis. 34, 479–492. doi: 10.1016/S0021-9924(01)00060-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Raichle, M. E., MacLeod, A. M., Snyder, A. Z., Powers, W. J., Gusnard, D. A., Shulman, G. L., et al. (2001). A default mode of brain function. Proc. Natl. Acad. Sci. U.S.A. 98, 676–682. doi: 10.1073/pnas.98.2.676

PubMed Abstract | CrossRef Full Text | Google Scholar

Rezaie, R., Simos, P. G., Fletcher, J. M., Cirino, P., Vaughn, S., and Papanicolaou, A. C. (2011). Temporoparietal brain activity as a longitudinal predictor of response to educational interventions among middle school struggling readers. J. Int. Neuropsychol. Soc. 17, 875–885. doi: 10.1017/S1355617711000890

PubMed Abstract | CrossRef Full Text | Google Scholar

Richards, T. L., Grabowski, T. J., Boord, P., Yagle, K., Askren, M., Mestre, Z., et al. (2015). Contrasting brain patterns of writing-related DTI parameters, fMRI connectivity, and DTI–fMRI connectivity correlations in children with and without dysgraphia or dyslexia. Neuroimage Clin. 8, 408–421. doi: 10.1016/j.nicl.2015.03.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Richman, J. S., and Moorman, J. R. (2000). Physiological time-series analysis using approximate entropy and sample entropy. Am. J. Physiol. Heart Circul. Physiol. 278, H2039–H2049.

PubMed Abstract | Google Scholar

Rigas, P., Adamos, D. A., Sigalas, C., Tsakanikas, P., Laskaris, N. A., and Skaliora, I. (2015). Spontaneous Up states in vitro: a single-metric index of the functional maturation and regional differentiation of the cerebral cortex. Front. Neural Circuits 9:59. doi: 10.3389/fncir.2015.00059

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosazza, C., and Minati, L. (2011). Resting-state brain networks: literature review and clinical application. Neurol. Sci. 32, 773–785. doi: 10.1007/s10072-011-0636-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Shanahan, M. (2010). Embodiment and the Inner Life: Cognition and Consciousness in the Space of Possible Minds. Oxford: Oxford University Press.

Google Scholar

Shew, W. L., Yang, H., Petermann, T., Roy, R., and Plenz, D. (2009). Neuronal avalanches imply maximum dynamic range in cortical networks at criticality. J. Neurosci. 29, 15595–15600. doi: 10.1523/jneurosci.3864-09.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Simos, P. G., Fletcher, J. M., Sarkari, S., Billingsley, R. L., Denton, C., and Papanicolaou, A. C. (2007). Altering the brain circuits for reading through intervention: a magnetic source imaging study. Neuropsychology 21, 485–496. doi: 10.1037/0894-4105.21.4.485

PubMed Abstract | CrossRef Full Text | Google Scholar

Simos, P. G., Rezaie, R., Fletcher, J. M., Juranek, J., Passaro, A. D., Li, Z., et al. (2011). Functional disruption of the brain mechanism for reading: effects of comorbidity and task difficulty among children with developmental learning problems. Neuropsychology 24, 520–534. doi: 10.1037/a0022550

PubMed Abstract | CrossRef Full Text | Google Scholar

Sporns, O. (2014). Contributions and challenges for network models in cognitive neuroscience. Nat. Neurosci. 17, 652–660. doi: 10.1038/nn.3690

PubMed Abstract | CrossRef Full Text | Google Scholar

Sporns, O., Chialvo, D. R., Kaiser, M., and Hilgetag, C. C. (2004). Organization, development and function of complex brain networks. Trends Cogn. Sci. 8, 418–425. doi: 10.1016/j.tics.2004.07.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Stam, C. J., de Haan, W., Daffertshofer, A., Jones, B. F., Manshanden, I., vanCappellen van Walsum, A. M., et al. (2009). Graph theoretical analysis of magnetoencephalographic functional connectivity in Alzheimer's disease. Brain 132, 213–224. doi: 10.1093/brain/awn262

PubMed Abstract | CrossRef Full Text | Google Scholar

Steriade, M., Contreras, D., Amzica, F., and Timofeev, I. (1996). Synchronization of fast (30–40Hz) spontaneous oscillations in intrathalamic and thalamocortical networks. J. Neurosci. 16, 2788–2808.

PubMed Abstract | Google Scholar

Swanson, J., Schuck, S., Mann, M., Carlson, C., Hartman, K., Sergeant, J., et al. (2005). Categorical and dimensional definitions and evaluations of symptoms of ADHD. History of the SNAP and the SWAN Ratings Scales. Inter. J. Edu. Psychol. Assess. 10, 51–70.

Google Scholar

Tagliazucchi, E., Balenzuela, P., Fraiman, D., and Chialvo, D. R. (2012). Criticality in large-scale brain FMRI dynamics unveiled by a novel point process analysis. Front. Physiol. 3:15. doi: 10.3389/fphys.2012.00015

PubMed Abstract | CrossRef Full Text | Google Scholar

Theiler, J., Eubank, S., Longtin, A., Galdrikian, B., and Farmer, J. D. (1992). Testing for nonlineaity in time series:the method of surrogate data. Physica D 85:77. doi: 10.1016/0167-2789(92)90102-S

CrossRef Full Text | Google Scholar

Tognoli, E., and Kelso, J. A. S. (2014a). Enlarging the scope: grasping brain complexity. Front. Syst. Neurosci. 8:122. doi: 10.3389/fnsys.2014.00122

PubMed Abstract | CrossRef Full Text | Google Scholar

Tognoli, E., and Kelso, J. A. S. (2014b). The metastable brain. Neuron 81, 35–48. doi: 10.1016/j.neuron.2013.12.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Tononi, G., Sporns, O., and Edelman, G. M. (1994). A measure for brain complexity: relating functional segregation and integration in the nervous system. Proc. Natl. Acad. Sci. U.S.A. 91, 5033–5037. doi: 10.1073/pnas.91.11.5033

PubMed Abstract | CrossRef Full Text | Google Scholar

Tort, A. B. L., Komorowski, R., Eichenbaum, H., and Kopell, N. (2010). Measuring phase-amplitude coupling between neuronal oscillations of different frequencies. J. Neurophysiol. 104, 1195–1210. doi: 10.1152/jn.00106.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Varela, F., Lachaux, J. P., Rodriguez, E., and Martinerie, J. (2001). The brainweb: phase synchronization and large-scale integration. Nat. Rev. Neurosci. 2, 229–239. doi: 10.1038/35067550

PubMed Abstract | CrossRef Full Text | Google Scholar

Vaughn, S., Wanzek, J., Wexler, J., Barth, A., Cirino, P. T., Fletcher, J. M., et al. (2010). The relative effects of group size on reading progress of older students with reading difficulties. Read. Writ. 23, 931–956. doi: 10.1007/s11145-009-9183-9

PubMed Abstract | CrossRef Full Text | Google Scholar

von Stein, A., and Sarnthein, J. (2000). Different frequencies for different scales of cortical integration: from local gamma to long range alpha/theta synchronization. Int. J. Psychophysiol. 38, 301–313. doi: 10.1016/S0167-8760(00)00172-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Voytek, B., Canolty, R. T., Shestyuk, A., Crone, N. E., Parvizi, J., and Knight, R. T. (2010). Shifts in gamma phase– amplitude coupling frequency from theta to alpha over posterior cortex during visual tasks. Front. Hum. Neurosci. 4:191. doi: 10.3389/fnhum.2010.00191

PubMed Abstract | CrossRef Full Text | Google Scholar

Wechsler, D. (1999). Wechsler Abbreviated Scale of Intelligence (WASI). San Antonio, TX: The Psychological Corporation

Google Scholar

Woodcock, R. W., McGrew, K. S., and Mather, N. (2001). WJ-III Tests of Achievement. Itasca, IL:Riverside Publishing.

Google Scholar