Delta-Band Neural Responses to Individual Words Are Modulated by Sentence Processing

Sophie Slaats; Hugo Weissbart; Jan-Mathijs Schoffelen; Antje S. Meyer; Andrea E. Martin

doi:10.1523/JNEUROSCI.0964-22.2023

Abstract

To understand language, we need to recognize words and combine them into phrases and sentences. During this process, responses to the words themselves are changed. In a step toward understanding how the brain builds sentence structure, the present study concerns the neural readout of this adaptation. We ask whether low-frequency neural readouts associated with words change as a function of being in a sentence. To this end, we analyzed an MEG dataset by Schoffelen et al. (2019) of 102 human participants (51 women) listening to sentences and word lists, the latter lacking any syntactic structure and combinatorial meaning. Using temporal response functions and a cumulative model-fitting approach, we disentangled delta- and theta-band responses to lexical information (word frequency), from responses to sensory and distributional variables. The results suggest that delta-band responses to words are affected by sentence context in time and space, over and above entropy and surprisal. In both conditions, the word frequency response spanned left temporal and posterior frontal areas; however, the response appeared later in word lists than in sentences. In addition, sentence context determined whether inferior frontal areas were responsive to lexical information. In the theta band, the amplitude was larger in the word list condition ∼100 milliseconds in right frontal areas. We conclude that low-frequency responses to words are changed by sentential context. The results of this study show how the neural representation of words is affected by structural context and as such provide insight into how the brain instantiates compositionality in language.

SIGNIFICANCE STATEMENT Human language is unprecedented in its combinatorial capacity: we are capable of producing and understanding sentences we have never heard before. Although the mechanisms underlying this capacity have been described in formal linguistics and cognitive science, how they are implemented in the brain remains to a large extent unknown. A large body of earlier work from the cognitive neuroscientific literature implies a role for delta-band neural activity in the representation of linguistic structure and meaning. In this work, we combine these insights and techniques with findings from psycholinguistics to show that meaning is more than the sum of its parts; the delta-band MEG signal differentially reflects lexical information inside and outside sentence structures.

Introduction

During language comprehension, listeners recognize words, retrieve stored information about them, and use this knowledge to combine the words into phrases and sentences. Psycholinguistic experiments have long shown that the behavioral responses to words change under the influence of the syntactic and sentential context that the words appear in (Marslen-Wilson and Welsh, 1978; Tyler and Wessels, 1983; Katz et al., 1987). In a step toward understanding how the brain builds sentence structure, the present study concerns the neural readout of this process. We ask (1) whether low-frequency neural readouts associated with words systematically change as a function of being or not being in a sentence context and (2) whether neural readouts are modulated by purely lexical properties over and above sensory and distributional variables. We do this by contrasting MEG responses to words in sentences with word lists, the latter lacking any syntactic structure or coherent lexical and combinatorial meaning.

In psycholinguistic models, language comprehension is instantiated as a cascaded process in which information can flow bidirectionally (Marslen-Wilson and Welsh, 1978; Martin, 2016, 2020). Put simply, this means that speech sounds cue stored representations of words, and while the next words are being recognized, the retrieved information about words cues representations of phrase and sentence structure. At the same time, the already formed representations of sentences, phrases, and words cue lower-level representations: the information flows in two directions (Schoffelen et al., 2017).

As words are being combined into phrases and sentences, then, responses to words change as a consequence of the top-down information flow. Indeed, a long tradition of research in psycholinguistics has shown that words in sentences are recognized faster than those same words appearing in isolation (Marslen-Wilson and Welsh, 1978; Tyler and Wessels, 1983). This effect is so powerful that it reduces effects of properties of the words themselves, such as word frequency. In isolation, highly frequent words are recognized faster than low-frequency words. In sentence context, this effect tends to be reduced: low-frequency words are recognized faster in sentence context than in isolation, although there is little change in recognition times for the high-frequency words (Schuberth and Eimas, 1977; Simpson et al., 1989).

To gain a full understanding of human sentence comprehension, those in the field currently face the challenge of integrating these findings with knowledge of neural processing. Although previous studies provide insight into the neural correlates of sentence structure (Ding et al., 2016; Meyer et al., 2017; Nelson et al., 2017; Ding et al., 2018; Brennan and Martin, 2020; Kaufeld et al., 2020; Bai et al., 2022; Coopmans et al., 2022; Tavano et al., 2022; ten Oever et al., 2022), much about the process of building these structures remains unknown (ten Oever et al., 2022). Furthermore, although we know that the neural signal is sensitive to lexical information (Brodbeck et al., 2018a,b; Armeni et al., 2019; Weissbart et al., 2020; Heilbron et al., 2021) we do not know how neural responses to words are transformed in the process of comprehension.

In this study, therefore, we aim to add to our understanding of how the brain leverages linguistic information when building sentence structure by finding a neural readout of the context effect on responses to words above and beyond statistical predictability effects as quantified through entropy and surprisal. To this end, we analyzed a published MEG dataset by Schoffelen et al. (2019) of participants listening to sentences and word lists. Despite these conditions being the main experimental manipulation in this open dataset, they have not previously been directly compared. Using temporal response functions (TRFs), we disentangled delta- and theta-band responses to individual words from responses to the speech envelope and word onsets, as well as entropy and surprisal. This method allowed us to model any differences between the conditions that go beyond our difference of interest (structured/unstructured), and, as such, control for them. We compared the responses to individual words between word lists and sentences. Any differences between the lexical responses in these conditions reflect the effect of structure building on the processing of words.

The lexical response was modeled using word frequency. We chose this feature because word frequency is a proxy for the likely familiarity of the listener with the word and relatedly of ease of processing. Any modulation as a consequence of word frequency, therefore, captures the presence of word identity information in the signal. Furthermore, word frequency is unigram; in other words, it does not depend on the context. Therefore, the value corresponding to a given word is the same in a sentence and a word list. Differences between the neural readout of both conditions will therefore be because of the sentence context supplying structure and meaning and not the predictor itself.

We hypothesized that the delta-band responses to word frequency would be different in word lists and sentences as a consequence of the (in)availability of sentence context (Huizeling et al., 2022; Meyer, 2018; Meyer et al., 2020a,b). Studies that investigated the presence of lower-level features in the neural signal as a function of the availability of linguistic information suggest that lower-level features are represented by the delta-band neural signal more reliably when higher-level information is available. For example, mutual information between the speech signal and the neural signal is higher in the presence of structure and meaning (Kaufeld et al., 2020; Coopmans et al., 2022; ten Oever et al., 2022), and the strength of speech tracking is dependent on the listener's knowledge of the language (Molinaro and Lizarazu, 2018; Blanco-Elorrieta et al., 2020) and general comprehension (Keitel et al., 2018). Following these results, we expected a stronger presence of the word frequency response (the lower-level feature) in the sentence condition than in the word list condition (the higher-level information) in the delta band specifically. Theta-band effects tend to be found as a function of acoustic rather than abstract linguistic manipulations (Sohoglu et al., 2012; Molinaro and Lizarazu, 2018; Etard and Reichenbach, 2019; Blanco-Elorrieta et al., 2020). In this study, we expected to observe this distinction between delta- and theta-band activity through an absence of effects in the theta band.

Materials and Methods

To answer our research question, we analyzed a part of the open-access large multimodal MEG dataset (N = 204) Mother of all Unification Studies published by Schoffelen et al. (2019). In addition, we performed two types of control analyses, an analysis of a dataset published by ten Oever et al. (2022) and a set of simulations. Methods for all analyses are described below.

Participants

A total of 102 native speakers of Dutch (51 men, 51 women) with a mean age of 22 years (range, 18–33) were included in this analysis. In this half of the dataset, participants were presented with the stimuli auditorily (as opposed to the other half, where stimuli were presented visually). All participants were right-handed, reported normal hearing, had normal or corrected-to-normal vision, and had no history of neurologic, developmental, or linguistic deficits. All participants provided informed consent, and the study was approved by the local ethics committee (Committee on Research Involving Human Subjects in the Arnhem-Nijmegen region, The Netherlands) and followed guidelines of the Helsinki Declaration. Participants took part in an fMRI and an MEG session, during which they listened to sentences and word lists. Only the MEG data are included in the present study.

Materials

The complete set of stimuli consisted of 360 natural Dutch sentences of 9–15 words (mean, 11.6), with varying syntactic structures, and 360 word lists. To create the word lists, the words from the sentences were scrambled so that more than two consecutive words did not form a coherent fragment. The stimuli were recorded by a female native speaker of Dutch. The sentences were pronounced naturally. The word lists were pronounced with neutral prosody and with a clear pause between each word. The files were recorded in stereo at 44 100 Hz. The sentences had an average duration of 4.27 s (SD 0.61), and the word lists of 7.67 s (SD 1.04). During the postprocessing, the audio files were low-pass filtered at 8500 Hz and normalized so that all the audio files had the same peak amplitude and peak intensity. In the word list condition, the individual words were spliced together with variable silence between them. This created conditions with different acoustic properties. We address this issue in the sections below, beginning with MEG preprocessing. In both conditions, the transition from silence to speech was ramped at the onset and offset with a rise/fall time of 10 ms. Word onsets and offsets were determined manually for each audio file using Praat software (Boersma and Weenink, 2018).

The stimuli were divided over two sets, A and B. During the MEG session, participants were presented with 120 sentences from set A and 120 word lists from set B (or the reverse). Across participants, all stimuli were presented the same number of times in the sentence and word list condition.

Procedure

Before the task, participants read written instructions and were allowed to ask clarification questions. The experimenter emphasized that the sentences and word lists should be attended to carefully and discouraged attempts to integrate the words in the word list condition. To familiarize the participants with the task, all participants performed a practice block with stimuli not included in the study. During the MEG measurement, the stimuli were presented in 24 blocks, alternating between sentence blocks (each containing five sentences) and word list blocks (each containing five word lists). The starting block type (either sentences or word lists) was randomized across participants. At the start of each block there was a 1500 ms presentation of the block type: zinnen (sentences in Dutch) or woorden (words in Dutch). The intertrial interval was jittered between 3200 and 4200 ms. During this period, an empty screen was presented, followed by a fixation cross.

To ensure participants paid attention to the stimuli, 20% of the trials were followed by a Yes/No question about the content of the preceding sentence/word list. Half the questions on the sentences addressed the content of the sentence (e.g., Did grandma give a cookie to the girl?), whereas the other half and all the questions about the word lists addressed one of the main content words (e.g., Was a grandma mentioned?). Participants answered the question by pressing a button for Yes/No with their left index and middle finger, respectively. Although the tasks were not identical between the conditions, the randomized order of appearance of question types ensured that participants could not approach the sentences any differently from the word lists; any sentence or list trial could be followed by the word monitoring task.

The stimuli were presented via plastic tubes and ear pieces in both ears. The hearing threshold was determined individually for each participant before the experiment, and the stimuli were presented at an intensity of 50 dB above the hearing threshold.

The experiment was run using Presentation software (version 16.0, Neurobehavioral Systems, www.neurobs.com). MEG was continuously recorded with a 275-channel axial gradiometer system (CTF) at a sampling frequency of 1200 Hz (cutoff frequency of the analog antialiasing low-pass filter was 300 Hz). Three head localizer coils were attached to the participant's head (nasion, left- and right ear canals) to determine the position of the head relative to the MEG sensors. The head position was monitored throughout the measurement. If needed, the participant was asked to reposition to correct for head position changes during breaks. The audio signal of the stimuli presented in the scanner was recorded along with the MEG data using an analog-to-digital converter channel.

Structural MRI images for source reconstruction were acquired using a T1-weighted magnetization-prepared rapid gradient echo pulse sequence with the following acquisition parameters: volume TR = 2300 ms, TE = 3.03 ms, flip angle = 8 degrees, 1 slab, slice matrix size = 256 × 256, slice thickness = 1 mm, field of view = 256 mm, isotropic voxel size = 1.0 × 1.0 × 1.0 mm. A vitamin E capsule was placed as a fiducial behind the right ear to allow visual confirmation of left–right consistency.

MEG preprocessing

The MEG data were preprocessed with custom-written MATLAB scripts using the FieldTrip toolbox (Oostenveld et al., 2011; Donders Institute for Brain, Cognition and Behavior, Radboud University, The Netherlands; http://fieldtriptoolbox.org). Before filtering, the data were epoched from audio onset to audio offset. The epochs were baseline corrected and bandpass filtered into the designated frequency band using a windowed-sinc finite impulse response (FIR) filter (15 s data padded), after which they were resampled to 120 Hz for TRF estimation.

The frequency band of interest was defined on the basis of the rate of occurrence of words in the stimuli, the differences in speech–brain coherence between conditions, and the literature (Blanco-Elorrieta et al., 2020; Donhauser and Baillet, 2020; Molinaro and Lizarazu, 2018; Weissbart et al., 2020). The word rate in the word lists was 1.5 Hz (SD 0.1), and in the sentences 2.7 Hz (SD 0.3). To compute speech–brain coherence, we first computed the broadband speech envelope by taking the absolute value of the Hilbert transform of the speech signal, low passing it at 20 Hz, and scaling the output between zero and one. We computed the magnitude squared coherence estimate of the broadband speech envelope and the MEG signal using Welch's method. The differences between word lists and sentences were estimated using a cluster-based permutation test. This revealed three peaks in the low-frequency signal—one between 1 and 3 Hz, one between 4.5 and 7 Hz, and one between 9.5 and 12 Hz (Fig. 1; Lam et al., 2018). On the basis of these clusters and frequency bands analyzed in the literature (Donhauser and Baillet, 2020), we analyzed two frequency windows, delta (0.5–4 Hz) and theta (4–10 Hz). To account for differences in speech–brain coherence that were exclusively because of acoustic differences between the conditions, we included the speech envelope as a predictor in all the models of the data (Fig. 1B, modulation spectra). Details of the models are below in Temporal response functions and Stimulus representation.

Figure 1.

A, Speech–brain coherence. Shaded area indicates SD. Black bars indicate frequencies that were part of clusters that contributed to the significant difference between sentence and word list coherence. B, Modulation spectra of the broadband speech envelopes (part of the TRF base model). The modulation spectra were obtained by concatenating the stimuli per stimulus type and performing a fast Fourier transform on snippets of 5 s. The resulting spectra were averaged.

Source reconstruction

MRI images were coregistered to the MEG headspace coordinate system by aligning the positions of the preauricular points and the nasion MEG coil to the MRI images using the MNE-Python coregistration GUI (Gramfort et al., 2013). For each participant, we reconstructed the cortical surface using the watershed algorithm from FreeSurfer. We created a surface-based source space with oct6 spacing, meaning ∼5 mm was between the source points. This generates 4098 sources per hemisphere. We created a single-layer Boundary element model (BEM) model with surface ico downsampling of 5120, from which the lead field was computed. The sources were reconstructed using a scalar Linear-constraint minimum-variance (LCMV) beamformer approach with a unit-noise gain beamformer to deal with depth bias. The data covariance used for computing LCMV filters was whitened using the covariance matrix of resting-state data. The resting-state data were bandpass filtered into the appropriate frequency band (i.e., 0.5–4 Hz for the delta band, and 4–10 Hz for the theta band). After application of the LCMV beamformer filters to the epoched MEG data, the source-localized epochs were morphed to fsaverage for group statistics. These source-localized, morphed epochs were then entered into the pipeline for temporal response function estimation. Source localization failed for 11 participants because of convergence issues for the noise covariance matrix or missing resting-state data (N_source = 91).

Temporal response functions

To characterize the effect of linguistic structure and meaning on the neural response, we estimated TRFs with different acoustic and linguistic features. This approach has been used to determine responses to different linguistic features, ranging from the speech envelope and phonemic information (di Liberto et al., 2015; Donhauser and Baillet, 2020) to lexical information (Broderick et al., 2018; Weissbart et al., 2020) and even syntactic embedding (Nelson et al., 2017). The response function of interest here is the response to word frequency as this is a unigram feature and therefore has the same per-word values in both conditions.

The TRFs were estimated using linear regression. We modeled the neural response by convolving the TRF kernel with the stimulus representation signal. In summary, this method reduces to a multivariate multiple linear regression, where we used lagged time series of stimulus features as predictors. The model equation reads as follows: yc(t)=∑f∑kxf(t)βf(t−τk) + η(t),(1) where {yc}t,{xf}t,{βf}t represent the recorded MEG signal of channel c, the input feature f, and its temporal response function respectively; {η}t is a Gaussian noise process accounting for measurement noise. This linear model can be easily rewritten in its vectorized form and further concatenated such that we model at once all channel equations independently. We estimate the coefficients of the TRFs βf̂ by minimizing the squared error between the measured MEG signals and the reconstructed signal obtained from Equation 1 while keeping the norm of TRFs coefficients, ||β||2 low to avoid overfitting. This minimization problem is solved in a closed form by the following: β̂=(XTX +λId)−1XTY,(2) where Y∈ℝN×C is the matrix representation of the measured MEG signal (for C channels arranged columnwise, each with N data samples); β̂∈ℝ(K.F)×C contains the estimated TRFs with K lags, F features for all C channels; X∈ℝN×(K.F) is a matrix containing all lagged feature time series of length N; λ is a regularization coefficient; and Id is the identity matrix. The regularization coefficient is needed to avoid overfitting, which in this case translates to the square matrix XTX not being full rank. Numerically, small eigenvalues or simply ill-conditioned matrices suffice to make the inversion unstable and thus will require regularization. In our case, this happens when features present some amount of autocorrelation (as columns of X are a time-lagged version of other columns). Continuous regressors such as the acoustic envelope (see below, Stimulus representation) will present strong autocorrelation and thus call for regularization.

In Equation 1, the vector of weights βf(t) represents the coefficients parametrizing the temporal response functions. They form a time course reminiscent of an event-related potential that tells us at which point in time (and, potentially, where) a feature modulates the neural signal. Thus, an increase at a certain lag for a given feature reflects an increase in the associated brain response to this feature at that given sensor and at the given time lag after stimulus onset. The concept of stimulus onset, especially for a continuous regressor such as the envelope, here reduces to a situation where the brain would be stimulated by an impulse of sound. Eventually, we estimate, from a system identification perspective, the transfer function mapping input to output when the brain is considered as a linear time-invariant system.

To evaluate how our models perform at reconstructing the neural data, we computed the Pearson's correlation coefficient between the true data and data reconstructed using the estimated TRFs. The correlation between the reconstruction and the original MEG indicates how much of the variance in the neural signal is explained by the features. The TRFs were not estimated on the same portion of data used to score the model. As further explained (see below, Model fitting), we used a nested cross-validation procedure to tune the regularization parameter, estimate the TRF coefficients, and finally score the resulting model. Unless specified otherwise, all analyses described below were done with custom-made Python scripts using MNE-Python (Gramfort et al., 2013). The whole analysis was conducted both in sensor and in source space.

Stimulus representation

Its multivariate character makes the TRF especially suitable for the current analysis: it allows for controlling for differences between conditions that are not currently under discussion by modeling them. To characterize the speech signal and part of its linguistic content, we constructed the following five different features: word frequency (the feature of interest) and four control features consisting of the speech envelope, word onsets, entropy, and surprisal.

The speech envelope feature was computed for each stimulus by taking the absolute value of the Hilbert transform and downsampling it to 120 Hz to match the downsampled MEG sampling rate. The envelope feature was added to represent the acoustic response and as such captures the difference between conditions observed in the cerebro-acoustic coherence that was caused by differences in the acoustic input (Fig. 1A,B).

The word onset feature was added to capture broadly any time-locked response to word onset for which the variance is not already explained by other features. As such, this feature can also capture any effects of segmentation that were different between the conditions. The word onsets and offsets were transcribed manually for each stimulus. We used a train of unit impulses, where the feature signal is one at the word onset sample and zero otherwise as follows: x(t)=∑wordsδ(t−tonset).(3)

These impulse trains were convolved with a Gaussian kernel with an SD of 15 ms. Such temporal smoothing has the effect of inflating the autocorrelation of the signal. We designed the width of this smoothing so that the smoothed impulses end up with energy spanning a frequency band comparable to our continuous regressor (envelope). The Fourier transform of a Gaussian is also a Gaussian, and the 15 ms SD of the temporal smoothing kernel equates to a spectral SD of 21.22 Hz. This ensured that all features required a similar degree of regularization in the regression analysis and made it possible to include impulse-like features such as word onsets and the envelope in the same regularized regression. Notably, this also translates into some uncertainty about or knowledge of the exact word onset timings.

Like the word onset feature, the word frequency feature was constructed as an impulse train of zeros everywhere but at word onset. Here, we used the respective word frequency value to modulate the height of the impulses. We used the log-transformed value of occurrence per million words, obtained from the SUBTLEX-NL corpus (Keuleers et al., 2010), as follows: xwf(t)=∑words−log(P(w))×δ(t−tonset),(4) where P(w) represents the unigram probability estimated from occurrence per million words.

If a word did not exist in the corpus, the fallback value of 0.301 (log/million) was used, corresponding to the lowest word frequency in the corpus. The values were z-scored across all stimuli. The resulting signal was convolved with the same Gaussian kernel as the word onset feature.

The entropy feature consists of lexical entropy, a weighted probability measure that quantifies the uncertainty about the upcoming word on the basis of the previous words. It provides a numeric answer to the question, Given the n previous words, with what degree of certainty can we predict the upcoming word? as follows: H(wi)=−∑kP(wk|wi−1...wi−n)×log2(P(wk|wi−1...wi−n)).(5)

The value was derived from a trigram model trained on the NLCOW2012 corpus using WOPR (van den Bosch and Berck, 2009). If a value was missing, the average of all entropy values was used. Like the word frequency feature, the entropy values were z-scored relative to all stimuli and inserted in a stick function, after which the stick function was convolved with the same Gaussian window. This feature was added to ensure that any effects on the word frequency feature were of a compositional semantic and structural nature rather than a probabilistic one.

The surprisal feature reflects how surprising a given word is in its immediate context. From an information-theoretic perspective, this reflects the information content, or self-information, of a word. It was calculated as the log 10 transformation of the conditional probability of a word, which was taken from the same trigram model as the entropy values. This means that surprisal is always based on the two preceding words; given the two preceding words, how high was the chance that the observed word would indeed appear? If the chance was low, surprisal is high. The feature was constructed in the same way as the word frequency and entropy features; the values were z-scored across all stimuli, inserted in a stick function at word onsets, and convolved with the Gaussian window as follows: I(w)=−log10(P(wi|wi−1...wi−n)).(6)

Because the three numerical lexical features (frequency, entropy, surprisal) might be correlated to some extent, we need to assert that the degree of multicollinearity present in our stimulus representation will not hinder the TRF coefficient interpretation. We checked whether the variance inflation factor (VIF) was below five (considered a relatively conservative measure of multicollinearity; Sheather, 2009; Tomaschek et al., 2018). The VIF was computed by correlating the z-scored entropy, surprisal, and word frequency values and by taking the diagonal of the inverted correlation matrix. This was done for all the stimuli and for both conditions separately. The VIF was never higher than five; the highest VIF was for surprisal at 4.8 in the word list condition.

Model fitting

The features were fitted in a cumulative manner to assess the contribution of each feature. This led to a total of seven models per frequency band, an Envelope model, consisting of only the speech envelope; an Onset model, consisting of the speech envelope and the word onset features; a Frequency model, consisting of the speech envelope, word onset, and word frequency features; an Entropy model, containing the speech envelope, word onset, and entropy features; a Surprisal model, consisting of the speech envelope, word onset, and surprisal features; and cross-combinations of those with and without the word frequency feature. An overview of all models and the corresponding features is provided in Table 1.

View this table:

Table 1.

The fitted encoding models

Before model fitting, the data were split pseudorandomly into a training and testing set at a 80/20 ratio. Care was taken that the sentences and word lists were evenly divided across the training and test sets. The sentence and word list models were each trained on 96 of 120 trials. The regularization parameter was optimized individually per participant, frequency band, and model (but not per condition) using an eightfold cross-validation procedure with 20 log-spaced values around the eigenvalues of the covariance matrix of the lagged speech envelope (λ = 60470.9) ranging from λ × 10⁻³ to λ × 10³. The best regularization parameter was determined as the value for which the average (across sensors) reconstruction accuracies were highest. Occasionally, reconstruction accuracies would not increase with a higher degree of regularization; instead, increasing the regularization would leave the reconstruction accuracy at the same value until overregularization occurred and reconstruction accuracy went down. In this case, the highest lambda value before a drop in accuracy occurred was chosen to ensure some degree of regularization. Each model was fitted on the complete training set using the regularization parameter from the cross-validation procedure, yielding the TRFs.

In the analysis of the source-localized MEG data, the manipulations were simplified because of computational limitations. The two maximal models were fitted, with word frequency as the only difference as follows: the Entropy/Surprisal model, consisting of the speech envelope, word onsets, entropy, and surprisal features; and the full model, consisting of all features. The cross-validation procedure was brought down to fivefold with 10 log-spaced values around the eigenvalue of the stimuli (60470.9) ranging from λ × 10⁻² to λ × 10².

Model evaluation

Each model was evaluated by convolving the estimated TRFs with the unseen stimuli from the test dataset. This yields, in essence, a prediction of the neural signal according to the model. The predicted neural signal was then correlated with the original neural signal from the test set using the Pearson product-moment correlation on a sensor-by-sensor or source-by-source basis. For every individual participant, this yielded a set of sensor- or source-based reconstruction accuracies for each model.

Statistical analysis

The TRF analysis has two deliverables. First, the TRF (the development of the estimated coefficients across time) is an ERP-like waveform that captures how the neural signal changes as a function of, for example, word frequency, and, second, the reconstruction accuracy, which is a metric of model fit. Here, we wanted to know (1) whether the responses to word frequency differ between sentences and word lists in time and space, so we compared the TRFs between conditions, and (2) whether the presence of the word frequency response differed between sentences and word lists, so we tested whether the word frequency predictor contributed differently to the reconstruction accuracy of a model in the two conditions.

Throughout, evaluation for statistical significance of the difference between TRFs was done using cluster-based permutation tests. Cluster-based permutation tests address the null hypothesis of exchangeability across conditions by a Monte Carlo estimate of the randomization distribution of a cluster-based test statistic, optimizing statistical sensitivity while controlling the false alarm rate. Here, we used the t statistic as the test statistic. In these tests, we create matrices of all sensors and samples. Then, we compute the difference between two conditions and express it as a t statistic for each of these data points. The t values are thresholded at an a priori threshold, and the thresholded t values are summed across clusters on the basis of spatial and temporal adjacency. The significance of the test statistic of the resulting largest cluster is compared with 1024 of similarly obtained test statistics, after random permutation of the condition labels. We used the function spatio_temporal_cluster_test from the MNE-Python library (Gramfort et al., 2013) with the t statistic as the test statistic and 1024 permutations.

To assess whether the responses to word frequency differed qualitatively between conditions in sensor space, the difference between the word frequency TRFs for the sentence and word list conditions was evaluated using a cluster-based permutation test. In addition, to characterize the response in each condition separately, we performed two cluster-based permutation tests with the same methods in which we contrasted the response against zero in each condition separately. In total, we performed three cluster-based permutation tests on the sensor TRFs, one on the difference between conditions and one on the TRF for each of the two conditions separately (against zero). In all cases, we calculated the threshold on the basis of the t distribution with a significance level of 5 × 10⁻⁸ with 101 (number of participants minus one) degrees of freedom. This equals three times the recommended threshold for the number of participants. The threshold was increased to yield the most informative results (i.e., to ensure not every sensor and time lag would be significant). Subsequent comparisons were done with a threshold calculated using a Bonferroni adjusted significance level (i.e., divided by two) to correct for multiple comparisons; everything else was the same.

In addition, we wanted to evaluate whether there was a latency difference between the responses in the two conditions. To this end, we compared the responses from the sentences and word list conditions in a cross-correlation. The cross-correlation was done on the grand-average TRF waveforms of overlapping sensors between conditions from the clusters resulting from the one-sample tests. We sequentially cross-correlated each sensor and normalized the values by dividing them by the maximal value from the cross-correlation for that sensor. We then obtained the peaks for every sensor. This number corresponds to the lag at which the two signals had the highest correlation and shows how different the responses are in time. Subsequently, we shifted the sentence response in time by the number of samples of the peak. We then correlated the shifted sentence response and the original word list response. To check for significance, we performed the same procedure for randomly selected channels and repeated this process 10,000 times.

In source space, we compared the TRFs for word lists and sentences using a cluster-based permutation test in two time windows on the basis of the results from the analysis in sensor space, 200–400 and 500–700 ms post stimulus onset (PSO), respectively. We did this to get a more reliable estimate of the spatial distribution of the effects, although cluster-based permutation tests account only for a difference between the distribution overall, therefore any spatial or temporal differences are approximations and inconclusive (Maris and Oostenveld, 2007; Sassenhagen and Draschkow, 2019). The threshold was set to the t distribution with an alpha of 0.025 (98.75th and 1.25th percentile) to correct for multiple comparisons, with 90 (number of participants minus one) degrees of freedom. Sources along the medial wall were excluded.

In the sensor space analysis, the reconstruction accuracies were averaged over sensors and submitted to a linear mixed model using lme4 in R software (Bates et al., 2015). The model had the factor condition (two levels, sentence and word list) and a random intercept for participant. In addition, the model contained three binomial factors, frequency, entropy, and surprisal, describing whether a feature was (1) or was not (0) in the model to calculate a slope for each feature separately as follows: accuracies ∼ condition *(frequency + entropy + surprisal)+(1|participant).

We used a stepwise variable selection to evaluate the contribution of each of these factors. To evaluate the contribution of a given factor (or interaction), a model with the factor was compared with a model without it, and the goodness-of-fit statistics were compared using a chi-square test. If the removal of a factor did not decrease goodness of fit, the next factor was removed. When the removal of a given feature or interaction significantly decreased model fit, the removal of features was stopped. The prefinal model should then describe the data best. As a final check, the Akaike information criterion (AIC) of the models was compared using the R package AICcmodavg (Mazerolle, 2020). Post hoc t tests were done between the Entropy/Surprisal and Full model to evaluate whether the effects held between the largest models.

In source space, a cluster-based permutation test was done to localize the interaction effect using the function permutation_cluster_test from the MNE-Python library. The test statistic was an F statistic from a two-way ANOVA with factors Condition (levels: word list, sentence) and Model (levels: Entropy/Surprisal, Full). The data were permuted 1024 times.

Control analysis I: data

The word lists were presented with variable silences between words. The sentences, on the other hand, were natural, with pauses occurring sparingly. This caused differences of word rate and signal length between the conditions that may affect our results. To examine potential effects of the pauses in the word list condition, we analyzed a second dataset of 16 participants listening to word lists and sentences using the same methods. Importantly, the word lists in this condition were naturally spoken, as were the sentences. This means that there were no pauses between the words in the word list condition, and there was coarticulation between words (Kaufeld et al., 2020). The data were supplied by ten Oever et al. (2022).

Control analysis I: Participants.

A total of 20 native speakers of Dutch (4 men, 16 women with a mean age of 39.5 years) participated in the experiment. Four participants were excluded from this analysis for a variety of reasons (e.g., session was not completed). All participants were right-handed, reported normal hearing, had normal or corrected-to-normal vision, and had no history of neurologic, developmental, or linguistic deficits. All participants provided informed consent. The study was approved by the ethical Commission for Human Research Arnhem/Nijmegen (project number CMO2014/288). Participants were remunerated for their participation.

Control analysis I: Materials.

The stimuli were identical to the stimuli used in Kaufeld et al. (2020). The experiment consisted of three conditions in total, sentences, jabberwocky, and word lists. Only the sentences and the word lists are analyzed here. The stimuli consisted of 10 words, which were all disyllabic except for de (the in Dutch) and en (and in Dutch). Sentences had a fixed syntactic structure of two coordinate clauses: [Adj N V N conj Det Adj N V N], for example, timid heroes pluck flowers and the brown birds gather branches. The word lists were scrambled versions of these sentences, and care was taken so there were no plausible internal combinations of words. The stimuli were recorded by a female native speaker of Dutch at a sampling rate of 44.1 kHz (monophonic). After recording, any pauses were normalized to ∼150 ms in all stimuli, and the intensity was scaled to 70 dB using Praat voice analysis software (Boersma and Weenink, 2018).

Participants were asked to perform four different tasks on these stimuli—a passive listening task, a syllable recognition task, a word recognition task, and a word combination recognition task. In this analysis, we did not distinguish among tasks. ten Oever et al. (2022) describes the tasks performed.

Control analysis I: Procedure.

At the beginning of each trial, participants were instructed to look at a fixation cross presented at the middle of the screen on a gray background. The audio was presented binaurally through tubes after an interval randomly jittered between 1.5 and 3 s. One second after audio offset, the task prompt (e.g., the syllables or words for recognition) was presented, which required participants to press a button on a button box. There were eight blocks of ∼8 min. After each block, participants could take a break, during which the head position was corrected. MEG was recorded using a 275-channel axial gradiometer CTF MEG system at a sampling rate of 1200 Hz. After the session, head shape was collected using the Polhemus digitizer (using as fiducials the nasion and the entrance of the ear canals as positioned with ear molds).

Control analysis I: MEG preprocessing

The MEG data were processed with custom-written Python scripts using MNE-Python (Gramfort et al., 2013). As in the main analysis, the raw MEG data were filtered using a windowed-sinc FIR filter between 0.5 and 4 Hz for the delta band, and 4 and 10 Hz for the theta band, after which the data were epoched from audio onset to audio offset and resampled to 120 Hz for TRF estimation.

Control analysis I: Stimulus representation

In this analysis, we used the envelope, word onset, and word frequency representations from the main analysis (see above, Stimulus representation).

Control analysis I: Model fitting

We used the model-fitting approach described earlier (see above, Model fitting). We fit three models, Envelope (with only the envelope feature), Onset (envelope and word onset features), and Frequency (envelope, word onset, and word frequency features). The data were split pseudorandomly into a training and a testing set at an 80:20 ratio, ensuring that the sets contained 50% of items from each condition. The regularization parameter was optimized individually per participant and model, using an eightfold nested cross-validation procedure with 20 log-spaced values around 60,000 (λ = 60,000) ranging from λ × 10⁻² to λ × 10².

Control analysis I: Model evaluation

For model evaluation, we used the procedure described earlier (see above, Model evaluation).

Control analysis I: Statistical analysis

Like in the main analysis, we assessed whether the responses to word frequency qualitatively differed between conditions by evaluating the difference between the word frequency TRFs for the sentence and word list conditions using a cluster-based permutation test. In addition, to characterize the response in each of the conditions separately, we performed two additional cluster-based permutation tests with the same methods in which we contrasted the response against zero in each condition separately. In total, we performed three cluster-based permutation tests on the TRFs, one on the difference between conditions and one on the TRF for each condition separately (against zero). In all tests, we calculated the threshold on the basis of the t distribution with a significance level of 0.05 with 16 (number of participants minus one) degrees of freedom. Only clusters with a p value smaller than 0.01 were considered. Subsequent comparisons were done with a threshold calculated using a Bonferroni-adjusted significance level to correct for multiple comparisons; everything else was the same. For comparison to the main analysis, we also compared the word onset response between conditions with the methods described above.

To evaluate the effect of word frequency in each condition, we compared the reconstruction accuracies from the Onset and Frequency models in interaction with condition. The reconstruction accuracies were averaged over all sensors (conservative measure). After checking for normality and sphericity through (1) visual inspection of Q-Q plots and histograms; (2) statistical testing using the Shapiro–Wilk test, Anderson–Darling test, and D'Agostino's K² test for kurtosis and skewness as implemented in SciPy algorithms; and (3) the Mauchly test for sphericity as implemented in the Pingouin package (Vallat, 2018), the averaged reconstruction accuracy values were submitted to a repeated-measures ANOVA using the Statsmodels package.

Control analysis II: simulations

Using simulations, we evaluated whether the interword interval has an impact on TRF model evaluation. We did this by simulating raw MEG data consisting of a signal (different impulse responses) and a variable amount of noise.

The simulated response was equivalent to the forward model, namely a noisy output of a convolution between a predefined kernel (the ground truth for the TRF estimate) and an impulse train (for the input signal). We generated those data with a variable amount of noise (i.e., explicitly manipulating the broadband signal-to-noise ratio) and with varying the interstimulus interval (ISI) while keeping the signal length the same and the number of impulses, or events, constant (in which case a shorter interstimulus interval results in the end portion of the output signal containing only noise).

We then scored the forward model by computing both the R² score and the Pearson's correlation coefficient between the reconstruction ŷ and the true signal using a test portion of the data, not used to estimate the coefficients β. Importantly, we then computed the scores in two ways, (1) from the fixed signal length data described above, as we also used a fixed number of impulses, or events, this resulted in a portion of the stimulated output signal to contain only noise (2) or from a shortened signal, where we truncated all signals to the last stimulus event. This resulted in shorter signals for a shorter ISI.

Data availability

The code is available at https://osf.io/ky9bj/, with the exception of the preprocessing scripts. The preprocessed data are available on request. The raw data can be downloaded from the Donders Institute repository at https://data.donders.ru.nl/collections/di/dccn/DSC_3011020.09_236?0.

Results

Behavioral results

We compared participants' responses to the task that was present in both conditions, which targeted one of the main content words (e.g., Was a grandma mentioned?). To balance the number of trials included in the accuracy scores, we took a random subset of questions from the word lists (12 or 13 trials). The average proportion of correct responses was higher in the sentence condition (mean_sent = 0.88; sd_sent = 0.08) than in the word lists (mean_list = 0.72; sd_list = 0.14; t = 10.08, p < 0.001), meaning that participants remembered the words from the sentences better than the words from the word lists (Fig. 2).

Figure 2.

Accuracy scores for the behavioral task performed during the MEG recording. The accuracy scores include responses to word monitoring only. The word list accuracy scores are a random subset of the full set of responses to balance the number of trials (n = 12) in the word list and sentence conditions.

Delta band

Sensor-level analysis

The cluster-based permutation test revealed differences between word lists and sentences in three clusters between 0 and 700 ms. Figure 3A suggests that the peak of the response to word frequency was delayed by ∼300 ms in the word list condition. To evaluate whether this was the case, we conducted one-sample cluster-based permutation tests and computed the cross-correlation between the two conditions for overlapping sensors from the clusters in both conditions. The one-sample cluster-based permutation test revealed a response in temporal areas in both conditions that peaks ∼250 ms in the sentence condition, and ∼600 ms in the word list condition (Fig. 3B,C).

Figure 3.

A, The word frequency TRF in both conditions in the delta band. Shown here is the mean of the sensors that were included in clusters that were different between the two conditions. Black bars indicate time points that contributed to clusters that allowed us to reject the null hypothesis. Shaded area indicates SD. B, Word frequency TRF in the sentence condition. Individual lines represent sensors. Sensors in bold contributed to the clusters that allowed us to reject the null hypothesis. C, Word frequency TRF in the list condition. Individual lines represent sensors. Sensors in bold contributed to the clusters that allowed us to reject the null hypothesis. D, The word onset TRF in both conditions in the delta band. Shown here is the mean of the sensors that were included in clusters that were different between the two conditions. Black bars indicate time points that contributed to clusters that allowed us to reject the null hypothesis. Shaded area indicates SD. Vertical gray lines indicate the time points of the scalp maps.

The cross-correlation on overlapping sensors between the two conditions (time courses and sensors; Fig. 4A) revealed a high correlation between the word list and the sentence responses at a delay of 330 ms (mean r = 0.9). Random sampling of sensors and lags revealed the distribution shown in Figure 4D; the observed values are in the upper 0.05% percentile, indicating that the observed correlation is likely not caused by chance.

Figure 4.

A, TRF time courses for shared sensors between the sentence (solid lines) and word list (dashed lines). Colors indicate sensor position. B, Cross-correlation between the sentence and word list responses for overlapping sensors between conditions from the clusters (scaled between −1 and 1). Colors indicate sensor position. C, The shifted response from the sentence condition (solid lines) to overlap with the word list condition (dashed lines). Colors indicate sensor position. D, Kernel density plots of means and SDs from correlations between randomly selected sensors at shifted randomly selected lags; the red bar indicates the values observed from the sensors selected after the cluster-based permutation test shifted at the lags from the cross-correlation. Coeff: coefficient.

Because we wondered whether the delay could be because of the differences in the presentation rate, we examined differences between the TRFs for the other word-level feature that was numerically identical between conditions, word onsets (unit-spike-train in both conditions). We compared the word onset response from a model with only the envelope and word onset features. This model is equivalent to an ERP analysis that corrects for overlapping event windows (as is the case in the sentence condition) and controls for acoustic differences. A small delay of ∼100 ms appears in this model. This delay is in accordance with findings of an ERP-analysis on high- versus low-constraining contexts (Liu et al., 2006; León-Cabrera et al., 2017). Importantly, this model collapses over variance caused by the lexical features included in the full model (word frequency, entropy, and surprisal). In other words, this underspecified model attributes variance that is in fact because of word frequency, entropy, or surprisal to the word onset predictor. When we include the other lexical predictors in the model and compared the conditions again, no such difference between the word onset responses is observed (Fig. 3D). In this response, there were some differences around time point zero before as well as slightly after; these differences may indicate differences in temporal expectancy of word onset between conditions.

The reconstruction accuracies were evaluated with the model accuracies ∼ condition * (frequency + entropy + surprisal) + (1/participant). The explanatory value of the interaction between condition and each of the lexical factors was evaluated; each interaction significantly improved model fit [frequency, χ²(1) = 6.88, p < 0.01; entropy, χ²(1) = 4.48, p < 0.05; surprisal, χ²(1) = 7.24, p < 0.01], so the full model was interpreted. The results of this model are summarized in Table 2.

View this table:

Table 2.

Results of the LME on the reconstruction accuracies in the delta band

View this table:

Table 3.

Results of the LME on the reconstruction accuracies in the theta band

Reconstruction accuracies were higher in the word list condition than in the sentence condition (β = 1.67 * 10⁻², SE = 9.43 * 10⁻⁴, t₍₁₅₃₀₎ = 17.69, p < 0.01). As can be seen in Figure 5A, each feature contributed positively to the reconstruction of the neural signal in the sentence condition, less so in the word list condition, hinting at an interaction effect. Indeed, the factor frequency interacted with condition (β = 2.47 * 10⁻³, SE = 9.43 * 10⁻⁴, t₍₁₅₃₀₎ = 2.63, p < 0.01), showing that reconstruction accuracies improved more from the addition of the word frequency predictor in the sentence condition than in the list condition (Fig. 5B). Further, although we do not discuss these effects, entropy and surprisal interacted with condition as well (entropy, β = 2.00 * 10^–3, SE = 9.43 * 10⁻⁴, t₍₁₅₃₀₎ = 2.12, p < 0.05; surprisal, β = 2.54 * 10⁻³, SE = 9.43 * 10⁻⁴, t₍₁₅₃₀₎ = 2.69, p < 0.01).

Figure 5.

Reconstruction accuracies in the delta band. A, Reconstruction accuracy difference with the envelope model for each model in the sentence condition. Middle line indicates the median, the white diamond indicates the mean. B, Reconstruction accuracy difference with the envelope model for each model in the word list condition. Middle line indicates the median, the white diamond indicates the mean. C, The interaction between condition and frequency on the reconstruction accuracies. Values on the y-axis are the difference with the envelope (as in A, B). Error bars represent the 95% confidence interval. Entr.: entropy, surp: surprisal, freq: frequency.

To gain more insight into the effect of frequency, we performed a post hoc t test comparing the two largest models (Entropy/Surprisal and Full). These tests confirmed that the word frequency predictor enhanced reconstruction accuracy in the sentence condition (t₍₁₀₁₎ = 5.35; p < 0.01), but not in the word list condition (t₍₁₀₁₎ = −0.15, p = 1; Bonferroni corrected).

Finally, we hypothesized that the higher reconstruction accuracy in the word list condition was because of the salience of isolated words, possibly evoking a larger auditory response. If this is true, a model with only the envelope predictor, and no word-level feature, should also fit the list condition better. To evaluate this hypothesis, we compared the reconstruction accuracies (averaged over all sensors) for the Envelope model between conditions. This model was not included in the analyses of the word frequency effect. And indeed, this was the case; reconstruction accuracies were higher for word lists than sentences using only the envelope as predictor (t₍₁₀₁₎ = 13.40, p < 0.01).

In sum, the response to word frequency differed between word lists and sentences. The TRFs in sensor space revealed a left-lateralized frontotemporal response to the feature that peaked ∼250 ms after word onset in the sentence condition, and ∼600 ms in the word list condition. The sentence effect is in line with other studies that used word frequency as a feature in TRF models of natural language comprehension (Brennan and Hale, 2019; Weissbart et al., 2020). A cross-correlation analysis between a set of left (and one right) temporal and frontal sensors that were involved in the response in both conditions suggested that the word list response peaks ∼300 ms later. The reconstruction accuracies in sensor space suggests that the word frequency predictor explains more variance over and above acoustics, entropy, and surprisal in the sentence condition, but not in the word list condition.