2.1.3. Procedure

Upon arrival at the research center, each child's active assent and signed parental consent was obtained according to protocols approved by the Research Ethics Boards of all participating institutions. Upon parental completion of the in-take assessments, children were tested in a sound-proof electromagnetically shielded EEG booth. Each child was seated at a distance of 58 cm from a 19" flat screen monitor. Children were instructed to respond by pressing a button when they saw a duck and to refrain from pressing the button if they saw any other image.

The entire session consisted of several practice blocks followed by the experimental block. In each block, the duck was shown 25% of the time, while the turtle was shown 75% of the time. Participants proceeded to start the experiment once that they attained 100% accuracy on three consecutive practice blocks. In the experimental block, there were 37 target images and 113 distractor images for a total of 150 trials.

#### 2.1.4. Data Collection and Processing

An elastic cap (Quik Cap, Neuroscan, El Paso, TX, USA) adhering to the standard ten–twenty international system of electrode placement, with recessed Ag–AgCl electrodes (10 mm each) was used for the EEG recordings. The cap montage included 12 electrodes corresponding to electrode reference points at Frontal (F3, Fz, F4), Central (Cz), Temporal (T7, T8), Parietal (P7, Pz, P8), and Occipital (O1, Oz, O2) sites. Horizontal and vertical electro-oculograms (EOG) recorded eye movements with two split bipolar electrodes positioned at the outer canthi for the horizontal EOG and on the suborbital ridge of each eye for the vertical EOG. Previous work and pilot studies in children of similar ages [43–47] sugges<sup>t</sup> no critical loss of reliability in source analysis performed with the present set up.

All impedances were kept below 5 kOhms. Low-pass and high-pass filtering (0.5 to 250 Hz) were applied to the signal prior to digitization. Trials from non-ocular electrode sites that were contaminated by excessive peak-to-peak deflection (i.e., >100 μV or <−100 μV) due to non-stereotypical noise were manually excluded. Brain Electric Source Analysis (BESA v.5.4.28; http://www.besa.de/), an electroencephalographic analysis software package was used to analyze EEG recordings and calculate ERP averages for each of the twelve electrode locations (F3, Fz, F4, Cz, T7, T8, P7, Pz, P8, O1, Oz, O2). Ocular correction was performed using the integrated BESA *adaptive artifact correction* [48] and the *surrogate model* [49]. Principal component analysis decomposition was used to correct for ocular artefacts by selecting components for the horizontal and vertical eye movements as well as eye-blinks. The proportion of rejected trials was less than 15% after artefact correction and removal. ERPs were averaged o ffline separately for each stimulus type (i.e., target and distractor) at each electrode with all epochs time-locked to the onset of the image and ending at 1000 ms after the onset of the image.

#### 2.1.5. ERP Data Analysis

One known issue in customary peak analysis techniques is that the assumed correspondences are based on subjective selections of time windows which are made post hoc, after examining the data [50]; this issue is exacerbated when young children are compared to adults since analyses often rely on many untested a-priori assumptions on alleged morphological correspondences between peaks of ERPs of young children and the adults. An alternative is to standardize time windows by latency and to examine the entire single ERP waveform across the entire epoch by using automatic blinded analysis procedures such as *waveform binning* [40,51–53]. We adopted a hybrid approach in which automatic blinded binning includes traditional peak analysis.

EEGLab [54] was used to analyze EEG recordings and calculate ERP averages for each data point for the twelve electrode locations (F3, Fz, F4, Cz, T7, T8, P7, Pz, P8, O1, Oz, O2). ERPs were averaged offline separately for each stimulus type (i.e., targets and distractors). Each epoch ranged from 200 ms prior to the presentation of the image to 1000 ms preceding the presentation of the image.

EEG sampling rate was 1000 Hz, resulting in a total of 1200 data points for each epoch. Time interval analysis of data points obtained from ERP data were simplified using a standard binning procedure that divided each epoch into bins of equal time intervals. Each epoch was divided into twelve bins of 100 ms time interval. To ensure equivalent data density across bins, the ERP averages across subjects

were transformed to vincentized quantile bins [55–57]. The computation of the quantile bins was performed assuming the following empirical distribution for the time series of the averaged ERPs:

$$F\_1^{-1}(\alpha) = \inf \{ -200 < t < +1000; \ F\_1(t) \ge \alpha \}, \text{ with } 0 < \alpha < 1\tag{1}$$

The Vincent average of the Fis mid-points was then computed as:

$$\sum w\_i F\_i^{-1}(a), \text{ where } i = 1 \dots n \text{, and } w\_1 + \dots + w\_n = 1 \tag{2}$$

Due to its shape-preservation property, this procedure minimizes the effects of distortions of individual differences in the distributions of averaged ERP peak mid-points, and given that the skewness in our data was modest, it partially offset biases in parametric testing associated with our relatively small sample size (see [57]). Amplitudes sorted within each bin for each epoch were assessed for each of the twelve electrode sites and averaged separately for target and distractor stimuli. To draw a comparison between ERPs for target and distractor stimuli over time, the number of bins for distractors was kept consistent with the target even though distractors did not require a response.

Vincentized amplitudes were analyzed in a mixed model design with interval bins (12 levels) as a between subjects factor and electrode (12 levels) and condition: (2 levels: target vs. distractor) as within subjects factors. To compare differences between target and distractor, focused ANOVA contrasts between mean amplitudes for each electrode and time interval were conducted to obtain the minimum significant standardized absolute difference using the following formula:

$$\text{Amp } \text{diff}\_{\mu V} = t\_{\text{crit}} \left( \sqrt{\text{MSE}\_{\text{within}} \Sigma (\lambda\_i^{'} / n\_j)} \text{I} \right) \tag{3}$$

where Abs Amp diffμ<sup>V</sup> indicates the significant difference in amplitude peaks between target and distractor conditions for a given electrode within a given interval bin, *t*crit represents the *t*-value corresponding to the critical *p*-value for determining significance threshold after using a Bonferroni correction for 12 electrodes for multiple comparisons. MSEwithin represents the error factor for focused *t*-contrasts across all comparisons. This was the within subjects Mean Square error for the interaction between electrode × condition from the omnibus ANOVA. λi represents the contrast weights [58].

#### 2.1.6. ERP Activity Paths

The differences in ERP waveforms between target and distractor conditions can be more easily interpreted via ERP activity path analysis [40,53]. This graphic summary method illustrates the temporal sequence of neural activity between electrode regions while considering differences between target and distractor conditions.

To obtain activity paths, we used the following procedure. All significant average ERP differences for target and distractor conditions observed were compared. The electrode corresponding to the maximum average ERP difference between target and distractor conditions at each interval bin was noted and the brain area corresponding to this electrode was plotted on head maps for target or distractor conditions at each interval bin. If similar (non-significantly different) maximum average ERP differences were observed for two or three electrodes within the same interval bin, areas corresponding to both or all three electrodes were assumed to be simultaneously activated for that interval and were illustrated on the head map.

The principle of *neural wiring minimization* [59] was used as a rationale for plotting the least-path sequence between simultaneous ERP activity electrodes/areas within same interval bin. Accordingly, the following restrictions were applied due to neuroanatomical boundaries: neural activation from occipital electrodes (O1, Oz, O2) could only move forward or laterally, neural activation from frontal electrodes (F3, Fz, F4) could only move backward or laterally. Likewise, lateral activation from left electrodes (F3, T7, P7, O1) could only occur toward the right, whereas lateral activation from right electrodes (F4, T8, P8, O2) could only occur toward the left.

#### *2.2. Part 2: Modeling of Children's ERP Data*

In the second part of the study, we run a type of Independent Component Analysis (ICA; [60]) on which we built models of the children's data. Figure 2 describes the steps in the analysis pipeline. ICA (Step 1) was used to build ERP topographical mappings (Step 2), and then to identify the dipoles corresponding to target and distractor ERPs; this information was then mapped on a pediatric structural MRI template from the Talairach coordinate system (Step 3). Using corresponding neuroanatomical labels and a visual matching procedure, we confirmed and translated the pediatric Talairach coordinates into an adult structural MRI template (Step 4).

#### 2.2.1. Independent Component Analysis (Step 1) and Topographic Mapping (Step 2)

A second Independent Component Analysis (ICA) was performed to identify ERP components using the EEGLab FASTICA algorithm [61]. While the ICA method can estimate location and timings of components, it cannot estimate an absolute magnitude for them since there is an inherent ambiguity between the strength of the component and the attenuation from it to the measurement point. Therefore, the results were converted to topographic maps.

EEGLab includes an editing graphic utility which allows to represent the epochs of averaged ERP onto topographic maps as clips of 10 ms, and then it permits to put these together in sequence resulting in a capture of time course of the dynamic ERP activity. The specific single topographic maps selected for further analysis reflected scalp activity at the mid-points of the time intervals corresponding with vincentized bins corresponding to those described in the ERP activity paths.

#### 2.2.2. Dipole Analysis (Step 3)

The DIPFIT component of EEGLab was used to estimate a set of dipoles in the averaged ERP data that would explain the independent components extracted. Each dipole is assumed to be a region of cortex where several thousand neurons act together in parallel so that their combined electric field is responsible for the EEG signal measured at the scalp. The DIPFIT software usually finds one or sometimes two dipoles for each of the specific regions that appear to produce the independent components.

The EEGLab MRI-based spherical head model with standard age-appropriate (pediatric) Talairach coordinates was selected. The labels of the brain regions which the locations corresponded to were found using the most recently updated version of the Talairach database [62]. We used the built-in function of this software which permits searching for the nearest grey matter within concentric cubes (voxels) from a minimum of ±1 mm up to a maximum range within ±5 mm to exact dipole origin. That is, nearest gray searches involve concentric cube searches with varying diameters. In general, it searches consecutively larger cubes until it finds a gray matter label, with the same outer limit of a 11 mm-wide cube, so it is also possible to find no gray matter labels.

#### 2.2.3. Translation to Adult MRI Template (Step 4)

The pediatric MRI coordinates obtained through the Talairach database [62] were converted to the Yale Bioimage Suite [63] by entering the pediatric coordinates and by matching the anatomical labels. This process was confirmed by visual analysis based on consensus among two independent anonymous judges with expertise in pediatric and adult neuroanatomy.

#### *2.3. Part 3: Simulation of Adult ERP Data*

In what follows, the descriptions of the components of the analysis pipeline for the third part of the study are organized according to the sequence of steps illustrated in Figure 1. From the information derived from the ICA (Step 1) and the estimated dipoles (Step 3), we derived simple single spiking representations adopting electric-fields estimation modeling (Step 5). Next, the results of Step 5 were used to implement an ERP-based Adaptive Control of Thought—Rational (ACT–R) modeling approach previously validated on the same task [43,44]. We organized the simulated dipoles and the simulated ERP spikes in patterns or chunks of activity corresponding to cortical areas postulated in the adult ACT–R (Step 6). The results of Step 6 allowed us to simulate the sequence of spatiotemporal activity as dipoles mapped onto the same adult structural MRI template as the children's data (Step 7). This was the basis for comparing children's and simulated adult estimated localization, therefore, testing for structural homology.

The results of Steps 6 and 7 were also used for aggregate series of simulated spikes to build polyspiking patterns for the simulated dipoles (Step 8). Finally, in Step 9, the simulated data obtained from Step 8 were converted to ERP topographic maps. This final step made possible to contrast children's and simulated ERP topographic mappings, therefore, testing for functional homology.

#### 2.3.1. Electric-Field Spiking Modeling (Step 5)

To simulate the electrical activity, each module in the neurocognitive model was assumed to be generating one or two dipoles in the location identified in the dipole-fitting stage. The module was assumed to produce its electrical energy in a rising and falling spike. For modeling purposes, a simple triangular wave was assumed, which peaked at the center of the module. The resulting electric field (voltage) was then calculated at the surface of the head for each electrode as the sum of the individual dipole contributions. Elsewhere, we have shown that this method generates reliable, valid and consistent descriptions of actual ERP activity [43,44]. Since the spiking activities of the components occurred at di fferent times in the observed data, it was not necessary to add the e ffects of more than one dipole at a given time.

The e ffect of each dipole was estimated in the simulation by following three steps (see Figure 3): (i) The square of the distance r from the dipole to an electrode was calculated by using Pythagoras. Next, (ii) the cosine of the angle θ between the electrode and the dipole was calculated by using vector dot product. Successively, (iii) the electric potential from the dipole at the electrode was derived from Coulomb's law (=k.p.cos(θ)/*r*2), where p is the strength of the dipole and k is a constant. (It was not necessary to know the value of the constant since relative magnitudes were used in the model).

**Figure 3.** Schematic representation of the calculation of an electric dipole field spiking by simulation; the left panel shows the output of the calculated electric-field potential represented as a simplified spike at fixed timing (determined by the simulated production schedule, here shown at an arbitrary time point for example sake; for the actual simulated timings in the present study see the production schedule shown in Table 1).


**Table 1.** Simulated locations and predicted timing of spiking occurrence according to the present adapted Python ACT–R architecture.

Note: The timings shown in the table include the base constant and the increments from the 50-millisecond cognitive cycle (described in Section 2.3.2) and each of them represent the approximate estimated moment in which a given production process in the ACT–R model is completed.

#### 2.3.2. ACT–R Simulation (Step 6)

To model and simulate our data, we used an adapted version of John R. Anderson's Adaptive (Control of Thought—Rational) ACT–R [64]. In the general architecture, cognition is considered to arise from the parallel interaction of several independent modules. However, top-down processes are directed by the Procedural Module, which is meant to model procedural memory. ACT–R models procedural memory as a production system. Specifically, procedural memory contains production rules (i.e., if/then rules). Communication to and from the Procedural Module is managed by buffers and chunks (see Figure 4). Chunks in ACT–R are lists of predicated information (for example, "duck" could be represented by the chunk: Is-an:animal, Name:duckling, Color:yellow, Size:small).

**Figure 4.** The organization of information in ACT–R (Adaptive Control of Thought–Rational). Adapted with permission from [42].

Each buffer can contain one chunk at a time. Each module has at least one buffer, so there is a visual buffer, a declarative memory buffer, and so on. Modules receive instructions from their buffers and place the results of their activity in their buffers. Collectively, the buffers can be thought of as working memory; they can also be thought of as representing the current context of the task. Productions "fire" when their "*if* condition" matches the contents of the buffers. The "*then*" part of a production then alters the content of the buffers. Productions can only fire one at a time.

In our version of ACT–R, the productions represent electric-field potentials. We programmed the simulation to fit the midpoint of the vincentized bins used to parametrize the ERP data series, so that each production was conditioned to occur at successive steps of approximately ti + 50 ms, with i = 0, 100 ... , 600. The 50-millisecond cognitive cycle is assumed in many realistic modeling architectures besides ACT–R (for example, Soar, EPIC, GOMS, see [65] for discussion). The neurobiological plausibility of this 50-millisecond cycle has been demonstrated by spiking neural network models simulating well known time constants for the GABA-A receptors in the Basal Ganglia, which in ACT–R (and in various other architectures) is assumed to be responsible for the working memory central executive [66]. Therefore, in the practical implementation of the model we assumed a base production-completion time constant which corresponded to the size of the vincentized bins (100 ms) plus the 50-millisecond cycle constant. Each module contained functions (specifically, ex-gaussian convolutions) to determine activation levels during rest, during the task and the decay rate after firing. E ffectively, these functions determined the spiking behavior associated with a given production: Faster productions corresponded to more intense and more quickly decaying spiking.

To implement the simulation, we adopted Python ACT–R [67]. In particular, the present version of the model assumed that the caudate in the basal ganglia acts as the central coordinator (executive) of productions. The hippocampus controls declarative memory while the cingulate cortex controls attention to conflicting stimuli. Frontal cortex supports declarative memory while visual processing takes place in the occipital with further processing in the parietal (see Figure 3). With the time constraints as shown in Table 1, we implemented an ACT–R model which predicted that initially the visual module (occipital) would be activated by the displayed pictures of the target (duck) and turtle (distractor) and would place a representation of the picture in the visual bu ffer (parietal). Next, the "parietal" representation would be used to retrieve the label of the object and the appropriate instruction about what to do in response to the image of that animal from declarative memory (temporal), which in turn would be placed in the planning bu ffer (frontal) and initiate the motor program for the manual response, or stop it (basal ganglia).

#### 2.3.3. Adult ACT–R Simulated Dipole Mapping (Step 7)

As in the case of the actual children's dipoles, the ACT–R simulated dipole data was mapped on a standardized MRI template provided by the Yale Bioimage Suite Web [63].

#### 2.3.4. Adult ACT–R Simulated Spike Series (Polyspiking) (Step 8)

In this step, the simulated electric-field spikes were chunked in ordered series corresponding respectively to the standard localizations and predicted timing of spiking assumed in the programmed schedule for the firing of ACT–R productions (see Table 1). The output of this step was a set of polyspiking patterns including a minimum of six spikes for each of the twelve electrodes. Within each pattern, the spike with maximal activation derived from the dipole of interest at a given bin interval, while all other spikes reflected resting or decaying background activations from the other dipoles. The latter configuration permitted to obtain *coarse coding* which could quantify the amplitude of each spike in the aggregated pattern as a color category (see example in Step 8, Figure 1) corresponding to a standard RGB value so that it could be ordered along an intensity scale.

#### 2.3.5. Adult ACT–R Simulated ERP Topographical Mapping (Step 9)

The simulated topographical maps were obtained through the same utility of EEGLab as the one used for the actual ERP data. To feed the simulated data in a compatible format, the values of the color intensities corresponding to the polyspiking patterns were previously converted (in Stage 8) into an arbitrary relative voltage scale (with range from blue/black, −6 μV, to red/orange, 6 μV).

#### 2.3.6. Actual vs. Simulated Data Comparisons: Homology Tests

#### MRI-mappings and Structural Test

To run the structural comparison, we first estimated the margin of localization error by computing the range di fferences based on the matches resulting from the search nearest grey area procedure; the measures were in millimeters. Successively, we computed *z*-scores of these ranges (henceforth called *z-Ranges*) which permitted to compare the extent of variations in localization in the actual children's data against those in the adult simulated estimations.

#### Topographical Mappings and Functional Test

ERP activity comparison between the actual children's topographic maps and the topographic maps derived from the simulated ERP activity were computed as linear regressions of vincentized averaged ERP amplitudes in the topographic map of children against the corresponding ACT–R simulated data; the fit was assessed by calculating the coe fficient of determination, *R*2.
