The data we used in this study were collected as part of a broader project aiming to identify the neural underpinnings of dyslexia in children. In this section, we briefly describe the key parameters of the RAN experimental task and the data collection procedure relevant to our analysis; we refer to Christoforou et al. (2021a) for full details on the data collection apparatus.

A particular challenge in analyzing EEG signals obtained during the serial RAN is the lack of experimenter-controlled, time-locked trials necessary to extract Event-related Potentials. As such, we opted to explore the neural activity time-locked to the onset of eye fixations; this activity is referred to as single-trial Fixation-related Potential (sFRP). To extract the sFRP, we integrate information from eye-tracking and EEG measurements. In particular, the onset of each eye-fixation’s timestamp in the dataset is used as a temporal marker to epoch the EEG segments. More precisely, given an EEG data segment EEG_p,c (i.e., the segment generated during the EEG pre-processing step above, which represents the EEG response of participant p while reading the entire letter-matrix of condition c), and given the corresponding set of fixations FIX_p,c, we epoch the EEG_p,c between −200 ms to 500 ms of the onset time of each f ∈ FIX_p,c and subtract the baseline amplitude of each epoch. This procedure results in a new set defined as

where sFRP_i corresponds to the EEG epoch on the onset time of the i-th fixation of FIX_p,c. We note that the two sets have the same cardinality. The generation of all sFRPs was implemented using custom python code.

Our objective was to isolate neural components in the extracted fixation-related potentials that were likely modulated by RAN tasks and were informative of differences between CAC and DYS. Our approach was motivated by the hypothesis that the neural activity of participants that have developed adequate reading skills would exhibit neural activation patterns congruent with other participants with adequately developed reading skills. While contrarily, participants who experienced reading difficulties would have neural responses that deviated from such stereotypical patterns. Toward this objective, we formulated an optimization procedure to isolate neural components congruent among participants with sufficiently developed reading skills and explore those components as potential differentiation neuro-markers between CAC and DYS. Here, we provide details of our approach to isolate such reading-related neural-congruency components.

We seek to identify components (i.e., spatial projections of the fixation-related potential) that capture neural activity that maximally correlates among a group of children with adequately developed reading skills (i.e., CAC group). For this, we formulate an optimization problem as follows: for a group of S participants, 𝒮 = {s₁, s₂, …, s_S}, where s_i ∈ ℤ⁺ denotes a participants index, representing a CAC group, we define the between-subject and within-subject cross-covariance matrices was :

where

where K is a normalizing scalar, δ_ij is the Kronecker delta¹, and X_s ∈ ℝ^{D × S.F} is the horizontally concatenated matrix comprised of the fixation-related potential of a participant s during a given condition (i.e., reading of a letter-matrix stimulus), defined as:

For a spatial projection vector w ∈ ℝ^D, the average Pearson Product Moment Correlation Coefficient between the fixation-related potentials, projected onto vector w, across every pair of participants in the group is then defined as:

We consider ρ as a measure of the degree of congruency in reading-related neural activity (projected onto component w) among participants with adequately developed reading skills. As such, we seek to find the component w that maximized ρ. Taking the derivative of ρ with respect to w and setting it to zero, we get the solution of the optimization given as the eigenvectors to the generalized eigenvalue problem:

where w_k is the k-th eigenvector of the matrix (Rw-1⁢Rb) and corresponds to the component (i.e., spatial projection vector) that captures the k-th largest correlation in neural activity, and λ_k is the corresponding eigenvalue and denotes the strength of the correlation. We note that since (Rw-1⁢Rb) is a D × D matrix (D being the number of channels), there are D solutions to the optimization problem (i.e. the D eigenvectors), each identifying a component at different correlation strength, with the first eigenvector (i.e., k = 1) having the strongest correlation, and subsequent components appearing in descending order of correlation strength. As such, the vector w₁defines a component (spatial projection) where neural activity is most strongly correlated among participants in the adequately developed reading skills, w₂ defines the component where neural-activity exhibits the second strongest correlation among the groups, and so on.

To determine the reading-related neural-congruency (RRNC) score of an individual participant s with respect to the kth component, we measure the correlation between the fixation-related potentials of the subject to the fixation-related potentials of each subject in the group, after projecting both onto the component w_k. Formally, we define the reading-related neural congruency score for a participant s and a component w_k as:

where

We calculate RRNC scores separately for each condition (i.e., Rime-confound, Rime-non-confound, Visual-confound, Visual-non-confound). Moreover, to avoid training bias during the component extractions, data from the subject to be tested were excluded from the component extractions step. Finally, we define the neural-congruency feature vector of each participant s as:

Each RRNC_s,k score measures the strength of the congruency in neural activity between participant s and the CAC group, with respect to the k-th component. Therefore, the neural-congruency feature vector c_s encapsulates the neural-congruency of participant s, across all K components. In essence, the feature vector c_s characterizes the overall congruency observed in the participant’s neural activity for each extracted component. We note that the dimensions of the vector c_s are indexed in descending order, according to the lambda score of each extracted component.

Our goal was to explore whether information captured in the feature vector of neural-congruency components c_s is predictive of differentiating between DYS and CAC groups during the RAN task. Toward this objective, we considered two approaches for aggregating the RRNC scores into determining markers of dyslexia. The two approaches are detailed below.

Given the solutions to the generalized eigenvalue problem, the temporal profile of each component was calculated as the product of each component w^k, with each single-trial response and taking the grant-average response of the projected components. Moreover, the topographical profile (i.e., the forward model) of each component was calculated as follows:

The forward model captures the covariance between each component’s activity as measured by each electrode.

To avoid training bias, all model parameters, including the extracted neural-congruency components and classifier weights, are trained using a leave-one-participant-out cross-validation procedure. The classifier’s generalization performance is calculated as the area under the Receiver Operator Characteristic curve (AUC) on cross-validated scores. A permutation test is used to determine statistical significance levels over AUC scores (10,000 repetitions). A two-way ANOVA was used for group and grade comparisons, with the neural congruency metrics as the dependent variable.

For each participant, Cumulative Neural-Congruency scores were calculated as the sum of the three components with the highest lambda values, corresponding to the components whose projection has the highest correlation. To avoid training biases during the neural-congruency component identification, data from the participant for whom the Neural-Congruency score was to be calculated was excluded from the component identification step. Neural-Congruency scores were obtained and analyzed separately for each of the four conditions (i.e., Rime-confusable; Rime Non-confusable; Visual-confusable; Visual Non-confusable). For each condition, a separate two-way ANOVA was performed to analyze the effect of group (i.e., CAC vs. DYS) and grade (grade 3 vs grade 6) on the Neural-Congruency scores. A two-way ANOVA on Rime-confusable Cumulative RAN-related Neural-Congruency scores shows a significant main effect of grade (p < 0.04), with grade 6 group showing higher neural-congruency than the grade 3 group. The analysis also releveled there was not a significant interaction effect between the group and grade, F(1,52) = 0.740, p = 0.39, nor a statistically significant effect for group (p = 0.19). A two-way ANOVA on Rime non-confusable Neural-Congruency scores revealed a significant main effect of group (p < 0.01). Participants in the control group showed a higher Cumulative RAN-Related Neural-Congruency scores than the participants with dyslexia group. There was no main effect of grade (p = 0.61), and there was no statistically significant effect observed between group and grade F(1,52) = 1.67, p = 0.20. On the Visual Confusable task, a two-way ANOVA revealed a statistically significant interaction effect between the group and grade F(1,52) = 4.22, p = 0.04. The analysis also revealed there was not a statistically significant main effect of either the group or the grade (p > 0.05). Finally, the Visual non-confusable task revealed a main effect of the group (p < 0.01), with participants in the control group showing a higher Neural-congruency scores than the children with dyslexia group (DYS). No grade or interaction effect between grade and groups were observed during the Visual non-confusable. Figure 2, shows the box plots for the four two-way ANOVA for each condition. Moreover, a two-factor repeated measures ANOVA was performed to compare the effect of modality condition (rime vs visual) and confusability (confusable vs non-confusable). The analysis showed a statistically significant difference in RAN-related Neural-Congruency scores between modality conditions (p < 0.001) with the Rime modality exhibiting higher scores.

We aimed to further explore the characteristics of the underlying neural activity captured by the neural-congruency components and determine the degree to which each neural-congruency component contributes toward inferring an individual’s group (i.e., DYS or CAC). We hypothesized that a weighted aggregation of individual neural-congruency components would carry predictive information about the participant’s condition. We employed a sparse logistic regression classifier using the ten components with the highest lambda values as independent variables, and an individual’s group (DYS or CAC) as the dependent variable. The classifier was modeled and trained according to the procedure described in Section “2.5.2. LASSO-weighted RAN-related neural-congruency metric.” The statistical significance levels over AUC scores were established using a permutation test (10,000 repetitions). In all four conditions, the cross-validated AUC scores of the classifiers show that the LASSO-weighted RAN-related Neural-congruency metric predicted an individual’s group. The prediction accuracy for all conditions, the Rime-confusable condition (AUC = 0.86, p < 0.00001), Rime Non-confusable condition (AUC = 0.86, p < 0.00001) the Visual Non-confusable (AUC = 0.81, p < 0.00001), and the Visual Confusable condition (AUC = 0.73, p < 0.005) was high and statistically significant. The ROC curves and corresponding AUC scores for each condition and the 95th-percentile envelop of the ROC curve under the null-distribution are depicted in Figure 3. The boxplot in Figure 4 shows the distribution of the lasso-weighted RAN-related neural-congruency for each grade and group and experimental condition.

The spatio-temporal profile of each RAN-related Neural Congruency for all conditions is shown in Figure 5 and the Supplementary material. Each spatio-temporal profile comprises the “Forward model,” which shows the spatial distribution of the correlated neural activity captured by the corresponding component, and the temporal profile–the time course of the FRP’s neural activity when projected onto that neural-congruency component. Visual inspection of the temporal profile provides insights into timeframe differences between groups and condition intensify. Similarly, visual inspection of the forward model alludes to potential brain areas from which the underlying neural activity originated from. The weights associated with each channel in the forward model capture the electrical coupling of the correlated components. The components are ordered based on their corresponding lambda scores, with component #1 reflecting the highest lambda, and component #10 the smallest lambda value. The LASSO-coefficients associated with each component are shown over each component’s forward model and indicate the weight used to aggregate the neural-congruency components.

Analysis of the behavioral data obtained during the experiment has been previously reported by Christoforou et al. (2021a) and it is outside the scope of this papers. However, to provide the context on our results on neurophysiological data, we include a summary of the behavioral data analysis of this experiment. A MANOVA analysis was performed on the behavioral data, with the naming speed performance time for each of the four RAN tasks as dependent measures and Group (2) as a fixed factor. The main group effects was significant, Wilks’ Lambda = 0.754, F(4,55) = 4.48, p < 0.01, η² = 0.20. Subsequent univariate analyses demonstrated that the group’s main effect was significant for all individual measures after Bonferroni adjustments (Supplementary Table 1). The DYS group performed significantly poorer than the chronological age controls in all naming speed measures.

Regarding the contribution of the Cumulative RAN-related neural-congruency components in differentiating between children with and without dyslexia, the results revealed significant differences between the groups in the both non-confusable conditions (i.e., Rime non-confusable, and Visual non-confusable). Specifically, on the one hand, cumulative RAN-related neural-congruency scores in typically developing children were significantly higher than their counterparts in the DYS groups. This finding denotes an increase in the synchronicity in neural activity in the CAC groups, which suggests that the CAC group has developed a more consistent stereotypical response in neural activations when engaging processes associated with the execution of the serial RAN task. Breznitz and colleagues (Breznitz, 2001, 2003; Breznitz and Misra, 2003) have proposed the ‘synchronization hypothesis’ to describe this phenomenon. According to this hypothesis, accurate information integration in decoding words can occur only when the modalities and brain systems are synchronized. This synchronization, therefore, requires that the processing speed and the accuracy with which content information is processed and transferred within and between the various activated neural systems are readily available. Our findings further confirm this hypothesis that this synchronicity in neural activity is evident in typically developing readers but to a lesser degree in children with dyslexia.

Indeed, on the other hand, the consistency in the synchronicity of the responses diminishes across the DYS groups, suggesting a lack of regularity in the processing the letter stimuli in the serial RAN. Visual inspection of the time course of the three neural-congruency components used in calculating the Cumulative RAN-related neural-congruency score suggests that the difference in congruency appears between 100 ms–200 ms following the fixation onset. The analysis also shows that the effect that emerged in the non-confusable tasks is not present in the confusable tasks (i.e., Rime-confusable, and Visual-confusable). That is, there were no significant differences observed in the Cumulative scores between DYS and CAC. We interpret these effects in the context of neural efficiency theory. Specifically, we argue that the CAC group has developed efficient mechanisms for recognizing, decoding, and reading letters as captured by the consistency in the neural responses across participants. In contrast, the DYS group does not show the same regularity in neural responses, suggesting the corresponding mechanisms for recognizing, decoding, and reading letters are less fine-tuned in the DYS group. This latter finding has additive value to Breznitz’s Asynchrony Theory (Breznitz, 2008) which proposes that dyslexia is an outcome of the failure to synchronize the various brain entities activated during reading-related processes. Nevertheless, this asynchrony is also evident for children with dyslexia and their typically developing counterparts when the stimulus’ complexity increases – as occur in the rime or visual confusability. This finding underscores the need for additional neural resources to resolve the stimulus’s confounding elements (rime or visual). To that end, the Cumulative RAN-related neural congruency components, fail to capture a consistent neural-response across either of the groups.

Although the simple aggregation of the top three neural-congruency scores revealed significant differences between groups under the non-confounding conditions, we hypothesized that an optimally weighted sum over all ten neural-congruency components would capture additional differences in neural activations between groups. Indeed, the sparse logistic regression classification revealed that an optimally weighted aggregation over the neural-congruency components differentiates between CAC and DYS across both confounding and non-confounding conditions. These results suggest that components beyond the three with the highest lambda values do capture neural activity relevant to the task. Moreover, the weights associated with each component differ in both amplitude and sign (i.e., they can contribute either positively or negatively to the sum). This finding suggests that the neural activations captured by each individual component might appear with different intensity and polarity in each group. For example, one component that might capture activity associated with character-disambiguation might exhibit stronger synchronicity in one group (i.e., contributing positively to the sum); however, another component might exhibit stronger synchronicity in the other group. Thus, simple aggregation of components might result in cancelling out this effect. Overall, our classification model suggests that individual components must be considered when analyzing neural activations.

Spatial (i.e., the forward model) and temporal profiles of the extracted components are depicted in Figure 4 for each condition. The forward model of the neural-congruency component #1 (i.e., the one with the highest lambda score) exhibits a similar topography across all four conditions; moreover, their corresponding temporal profiles show the neural activity is most strongly modulated at a time window of around 100 ms. The similarity suggests that component #1 captures neural activity common to all four conditions, albeit at different intensity levels. At the very least, this finding confirms previous evidence showing that processing complex features of textual stimuli is reflected in the electrophysiological responses around 100 ms after stimulus presentation (Hauk et al., 2006).

The temporal profile of components #10 shows similarity in waveform across the four conditions and a peak amplitude at around 100 ms following the fixations onset. Similarities in the spatial profiles are observed among some of the remaining components as well, although the indexing/ordering of those components varies among conditions. Such similarities suggest that the matching neural-congruency components likely capture neural activity originating from the same underlying source. The variation in the indexing is expected since the ranking of the components is established independently for each condition and depends on the relative strength of all the neural-congruency components in that condition. Interestingly, several projected temporal profiles in each condition display a stereotypical response consistent with neural activations often observed in traditional grant average Event-related Potential analysis (i.e., N/P100 and N170). For example, the temporal profile of component #9 shows an N170 waveform response with visible differences in amplitude at around 170 ms between CAC and DYS. Moreover, the forward model topography of this component shows it to emerge more strongly in electrodes over the left posterior-occipital regions. In the literature, the difference in N170 over the left-posterior occipital region is regarded as an electrophysiological marker of visual expertise (Varga et al., 2020). Particularly, children with a lower letter knowledge, as pre-readers, have shown reduced N170 amplitudes and delayed N170 latency compared to typical readers during letter-string presentations (Maurer et al., 2005). Therefore, component #9, extracted by our method, can be interpreted as potentially indicating a visual precursor to literacy resulting from familiarity with letter strings. Furthermore, negativity components at 170 ms have been associated with attention modulation (i.e., Kropotov, 2016), and in turn, attention represents a known latent common cognitive factor of RAN and reading (Papadopoulos et al., 2016). Therefore, part of the neural activity captured by the component could also reflect those distal processes. Moreover, several forward model scalp plots display topographies that often arise as the “forward problem” solutions to single-source dipole models, indicating that those components likely capture source-localized neural activity of different neural processes. Finally, components with high lasso coefficients (i.e., either positive or negative) capture underlying neural activity that more strongly differentiate between DYS and CAC. Thus, group differences among groups appear in the underlying sources modeled by those components.

Further studies could investigate these components in more detail and draw additional conclusions about the neural and cognitive processes that contribute to these differentiations. For example, an interesting expansion of these findings relates to examining the validity of the suggested neural-congruency components analysis against other more conventional EEG analyses. Given that the literature has only recently started to investigate the important information that the FRPs can provide to studying the RAN-reading relationship or other similar correlates of reading performance (e.g., Christoforou et al., 2021b; Fella et al., 2022), the present findings are considered a promising beginning of this quest. Another interesting expansion would be to examine the relationship of the extracted neural congruency components to eye-tracking-based metrics during RAN, such as the recently proposed entropy-based gaze time-series analysis on RAN (Wang et al., 2022).

In conclusion, the RAN-related neural congruency component, identified by our proposed method, carry information on the neural basis of naming speed that differentiates between children with dyslexia and their typically developing counterparts. The topographies of the resulting components suggest that each component likely captures source-localized neural activity corresponding to distinct neural processes. Moreover, neural differences appear to be distributed across several RAN-related Neural-congruency components but at different intensity levels. Therefore, optimally combining the RAN-related components using machine learning enhanced the power of analysis in identifying differences in both the confusable and non-confusable conditions, which have been missed by simple aggregation of the RAN-related Components. Our findings also support the Neural-congruency hypothesis (Christoforou et al., 2021b; Christoforou and Theodorou, 2021), indicating that neural activity elicited during cognitive tasks is more congruent among participants that have mastered the cognitive skills but less congruent otherwise. Finally, our proposed approach opens up new research directions in studying the neural underpinnings of naming speed and their association with reading performance and reading difficulties. For instance, until recently, evidence concluded that individuals with reading difficulties show increased reliance on inferior frontal regions of the reading network and right-hemisphere posterior regions (e.g., Richlan et al., 2011; Norton et al., 2015). The present findings show that several brain areas contribute to the execution of naming tasks which show similarities with reading tasks. Although simple at the surface level, RAN tasks are multi-componential, as is reading (Papadopoulos et al., 2016). Thus, we believe that the present method succeeded in better defining the properties of other processes (including those of reading) that RAN carries and how these are critical in determining naming speed’s influential role on reading performance. At the very least, capturing the RAN performance in the form of neural components helps us better understand the process involved in performing RAN tasks and explore some reasons for poor performance. Next, provided that the generated components will be further systematically tested against behavioral performance measures, the likelihood of deciphering issues relevant to the significant similarity of the RAN and reading measures, such as seriality, is possible (see Altani et al., 2020). Thus, in future work, we plan to further explore the spatiotemporal characteristics and brain sources of the RAN-related Neural-congruency components and the relationship between the neural underpinnings of naming speed - as captured by the Neural-congruency components- and reading difficulties, as well ex explore eye-tr.