Charlotte Caucheteux

Alexandre Gramfort

Jean-Rémi King

Model-based analysis of brain activity reveals the hierarchy of language in 305 subjects

Model-based analysis of brain activity reveals the hierarchy of language in 305 subjects

Charlotte Caucheteux, Alexandre Gramfort, Jean-Rémi

Introduction

One of the most successful paradigms to decompose the brain bases of language consists in correlating the brain responses of multiple subjects listening to the same carefully controlled stimuli [START_REF] Brennan | Syntactic structure building in the anterior temporal lobe during natural story listening[END_REF][START_REF] Fedorenko | Neural correlate of the construction of sentence meaning[END_REF][START_REF] Blank | Syntactic processing is distributed across the language system[END_REF][START_REF] Mollica | Composition is the core driver of the language-selective network[END_REF]. In particular, [START_REF] Lerner | Topographic Mapping of a Hierarchy of Temporal Receptive Windows Using a Narrated Story[END_REF] recorded subjects with functional magnetic resonance imaging (fMRI) while they listened to a story whose (1) sounds (2) words, (3) sentences or (4) paragraphs were scrambled, as well as ( 5) to the regular version of the story (Figure 1A). The authors then estimated the Inter Subject Correlation (ISC), i.e., the correlation between i) the brain activity of a voxel in response to one scrambling condition and ii) the brain activity of a voxel averaged across all other subjects, in response to the same scrambled stimulus (Figure 1B). While successful, this 'model-free' approach is costly: it requires n subjects ×n conditions acquisitions of brain activity in response to the same variably scrambled stimuli.

Here, we investigate whether and how a model-based approach can replicate Lerner et al.'s findings, even if we only have access to the recordings elicited by the regular story in a single subject. We further apply the method to extend Lerner et al's results to a large dataset of 305 individuals.

Methods

First, we formalize the 'model-free ' and 'modelbased' approaches in the context of narrative listening, and explicit the link between the two.

Definitions Let's define

• w = ('Once', 'upon', ... , 'The', 'end.') the regular story. Ω the story's vocabulary.

• w |sound , w |word , w |sent , w |parag the story scrambled at the acoustic, word, sentence and paragraph level, respectively, following the setting of Lerner et al. (cf. Appendix ?? for the scrambling paradigm).

• : Ω M → R T : the function returning the brain recordings of length T time samples (i.e., the number of fMRI pulses) induced by a sequence of M words. This method allows to decompose the hierarchy of language processing in the brain, from the acoustic to the paragraph level.

C. We aim to replicate the results of Lerner et al. using only the recordings induced by the regular story (black). To this aim, we scramble, not the stimulus of the subject, but the inputs of a deep language model (GPT-2). For each condition (word, sentence or paragraph), we extract the corresponding activations x * averaged over K random scrambles. We then compare the brain signals of the current subject y with the activations x * elicited by the scrambled texts, after a linear transformation f θ that maps x * onto a brain-like space. Because GPT-2 is not trained to process waveform, we use the phonemes, stresses and tones of the stimulus instead of x * for the acoustic condition.

• : Ω M → R M ×D the function returning the activations of a deep language model induced by a sequence of M words.

• y ∈ R T the brain recordings of one subject elicited by w, recorded at one voxel. Here, (w) = y.

• y |sound , y |word , y |sent , y |parag the recordings elicited by the scrambled versions of w.

• ρ : R T × R T → R, Pearson's correlation

For clarity, we describe below the model-free and model-based approaches for the sentence condition. The same methods can be used for the sound, word and paragraph conditions. Lerner et al. do not have a model of how the brain should react to sentences. Instead, they assume that the neural signature of sentence-level processing corresponds to the brain response shared across all subjects listening to scrambled sentences w |sent . They thus compute the 'ISC score' for each subject, i.e., the correlation between i) the brain response to the scrambled story w |sent of a given subject (y |sent ) and ii) the brain response to the same stimulus averaged across all other subjects (y |sent ): R = ρ y |sent , y |sent .

Model-free analysis

(1)

This approach boils down to a leave-onesubject-out cross-validation, using Pearson correlation as evaluation metric and the average population response as estimator.

Model-based analysis

Here, we propose a model-based analysis to circumvent the need for a) the population average y, b) the scrambled stimuli w |sent .

To eliminate the need for the population average, we capitalize on the recent findings that deep language models tend to linearly predict brain responses to language [START_REF] Jain | Incorporating Context into Language Encoding Models for fMRI[END_REF][START_REF] Gauthier | Linking artificial and human neural representations of language[END_REF][START_REF] Toneva | Interpreting and improving natural-language processing (in machines) with natural language-processing (in the brain)[END_REF][START_REF] Schrimpf | Artificial Neural Networks Accurately Predict Language Processing in the Brain[END_REF][START_REF] Caucheteux | Language processing in brains and deep neural networks: computational convergence and its limits[END_REF]. We can thus assume that the average brain response ( ) can be well approximated by f θ , a linear function that maps the deep language model to the brain response. i.e.,

i)

≈ f θ • .

In practice, the coefficients θ of f θ are estimated using ridge regression. Finite Impulse Response functions are employed to allow It is considered to process 'word '-level (yellow) if its word score is significant but not its acoustic oneand similarly for 'sentence' (green) and 'paragraph' (blue). A. Adapted from [START_REF] Lerner | Topographic Mapping of a Hierarchy of Temporal Receptive Windows Using a Narrated Story[END_REF]. Labels are based on the brains-to-brain correlation scores (Figure 1B) averaged over seven subjects listening to a 7 min story. B. Labels are based on the model-to-brain scores (Figure 1C), averaged over 75 subjects listening to the same 7 min story. Significance is inferred using a Wilcoxon test across subjects, corrected with False Discovery Rate (FDR) across the 465 brain regions in each hemisphere (cf. Appendix C), with a significance threshold of p < 10 -3 (cf. Appendix D). C. Same as B., but on the brain of 305 subjects listening to 4 hours of 15 audio stories (including the 7 min one). Because of the large number of subjects, the significance threshold is set to p < 10 -25 . the activations of the deep language model of length M (number of words) to map onto the slow and delayed brain recordings of length T (number of pulses) (cf. Appendix B).

To eliminate the need for the scrambled stimuli, we show below that equation (1) can be rewritten only as a function of w as opposed to w |sent .

First, we separate the representation of the sentence from that of its context. To this end, for each sentence s of w, we note Ω s the set of sequences ending with s, and whose preceding context is random. The representation of s without context, is, by construction, also the sentence representations of all sequences w ∈ Ω s . Thus, if we denote y * s this common representation, the brain response of one subject to a sequence w can be modeled as

∀w ∈ Ω s , (w ) = y * s + ε w , (2)
with ε w the context-dependent contribution to (w ). Assuming it is a zero-mean random perturbation we have:

E w (w )] = y * s , (3)
with w sampled uniformly in Ω s . Importantly, we do not assume that words are independent of their context but that the shufflings defined for each sentence are independent of one another. This statement is true by construction:

shuffled contexts are realizations of a uniform sampling of permuted texts. Furthermore, the assumption that activations of shuffled versions of the same context have a zero-mean is not critical: assuming a constant mean would not alter the methods and results, because the final metrics (Pearson correlation) is invariant to such constant. Similarly, we can retrieve x * s , the contextindependent representation of a particular sequence s in a deep language model

E w (w )] = x * s . (4)
In practice, it is approximated with an average over K i.i.d. samples:

x * s ≈ 1 K K k=1 (w k ) , (5)
where w 1 , . . . , w K are sentences uniformly sampled in Ω s . Given equations (3), (4) and hypothesis i).

ȳ * s = E w (w ) = E w [f θ • (w )] ⇒ ȳ * s = f θ (x * s ) , (6)
with w sampled in Ω s .

From now on, we note y * (resp. x * ) the context-free representation of the whole story w extracted from the brain (resp. network) activations. We obtain, ȳ * = f θ (x * ).

We now assume that random contexts do not actually affect the brain response to the current sentence in each subject at a given voxel, i.e., ii) y |sent = y * |sent = y * .

Under this condition and given equation ( 6),

y |sent = y * = f θ (x * ) , and
ρ y |sent , y |sent ) = ρ y * , f θ (x * ) = ρ y -ε w , f θ (x * ) ≈ ρ y, f θ (x * ) ,
with ε w the strictly contextual effects in y (y = y * + ε w ), independent from x * (context-free).

Finally, under the assumptions that i) the deep neural network approximates the average brain response and ii) random context is not maintained in memory, the brains-to-brain scores R = ρ y |sent , y |sent are equivalent to the model-to-brain scores R = ρ y, f θ (x * ) .

Experiment

To test our model-based approach, we first apply it to the fMRI responses of 75 subjects listening to the same 7 min story analysed in Lerner et al (Nastase et al., 2020) 1 . Thus, for each condition (word, sentence and paragraph), subject and voxel, we compute the model-tobrain correlation R = ρ y, f θ (x * ) .

The extraction of the fMRI signals y, and the estimation of the mapping function f θ are standard and thus detailed in Appendices A and B. To estimate context-free representations, we i) scramble the stimulus at the word, sentence or paragraph level, ii) extract the corresponding activations x from a deep language model, and iii) compute x * , as detailed below.

Scrambling the stimulus at the word, sentence and paragraph level Words and sentences of the stimulus are delimited using Spacy tokenizer [START_REF] Honnibal | spaCy: Industrial-strength Natural Language Processing in Python[END_REF]. Note that punctuation marks are not considered as words (e.g., 'time.' forms one token, not two). We define paragraphs as contiguous chunks of eight sentences. To 'scramble' a sequence at the word (resp. sentence, paragraph) level, we uniformly shuffle the indices of its words (resp. sentences, paragraphs) and form the new sequence accordingly.

Extracting deep models' activations

For each version of the scrambled stimulus, we extract the activations from GPT-2 ( ), a deep neural language model trained to predict a word given its past context. GPT-2 consists of 12 transformer layers of dimensionality 768, 8 heads, and has 1.5 billion parameters in total. We use the model provided by Huggingface [START_REF] Wolf | Transformers: Stateof-the-art natural language processing[END_REF], trained on a dataset of 8 million web pages.

To extract the activations elicited by a sequence w of M words from layer l, we proceed as follows: we tokenize the sequence into subwords called "Byte Pair Encoding" (BPE) [START_REF] Sennrich | Neural Machine Translation of Rare Words with Subword Units[END_REF] using the GPT-2 tokenizer provided by Huggingface. Then, we feed the network with the M BPE tokens (M ≥ M , up to 256 tokens in memory) and extract the corresponding activations from layer l, of shape (M × D) with D = 758. Then, we sum the activations over the BPEs of each word to obtain a vector of size (M × D).

All our analyses are based on the eighth layer of GPT-2. We choose GPT-2 because it has been shown to best encode the brain activity elicited by language stimuli [START_REF] Caucheteux | Decomposing lexical and compositional syntax and semantics with deep language models[END_REF][START_REF] Schrimpf | Artificial Neural Networks Accurately Predict Language Processing in the Brain[END_REF]. We choose its eighth layer because the intermediate layers of transformers have shown to encode relevant linguistic features [START_REF] Jawahar | What Does BERT Learn about the Structure of Language[END_REF][START_REF] Christopher | Emergent linguistic structure in artificial neural networks trained by self-supervision[END_REF] and to better encode brain activity than input and output layers [START_REF] Caucheteux | Language processing in brains and deep neural networks: computational convergence and its limits[END_REF][START_REF] Toneva | Interpreting and improving natural-language processing (in machines) with natural language-processing (in the brain)[END_REF]. Our results successfully generalize to two other architectures as well as to the other intermediate layers of GPT-2 (Appendix E).

Computing x * for the word, sentence and paragraph conditions For each of the word, sentence and paragraph conditions, we compute x * : a context-free representation of x. In short, x * are the activations of GPT-2, averaged over several scrambled contexts. For clarity, we focus on the sentence level to detail the approach.

To build the sentence-level representation

x * of the stimulus, we use the approximation introduced in equation ( 5). For each sentence s of one story w, we i) generate K=10 sequences ending with s, but with scrambled previous context. The scrambled context is uniformly sampled from the other sentences in the same story w. Then, ii) we extract the K corresponding activations from GPT-2 (as described in the previous section) and iii) average the activations across the K samples. GPT-2 activations are extracted for each word. Thus, for each of the M s words of sentence s, we obtain a vector x * s of shape M s × D. We concatenate these vectors to obtain x * , a sentencelevel representation of the whole story w, of shape M × D. This method is adapted from [START_REF] Caucheteux | Decomposing lexical and compositional syntax and semantics with deep language models[END_REF], in which we computed the average over GPT-2's activations to extract syntactic representations from the input sequence.

Acoustic features GPT-2 takes words as input and not sounds. To build x * at the acoustic level, we simply use non-contextual acoustic features: the word rate (D = 1), phoneme rate (D = 1) phonemes, stress, and tone (categorical, D = 117). For the latter, we use the annotations provided the original Narratives dataset [START_REF] Samuel | Narratives: fMRI data for evaluating models of naturalistic language comprehension[END_REF].

Results

The results are displayed in Figure 2B. The hierarchy of temporal receptive fields (TRFs) typically associated with acoustic, word, sentence and paragraph processing along the temporoparietal axis is remarkably well replicated in both hemispheres (Figure 2B). Notably, both the model-free and model-based methods evidence that the precuneus, the superior frontal gyrus and sulcus are characterized by sentenceand paragraph-level TRFs (Figure 2A andB).

Our results differ from Lerner et al.'s in several ways. First, the acoustic TRFs are slightly more inferior with the model-based method. Second, frontal regions are detected to be associated not only with sentences and paragraphs, but also with words (consistent with [START_REF] Huth | Natural speech reveals the semantic maps that tile human cerebral cortex[END_REF][START_REF] Caucheteux | Decomposing lexical and compositional syntax and semantics with deep language models[END_REF][START_REF] Goldstein | Thinking ahead: prediction in context as a keystone of language in humans and machines[END_REF]). Given that Lerner et al's dataset is not public, it is difficult to quantify these differences and determine whether they reflect an improved sensitivity, or, more simply, inter-individual differences.

Our model-based method can, in principle, be applied to any natural stories. To test this prediction, we extend our analyses to 305 subjects listening to 4.1 hours of fifteen narratives (Figure 2C). Our model-based approach recovers the hierarchy of TRFs, and further reveals additional word-and sentence-level representations in the precuneus and prefrontal regions.

Discussion

Here, we leverage the modeling power of deep language models to show that the seminal results of Lerner et al. can be retrieved without having subjects listening to multiple scrambled stimuli. Critically, we formalize the assumptions under which 'model-based' and 'modelfree' approaches can be linked [START_REF] Lerner | Topographic Mapping of a Hierarchy of Temporal Receptive Windows Using a Narrated Story[END_REF].

Our model-based method recovers the hierarchy of TRFs evidenced by Lerner et al., in the brain of an unusually large cohort of 305 subjects. Thus, our study complements the recent work of [START_REF] Jain | Incorporating Context into Language Encoding Models for fMRI[END_REF][START_REF] Toneva | Interpreting and improving natural-language processing (in machines) with natural language-processing (in the brain)[END_REF][START_REF] Toneva | Combining computational controls with natural text reveals new aspects of meaning composition[END_REF] who predict brain responses to speech from language models input with variably-long contexts. Specifically, we show that previous model-based results unravel the same mechanisms that was previously identified with model-free approaches.

The replication is not perfect: the acoustic and word TRFs slightly differ between the two methods. This may be explained by individual subject's variability, which is only captured by the model-based approach. Further research, using the non-public data from Lerner et al. should investigate these remaining differences.

In line with previous work [START_REF] Brennan | Naturalistic sentence comprehension in the brain[END_REF][START_REF] Brennan | Hierarchical structure guides rapid linguistic predictions during naturalistic listening[END_REF][START_REF] Gauthier | Linking artificial and human neural representations of language[END_REF][START_REF] Schrimpf | Artificial Neural Networks Accurately Predict Language Processing in the Brain[END_REF], our study demonstrates that deep neural networks build constructs that predict brain activity, accurately enough to recover the hierarchy of language processing in the brain. The success of replication thus reinforces the idea that naturalistic stimuli and deep neural networks form a powerful couple to study the neural bases of language (Hamilton and Huth, 2020). five delays. Thus, f θ = f θ • g, with f θ a linear function whose parameters θ are learned, and g a temporal alignment function.

To estimate θ, we fit an 2 -penalized linear regression to predict y given g(x * ) on a training set of time samples. θ thus minimizes

argmin θ ∈Θ y train -f θ • g(x * train ) 2 + λ θ 2 ,
with λ the regularization parameter. We assess the mapping with a Pearson correlation score evaluated on the left out times samples:

R = ρ y test , f θ • g(x * test ) .
In practice, x * and g(x * ) are standardized (0-mean, 1-std) and brain signals y are scaled based on quantiles using scikit-learn Ro-bustScaler [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] with quantile range (.01, .99). We use the RidgeCV implementation of scikit-learn with a pool of twenty possible penalization parameters between 10 -3 and 10 6 . We learn f θ on 90% of the T time samples, and compute the correlation scores R on the 10% left out data. We repeat the procedure on 10 test folds using a cross-validation setting, following the KFold implementation of scikit-learn without shuffling. Finally, we average the R over the 10 folds to obtain one model-to-brain correlation score per subject, voxel and feature space x * .

C Brain parcellation

In Figure 2, we use a subdivision of Destrieux' atlas [START_REF] Destrieux | Automatic parcellation of human cortical gyri and sulci using standard anatomical nomenclature[END_REF]. Regions of more than 200 vertices are split into smaller regions, so that each region contains at most 200 vertices. Thus, from the 75 regions of Destrieux' atlas (in each hemisphere), we obtain a parcellation of 465 brain regions per hemisphere.

D Significance

In Figure 2, we test whether the model-to-brain correlations (R) are significantly different from zero. To this aim, we use a two-sided Wilcoxon test across subjects (N = 75 in Figure 2B, N = 305 in Figure 2A), corrected using False Discovery Rate (FDR) across the 465 region of interests in each hemisphere.

  mɪat ə nɒpˈə snʌw [Once upon a time. # Merlin met Arthur. #] [Arthur found the sword. # The end. #] [Arthur found the sword # The end #] [Once upon a time # Merlin met Arthur #] Merlin met Arthur # The end # Arthur found the sword # Once upon a time # Merlin end found the. sword upon time The Arthur. Once a end Arthur met ρ[y, y]

Figure 1 :

 1 Figure1: Objective and methods A. In Lerner et al.'s seminal study, each subject is presented successively with i) a 7 min long story (black), ii) the same story after its paragraphs (blue) iii) sentences iv) words (orange) or iv) acoustics (red) has been scrambled. B. For each condition, subject and voxel, the authors compute the inter-subject correlation (ISC), i.e the correlation ρ between i) the brain of the current subject y and ii) the average brain signals of the other subjects ȳ. This method allows to decompose the hierarchy of language processing in the brain, from the acoustic to the paragraph level. C. We aim to replicate the results of Lerner et al. using only the recordings induced by the regular story (black). To this aim, we scramble, not the stimulus of the subject, but the inputs of a deep language model (GPT-2). For each condition (word, sentence or paragraph), we extract the corresponding activations x * averaged over K random scrambles. We then compare the brain signals of the current subject y with the activations x * elicited by the scrambled texts, after a linear transformation f θ that maps x * onto a brain-like space. Because GPT-2 is not trained to process waveform, we use the phonemes, stresses and tones of the stimulus instead of x * for the acoustic condition.

Figure 2 :

 2 Figure 2: Results. Following Lerner et al's, a brain region is considered to process 'acoustic' level information if its acoustic score (either brains-to-brain or model-to-brain correlation) is significant (red).It is considered to process 'word'-level (yellow) if its word score is significant but not its acoustic oneand similarly for 'sentence' (green) and 'paragraph' (blue). A. Adapted from[START_REF] Lerner | Topographic Mapping of a Hierarchy of Temporal Receptive Windows Using a Narrated Story[END_REF]. Labels are based on the brains-to-brain correlation scores (Figure1B) averaged over seven subjects listening to a 7 min story. B. Labels are based on the model-to-brain scores (Figure1C), averaged over 75 subjects listening to the same 7 min story. Significance is inferred using a Wilcoxon test across subjects, corrected with False Discovery Rate (FDR) across the 465 brain regions in each hemisphere (cf. Appendix C), with a significance threshold of p < 10 -3 (cf. Appendix D). C. Same as B., but on the brain of 305 subjects listening to 4 hours of 15 audio stories (including the 7 min one). Because of the large number of subjects, the significance threshold is set to p < 10 -25 .

Figure 3 :

 3 Figure 3: Replication to two other architectures. Same as Figure 2.C but using the intermediate layers of XLNet and Distilgpt2 causal architectures (l = 4 for Distilgpt2, out of 6 layers in total and l = 8 for XLNet, out of 12 layers in total). As in Figure 2.C, the significance threshold is set to p < 10 -25 .

http://datasets.datalad.org/?dir=/labs/ hasson/narratives

http://datasets.datalad.org/?dir=/labs/ hasson/narratives

https://huggingface.co/

Appendix

To replicate Lerner et al.'s findings, we compute the model-to-brain correlation (cf. Section 2): R = ρ y, f θ (x * )) , for the acoustic, word, sentence and paragraph level respectively. Here, we provide additional details on how to extract the brain signals y and estimate the mapping function f θ in order to reproduce the experimental setting used in Section 3.

A Brain signals

Functional MRI dataset We use the fMRI recordings of the Narratives dataset (Nastase et al., 2020) 2 , a publicly available dataset gathering the brain recordings of 305 subjects listening to narratives. We use the unsmoothed version of the fMRI recordings, already preprocessed in the original dataset. As suggested in the original paper, we reject subject / narrative pairs because of noisy recordings, resulting in 617 unique (story, subject) pairs and 4.1 hours of audio stimulus in total. To replicate the results of [START_REF] Lerner | Topographic Mapping of a Hierarchy of Temporal Receptive Windows Using a Narrated Story[END_REF], we restrict the analyses to the 75 subjects listening to the 'Pieman' story (7 min long), including the seven subjects analysed in the original paper (only the data for non-scrambled stimuli are publicly available). Then, we extend the analyses to the brain recordings of 305 subjects listening to fifteen narratives (from 3 min to 57 min), from the same dataset [START_REF] Samuel | Narratives: fMRI data for evaluating models of naturalistic language comprehension[END_REF]. For both analyses, we only have access and thus use the brain recordings elicited by regular -i.e non scrambled-version of the stimuli.

B Mapping x * onto the brain

The linear function f θ maps x * onto y, the fMRI recordings of one subject at one voxel. Vector y is of length T , the number of fMRI time samples, whereas x * is of length M , the number of words (or phonemes for acoustic features) in the story. To align the two time domains, we apply the function g : R M ×D → R T ×5D that i) sums the features x * between the successive fMRI time samples, and ii) uses a Finite-Impulse Response model (FIR) with


Charlotte Caucheteux

Jean-Rémi King

THE MAPPING OF DEEP LANGUAGE MODELS ON BRAIN RESPONSES PRIMARILY DEPENDS ON THEIR PERFORMANCE A PREPRINT

Keywords: Natural Language Processing, Encoding, Functional Magnetic Resonance Imaging, Magnetoencephalography

Brains and algorithms partially converge in natural language processing

Charlotte Caucheteux, Jean-Rémi

Introduction

Deep learning has recently made remarkable progress in harnessing capabilities hitherto considered unique to the human species (4; 5; 6). In particular, language transformers demonstrate unprecedented completion, translation, and summarization abilities (7; 8; 9; 10). Do these algorithms process words and sentences like the human brain? Preliminary evidence suggests that they might. First, word embeddings -high dimensional dense vectors trained to predict lexical neighborhood (11; 12; 13; 14) -have been shown to linearly map onto the brain responses elicited by words presented either in isolation (15; 16; 17) or within narratives (18; 19; 20; 21; 22; 23). Second, the "contextualized" activations of language transformers improve the precision of this mapping, especially in the prefrontal, temporal and parietal cortices (24; 25; 26). Third, specific computations of deep language models, such as the estimations of word surprisal (i.e. the probability of a word given its context) and the parsing of syntactic constituents have been shown to correlate with evoked related potentials (27; 28; 29; 30). However, the comparison between deep language models and the brain is fragmentary. First, most studies map the highdimensional activations of deep language models onto fMRI: yet, these slow brain signals are unable to determine the sequence of brain representations elicited as the sentence unfolds. Second, past studies are based on i) a small number of subjects and ii) on a small set of language models varying in dimensionality, architecture, training objective, and training corpus. These computational differences prevent a formal comparison between the brain and these algorithms.

Third, the mapping between brains and algorithms could be driven by factors largely independent from language Figure 1: Hypotheses and Methods. A. The three panels represent three hypotheses on the link between language transformers and the brain. Each dot represents one hypothetical language transformer. Language transformers can be considered to converge to brain-like computations if their language performance (x-axis: i.e. top-1 accuracy at predicting a word from its previous context) correlates with their ability to map onto brain responses to the same stimuli (i.e. y-axis: brain score), and vice versa for a divergence hypothesis. High-dimensional neural networks can, in principle, capture relevant information (1; 2), and thus lead to a fortunate similarity with brain responses. B. Using fMRI and MEG recordings in the same subjects (3), we compare both the language performance and the brain-mapping of the layer-wise activations ('embeddings') extracted from a large variety of language transformers. C. To compute the brain scores, we (1) fit a linear regression W from the model's activations X to predict brain responses Y and (2) evaluate this mapping with a correlation between the predicted and true brain responses to held-out sentences Y test . MEG scores and fMRI scores are computed independently.

processes: for instance, networks with random weights have been shown to significantly predict brain responses to sound, speech and language stimuli (31; 32; 33). Thus, we ask two questions: is the similarity between language algorithms and the brain driven by specific factors? If so, can these algorithms help reveal the spatiotemporal hierarchy of language computations in the human brain?

To address these issues, we analyze the brain responses of 102 healthy adults, scanned with both functional magnetic resonance imaging (fMRI) and source-localized magneto-encephalography (MEG) by Schoffelen et al (3) during two 1 h-long sessions, during which they read isolated Dutch sentences composed of 9 to 15 words. We then extract 32,400 embeddings from a variety of language transformers to compare their ability to linearly map onto these brain responses ('brain score' [START_REF] Schrimpf | Artificial neural networks accurately predict language processing in the brain[END_REF]). Finally, we assess how differences in training, architectures, and language performance independently contribute to these brain mappings (Figure 1).

The results demonstrate that the similarity between language models and the brain primarily depends on their language ability (i.e. accurately predicting words from the context).

Results

Intermediate layers predict brain responses best.

To evaluate how the activations (X) of a deep language model map onto the brain (y) in response to the same sentences, we fit, within each subject, an 2 -penalized linear regression (W ) to predict single-sample fMRI and MEG responses for each voxel/sensor independently. We assess the accuracy of this mapping with a Pearson R correlation (hereafter referred to as 'brain score' [START_REF] Daniel Lk Yamins | Performance-optimized hierarchical models predict neural responses in higher visual cortex[END_REF]) between true and predicted brain activations on held-out recordings of distinct sentences, using a five-fold cross-validation. Except if stated otherwise, we report the average brain scores across all voxels (or across MEG sensors and time samples) in the text, and refer the reader to the figures for a more complete description.

Finally, we assess the statistical significance of these (average or single-voxel/channel) brain scores with a two-sided Wilcoxon test across subjects. Average MEG scores (across subjects, time, and channels) of each of the embeddings (dots) extracted from 18 causal architectures, separately for the input layer (word embedding, green) and the middle layers (red). C. Zoom of (B), focusing on the best neural networks (i.e. accuracy >35 %), revealing a slight plateau and divergence of the middle and input layers, respectively. D. Permutation importance reveals how each property of the language transformers specifically contribute to the brain scores (∆R). All properties significantly contribute to the brain scores (∆R > 0, all p < 0.0001 across subjects). Ordered pairwise comparisons of the permutation scores are marked with a star ('*': p < .05, '**': p < .01, '***': p < .001). E-H. Same as A-D, but evaluated on fMRI recordings. All error bars are the 95% confidence intervals across subjects.

We evaluate the brain scores of 32 transformer architectures (varying from 4 to 12 layers, each ranging from 128 to 512 dimensions, and each benefiting from 4 to 8 attention heads), trained on the same Wikipedia dataset either with a 'causal' language modeling (CLM) or a 'masked' language modeling task (MLM). For each architecture, we input the model with the sentence read by the subjects in the MEG or fMRI scanner, extract the activations from every layer, and, finally, compute the corresponding fMRI and MEG scores.

The brain scores of all these trained language models are significantly above chance (all p < 10 -9 , Figure 2A andE). As detailed in supplementary analyses, the modest values of these brain scores reflect the notoriously high level of noise in single-sample single-voxel/channel neuroimaging data. Indeed, fMRI and MEG scores reach R = .048 and R = .041, respectively, for the best layer of a typical 12-layer CLM, which is close to and even exceeds the noise ceiling (fMRI: R = .060, MEG: R = .020, Figure S4).

Overall, the brain scores vary as a function of the relative depth of the embedding within the transformer. Specifically, both MEG and fMRI scores follow an inverted U-shaped pattern across layers for all architectures (Figure 2A andE): middle layers 1 systematically outperform output (fMRI: ∆R = .011 ± .001, p < 10 -18 , MEG: ∆R = .003 ± .0005, Figure 3: The similarity between language transformers and the brain reveals the spatio-temporal hierarchy of language representations. Lexical and compositional representations can be isolated from the word embedding (green) and one middle layer (red) of a typical language transformer (here, a 12-layer causal transformer). To account for low-level visual representations, we also compute the brain scores of a convolutional neural network trained on character recognition (blue). A. Mean (across subjects) fMRI scores obtained with the convolutional neural network (blue), the word embedding layer (green), and the ninth layer of a 12-layer transformer (red). All colored regions display significant fMRI scores across subjects after FDR correction for multiple comparisons. p < 10 -13 ) and input layers (fMRI: ∆R = .031 ± .001, p < 10 -18 , MEG: ∆R = .009 ± .001, p < 10 -17 ). This result confirms that the intermediary representations of deep language transformers are more "brain-like" than input and output layers [START_REF] Toneva | Interpreting and improving natural-language processing (in machines) with natural language-processing (in the brain)[END_REF].

2.2

The emergence of brain-like representations predominantly depends on the network's ability to predict missing words.

The above findings result from trained neural networks. However, recent studies suggest that random (i.e. untrained) networks can significantly map onto brain responses (31; 32; 33). To test whether brain mapping specifically depends on the language proficiency of the model, we assess the brain scores of each of the 32 architectures trained with 100 distinct amounts of data. For each of these training steps, we compute the top-1 accuracy of the model at predicting masked or incoming words from their contexts. This analysis thus results in 32,400 embeddings, whose brain scores can be evaluated as a function of language performance (Figure 2B andF).

We observe three main findings. First, random embeddings systematically lead to significant brain scores across subjects and architectures. The mean fMRI score across voxels is R = .19 ± .01, p < 10 -16 . The mean MEG score across channels and time sample is R = .18 ± .008, p < 10 -16 ). This result suggests that language transformers partially map onto brain responses independent of their language abilities. Second, brain scores strongly correlate with language accuracy in both MEG (R = .77 Pearson's correlation on average ± .01 across subjects) and fMRI (R = .57 ± .02, Figure 2B andC). The correlation is higher for middle (for fMRI: R = .81 ± 02 and MEG: R = .91 ± 01) than input (R = .39 ± 03) and output layers (R = .63 ± 03). Beta coefficients for each particular layer and architecture are displayed in Figure S5A andB. Furthermore, single-voxel analyses show that this correlation between brain score and language performance is driven mainly by the superior temporal sulcus and gyrus for the embedding layer (mean R = .52 ± .06) and is widespread for the middle layers, exceeding R = .85 correlation in the superior temporal sulcus, infero-frontal, fusiform and angular gyri (Figure S5C).

Overall, this result suggests that the better language models are at predicting words from the context, the more their activations linearly map onto those of the brain.

Third, the highest brain scores are not achieved by the very best language transformers (Figure 2C andG). For instance, CLM transformers best map onto MEG (R = .039) and fMRI (R = .056) when they reach a language performance of 43% and 32%, respectively. By contrast, the very best transformers reach a language accuracy of 46%, but have significantly smaller brain scores (Figure 2C andG).

Architectural and training factors impact brain scores too.

Language proficiency co-varies with the amount of training as well as with several architectural variables. To disentangle the contribution of each of these variables to the brain scores, we perform a permutation feature importance analysis. Specifically, we train a Random Forest estimator (35) to predict the average brain scores (across voxels or MEG sensors) of each subject independently, given the layer of the representation, the architectural properties (number of layers, dimensionality, attention head), task (CLM, MLM), amount of training (number of steps) and language performance (top-1 accuracy) of the transformer. Permutation feature importance then estimates the unique contribution of each feature in explaining the variability of brain scores across models (36; 35). The results confirm that language performance is the most important factor that drives the brain scores (Figure 2 D,H). This factor supersedes other covarying factors such as the amount of training, and the relative position of the embedding with regard to the architecture ('layer position'): ∆R = .56 ± .01 for fMRI, ∆R = .51 ± .02 for MEG. Nevertheless, these other factors contribute significantly to the prediction of brain scores (p < 10 -16 across subjects for all variables).

Overall, these results show that the ability of deep language models to map onto the brain primarily depends on their ability to predict words from the context, and is best supported by the representations of their middle layers.

2.4

The mapping between the brain and language models helps to automatically decompose the cortical hierarchy of language.

Where and when are the language representations of the brain similar to those of deep language models? To address this question, we extract the activations of the first (X CLM1 ,a.k.a "word embedding") and ninth layers (X CLM9 ) of a representative 12-layer transformer trained on causal language modeling (CLM). Unlike the 9 th layer, the word embedding layer represents each word as a unique vector independent of its context [START_REF] Mikolov | Distributed representations of words and phrases and their compositionality[END_REF]. On the contrary, deeper layers are compositional: they combine word representations to best predict incoming words and can thus capture sentence-level properties like syntax (38; 39). To control for sub-lexical features, we also extract the activations X CN N of the last layer of a convolutional neural network (CNN) trained on character recognition [START_REF] Baek | What is wrong with scene text recognition model comparisons? dataset and model analysis[END_REF] and input with the image of each word. By definition, such CNN is devoid of context and is thus unable to capture the meaning of words.

Consequently, we hereafter refer to these different sets of activations as visual, lexical and compositional embeddings (Figure 3A).

In fMRI, the brain scores of the visual embedding peak in the early visual cortex (V1) (mean brain scores across voxels: R = .022 ± .003, p < 10 -11 ). By contrast, the brain scores of lexical embedding peak in the left superior temporal gyrus (R = .052 ± .004, p < 10 -13 ) as well as in the inferior temporal cortex and middle frontal gyrus (R = .053 ± .003, p < 10 -15 ) and are significant across the entire language and reading network (Figure 3B). Finally, the brain scores of the compositional embedding are significantly higher than those of lexical of embeddings in the superior temporal gyrus (∆R = .012 ± .001, p < 10 -16 ), the angular gyrus (∆R = .010 ± .001, p < 10 -16 ), the infero-frontal cortex (∆R = .016 ± .001, p < 10 -16 ) and the dorsolateral prefrontal cortex (∆R = .012 ± .001, p < 10 -13 ). While these effects are lateralized (left hemisphere versus right hemisphere: ∆R = .010 ± .001, p < 10 -14 ), they are significant across a remarkably large number of bilateral areas (Figure 3B).

Overall, these results confirm that trained deep neural networks linearly map onto the brain (41; 26). However, the sequence of representations underlying these shared representations remains unknown.

2.5

The model-to-brain mapping reveals the unfolding of language representations over both time and space.

To characterize the unfolding of brain responses over both time and space, we perform the same analysis using sourcelocalized MEG recordings. The resulting brain scores are consistent with -although less spatially precise than -the fMRI results (Figure 3C, average brain score between 0 and 2 s). For clarity, Figure 3D and Video 2 (SI.6.4) plot the gain in MEG scores: i.e. the difference of prediction performance between i) word and visual embeddings (green) and ii) the difference between compositional and word embedding (red). The brain scores of the visual embedding peak around 100 ms in V1 (R = .008 ± .002, p < 10 -3 ), and rapidly propagate to higher-level areas (Figure 3D, Video 2, SI.6.4). The gain achieved by the word embedding can be observed in the left posterior fusiform gyrus around 200 ms and peaks around 400 ms and in the left temporal and frontal cortices. Finally, the gain achieved by the compositional embedding is observed in a large number of bilateral brain regions, and peaks after ≈ 1s (Figure 3C andD).

After that period, brain areas outside the language network, such as area V1, appear to be better predicted by word and compositional embeddings than by visual ones (e.g between visual and word in V1: ∆R = .016 ± .002, p < 10 -10 ).

These effects could thus reflect feedback activity [START_REF] Seydell-Greenwald | Spoken language comprehension activates the primary visual cortex[END_REF] and explain why the corresponding fMRI responses are better accounted for by word and compositional embeddings than by visual ones.

Together with Supplementary Figure S5, these results show with unprecedented spatio-temporal precision, that the brain-mapping of our three representative embeddings automatically recovers the hierarchy of visual, lexical, and compositional representations of language in each cortical region.

Discussion

Do deep language models and the human brain process sentences in the same way? Following a recent methodology (34; 43; 44; 45; 46; 47; 48; 49; 31; 41; 26; 31), we address this issue by evaluating whether the activations of a large variety of deep language models linearly map onto those of 102 human brains. Each subject was recorded by Schoffelen and colleagues (3) with both MEG and fMRI, while they read isolated Dutch sentences composed of 9 to 15 words. Our study provides two main contributions.

Language performance is the primary factor that drives brain-mapping. First, not only do language transformers linearly map onto brain responses (24; 25; 26; 23; 27; 30; 32), but this property primarily depends on their language performance: i.e. whether these models accurately predict words from their context. Our analysis indicates that language performance is the most contributing factor explaining the variability of brain scores across models, spanning 32 architectures, two tasks, and 100 training steps (Figure 2D andH). Overall, deep language transformers thus appear to mainly converge to brain-like representations during their training.

Language transformers help decompose the cortical hierarchy of language in space and time. Second, our model comparison decomposes the visual, lexical, and compositional representations in the cortex. Whereas the areas involved in language processing are well known (50; 51; 52; 41; 53), the precise nature, format, and dynamics of their lexical and compositional representations remain largely unknown (54; 55; 52). Here, we track these hierarchical representations with unprecedented spatio-temporal precision. Early visual responses (<150 ms) are quasi-entirely accounted for by visual embeddings, and then transmitted to the posterior fusiform gyrus, which switches from visual to lexical representations around 200 ms (Video 2 and SI.6.4). This finding strengthens the claim that this area is responsible for orthographic and morphemic computations (51; 56; 57). Then, around 400 ms, lexical embeddings predict a large fronto-temporo-parietal network which peaks in the left temporal gyrus; these representations are then maintained for several seconds (15; 24; 26; 17). This result not only confirms the wide spread distribution of meaning in the brain [START_REF] Huth | Natural speech reveals the semantic maps that tile human cerebral cortex[END_REF], but also reveals its remarkably long-lasting nature.

Finally, compositional embeddings peak in the brain regions associated with high-level language processing such as the infero-frontal and anterior temporal cortices as well as the superior temporal cortex and the temporal-parietal junction (58; 52; 28). We confirm that these left-lateralized representations are significant in both hemispheres (59; 60).

Critically, MEG suggests that these compositional effects become dominant and clearly bilateral long after word onset (>800 ms). This surprisingly late response may be due to the nature of the sentences, whose complex syntactic structure may slow down compositional computations.

Overall, our results suggest a convergence between the brain and language transformers. However, several factors qualify this conclusion.

Brain scores are limited by the signal-to-noise ratio. The mapping between these models and brain recordings is low. This phenomenon is expected: i) neuroimaging is notoriously noisy and ii) we analyze and model here singlesample responses of single-voxel/sensor. However, the resulting brain scores are i) highly significant (all p < 10 -9 on average across both all fMRI voxels and MEG sensors), and ii) in the same order of magnitude to our noise ceilings (Figure S4) as well as previous reports (e.g. [START_REF] Huth | Natural speech reveals the semantic maps that tile human cerebral cortex[END_REF], before correcting for the noise ceiling). Besides, we generally report brain scores averaged across all voxels or MEG channels, even though many brain areas do not strongly respond to language (S4). Yet, brain scores often reach R > .10 in the brain areas associated with language (Figure 3). Critically, the core of our study is the link between brain scores and language performance. This effect is very strong: the correlation between the language performance and brain scores is above R= .90 for MEG and R= .80 for fMRI (Figure S5).

Language performance is not the only variable modulating brain scores. Permutation feature importance shows that several factors such as the amount of training and the architecture significantly impact brain scores. This finding contributes to a growing list of variables that lead deep language models to behave more-or-less similarly to the brain.

For example, Hale et al. [START_REF] Hale | Finding syntax in human encephalography with beam search[END_REF] show that the amount and the type of corpus impact the ability of deep language parsers to linearly correlate with EEG responses. The present work complements this finding by evaluating the full set of activations of deep language models. It further demonstrates that the key ingredient to make a model more brain-like is, for now, to improve its language performance. This conclusion, however, should be qualified, because the brain scores of some layers of the very best models tend to ultimately decrease with language performance, especially in fMRI (Figure 2G). We speculate that this unexpected phenomenon rises because transformers may over-specialize in an objective that differs from the human brain's: predicting a word from its context, as opposed to generating the meaning of a sentence.

Structure of language is different from that of the brain. The training and the architecture of transformers [START_REF] Vaswani | Attention is all you need[END_REF] are in many ways not biologically plausible. On the one hand, the brain i) is a recurrent architecture ii) is trained on a relatively small amount of grounded sentences, and iii) presumably computes prediction errors at each level of the language hierarchy [START_REF] Friston | The free-energy principle: a unified brain theory[END_REF]. On the other hand, transformers are i) feedforward neural networks, ii) trained on huge but strictly textual corpora [START_REF] Tom B Brown | Language models are few-shot learners[END_REF], iii) can memorize and access a very large number of words, and iv) only minimize prediction errors at their final layer. Besides, language transformers are still far from human-level performance in a variety of tasks such as dialogue, summarization, and systematic generalization (62; 63). Thus, it is all-the-more remarkable to see that such algorithms partially map onto brain responses.

Input and output layers show a limited convergence. The input and output layers converge less than the middle layers (Figure S5). Why is there such a difference? We speculate that syntactic representations may drive the convergence of the middle layers. Indeed, unlike word embeddings, middle layers have been shown to encode syntactic trees [START_REF] Christopher | Emergent linguistic structure in artificial neural networks trained by self-supervision[END_REF] and co-references (39; 64). Our supplementary analyses support this possibility: middle layers best encode syntactic features and this information varies with language performance similar to brain scores (Figure S5). Studying the precise nature of the shared representations between brains and transformers is an exciting direction for future work.

Overall, this study provides evidence of shared language representations between the adult human brain and language transformers, suggesting a partial convergence between the two systems. This result is important for three reasons.

First, it suggests that there is a limited number of -and perhaps unique -solutions to process language. Second, this convergence provides a concrete framework to understand the computational bases of language: deep language networks can be used as meaningful models of language processing only if they process language like our brain. Similarly, these models allow the community to move away from factorial design, and capitalize on the incremental properties of uncontrolled settings [START_REF] Liberty | The revolution will not be controlled: natural stimuli in speech neuroscience[END_REF]. Third, current language models remain relatively poor at general understanding, and zero-shot generalization (although see [START_REF] Tom B Brown | Language models are few-shot learners[END_REF]. Our results thus provide a stepping stone to unravel the cognitive operations specific to the human species, and, ultimately, implement them in machine learning algorithms.

Methods

We assess the similarity between (i) the activations of deep neural networks and (ii) those of the brain of 102 subjects, recorded with magneto-encephalography (MEG) and functional magnetic resonance imaging (fMRI) by Schoffelen et al. [START_REF] Schoffelen | A 204-subject multimodal neuroimaging dataset to study language processing[END_REF], when these two sets of systems are input with the same 400 isolated sentences.

Deep Neural Networks

Language Transformers

To model word and sentence representations, we trained a variety of transformers [START_REF] Vaswani | Attention is all you need[END_REF], input them with the same sentences that the subject read, and extracted the corresponding activations from each layer. We always extract activation in a "causal" way: for example, given the sentence 'THE CAT IS ON THE MAT', the brain response to 'ON' would be solely compared to the activations of the transformer input with 'THE CAT IS ON', and extracted from the 'ON' contextualized embeddings. Word embeddings and contextualized embeddings were generated for every word, by generating word sequences from the three previous sentences. We did not observe qualitatively different results when using shorter or longer contexts. It is to be noted that the sentences were isolated, and were not part of a narrative.

In total, we investigated 32 distinct architectures varying in their dimensionality (∈ [128, 256, 512]), number of layers (∈ [START_REF] Turing | Computing machinery and intelligence[END_REF][START_REF] Devlin | BERT: pre-training of deep bidirectional transformers for language understanding[END_REF][START_REF] Mikolov | Efficient estimation of word representations in vector space[END_REF]), attention heads (∈ [START_REF] Turing | Computing machinery and intelligence[END_REF][START_REF] Devlin | BERT: pre-training of deep bidirectional transformers for language understanding[END_REF]), and training task ("causal" language modeling and "masked" language modeling). While "causal" language transformers are trained to predict a word from its previous context, "masked" language transformers predict randomly masked words from a surrounding context. We froze the networks at ≈ 100 training stages (log distributed between 0 and 4,5M gradient updates, which corresponds to ≈ 35 passes over the full corpus), resulting in 3,600 networks in total, and 32,400 word representations (one per layer). The training was early-stopped when the networks' performance did not improve after 5 epochs on a validation set. Therefore, the number of frozen steps varied between 96 and 103 depending on the training length.

The algorithms were trained using XLM implementation2 (9), on the same Wikipedia corpus of 278,386,651 words (in Dutch) extracted using WikiExtractor3 and pre-processed using Moses tokenizer (66), with punctuation. We restricted the vocabulary to the 50,000 most frequent words, concatenated with all words used in the study (50,341 vocabulary words in total). These design choices enforce that the difference in brain scores observed across models cannot be explained by differences in corpora and text preprocessing.

To evaluate the language processing performance of the networks, we computed their performance (top-1 accuracy on word prediction given the context) using a test dataset of 180,883 words from Dutch Wikipedia.

Visual Convolutional Neural Network

To model visual representations, every word presented to the subjects was rendered on a gray 100 x 32 pixel background with a centered black Arial font, and input to a VGG network pretrained to recognize words from images [START_REF] Baek | What is wrong with scene text recognition model comparisons? dataset and model analysis[END_REF], resulting in an 888-dimensional embedding. This embedding was used to replicate and extend previous work on the similarity between visual neural network activations and brain responses to the same images (e.g. (34; 45; 46)).

Neuroimaging

Protocol

For all the analyses, we used the open-source dataset released by Schoffelen and colleagues (3), gathering the functional magnetic resonance imaging (fMRI) and magneto-encephalography (MEG) recordings of 204 native Dutch speakers (100 males), aged from 18 to 33 years. Here, we focused on the 102 right-handed speakers who performed a reading task while being recorded by a CTF magneto-encephalography (MEG) and, in a separate session, with a SIEMENS Trio 3T Magnetic Resonance scanner [START_REF] Schoffelen | A 204-subject multimodal neuroimaging dataset to study language processing[END_REF].

Words (in Dutch) were flashed one at a time with a mean duration of 351 ms (ranging from 300 to 1400 ms), separated with a 300 ms blank screen, and grouped into sequences of 9 -15 words, for a total of approximately 2,700 words per subject. Sequences were separated by a 5 s-long blank screen. We restricted our study to meaningful sentences (400 distinct sentences in total, 120 per subject). The exact syntactic structures of sentences varied across all sentences.

Roughly, sentences were either composed of a main clause and a simple subordinate clause, or contained a relative clause. Twenty percent of the sentences were followed by a yes/no question (e.g. "Did grandma give a cookie to the girl?") to ensure that subjects were paying attention. Questions were not included in the dataset, and thus excluded from our analyses. Sentences were grouped into blocks of five sequences. This grouping was used for cross-validation to avoid information leakage between the train and test sets.

Magnetic Resonance Imaging (MRI)

Structural images were acquired with a T1-weighted magnetization-prepared rapid gradient-echo (MP-RAGE) pulse sequence. The full acquisition details, available in (3), are summarized here simplicity: TR=2,300 ms, TE=3.03 ms, 8 degree flip-angle, 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. Structural images were defaced by Schoffelen and colleagues. Preprocessing of the structural MRI was performed with Freesurfer (67), using the recon-all pipeline and a manual inspection of the cortical segmentations, realigned to 'fsaverage'. Region-of-interest analyses were selected from the PALS Brodmann's Area atlas [START_REF] David | A population-average, landmark-and surface-based (pals) atlas of human cerebral cortex[END_REF] and the Destrieux atlas [START_REF] Destrieux | Automatic parcellation of human cortical gyri and sulci using standard anatomical nomenclature[END_REF].

Functional images were acquired with a T2 * -weighted functional echo-planar blood oxygenation level-dependent (EPI-BOLD) sequence. The full acquisition details, available in (3), are summarized here for simplicity: TR=2.0 seconds, TE=35ms, flip angle=90 degrees, anisotropic voxel size=3.5×3.5×3.0 mm extracted from 29 oblique slices. fMRI was preprocessed with fMRIPrep with default parameters [START_REF] Esteban | fmriprep: a robust preprocessing pipeline for functional mri[END_REF]. The resulting BOLD times series were detrended and de-confounded from 18 variables (the 6 estimated head-motion parameters (trans x,y,z , rot x,y,z ) and the first 6 noise components calculated using anatomical CompCorr (71) and 6 DCT-basis regressors using nilearn's clean_img pipeline and otherwise default parameters [START_REF] Abraham | Machine learning for neuroimaging with scikit-learn[END_REF]. The resulting volumetric data lying along a 3mm "line" orthogonal to the mid-thickness surface were linearly projected to the corresponding vertices. The resulting surface projections were spatially decimated by 10, and are hereafter referred to as voxels, for simplicity. Finally, each group of 5 sentences was separately and linearly detrended. It is noteworthy that our cross-validation never splits such groups of five consecutive sentences between the train and test sets. Two subjects were excluded from the fMRI analyses because of difficulties in processing the metadata, resulting in 100 fMRI subjects.

Magneto-encephalography (MEG)

The MEG time series were preprocessed using MNE-Python and its default parameters except when specified [START_REF] Gramfort | Mne software for processing meg and eeg data[END_REF].

Signals were band-passed filtered between 0.1 and 40 Hz filtered, spatially corrected with a Maxwell Filter, clipped between the 0.01 st and 99.99 th percentiles, segmented between -500 ms to +2,000 ms relative to word onset and baseline-corrected before t=0. Reference channels and non-MEG channels were excluded from subsequent analyses, leading to 273 MEG channels per subject. We manually co-referenced (i) the skull segmentation of subjects' anatomical MRI with (ii) the head markers digitized before MEG acquisition. A single-layer forward model was generated with the Freesurfer-wrapper implemented in MNE-Python [START_REF] Gramfort | Mne software for processing meg and eeg data[END_REF]. Due to the lack of empty-room recordings, the noise covariance matrix used for the inverse operator was estimated from the zero-centered 200 ms of baseline MEG activity preceding word onset. Subjects' source space inverse operators were computed using a dSPRM. The average brain responses displayed in Figure 1D were computed as the square of the average evoked related field across all words for each subject separately, averaged across subjects, and finally divided by their respective maxima, to highlight temporal differences.

Video 1 displays the average sources without normalization (SI.6.4). Seven subjects were excluded from the MEG analyses because of difficulties in processing the metadata, resulting in 92 usable MEG recordings.

Noise Ceiling: Brain → Brain mapping

To estimate the amount of explainable signal in each MEG and fMRI recording, we trained and evaluated, through cross-validation, a linear mapping model W to predict the brain responses of a given subject to each sentence Y from the aggregated brain responses of all other subjects who read the same sentence X. Specifically, five cross-validation splits were implemented across 5-sentence blocks with scikit-learn 'GroupKFold' [START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF]. For each word of each sentence i, all but one subject who read the corresponding sentence were averaged with one another to form a template brain response:

x i ∈ R n with n the number of MEG channels or fMRI voxels, as well as a target brain response y i ∈ R n corresponding to the remaining subject. X and Y were normalized (mean=0, std=1) across sentences for each spatio-temporal dimension, using a robust scaler clipping below and above the 0.01 st and 99.99 th percentiles, respectively. A linear mapping W ∈ R n×n was then fit with a ridge regression to best predict Y from X on the train set:

W = (X T train X train + λI) -1 X T train Y train (1)
with λ the l2 regularization parameter, chosen amongst 20 values log-spaced between 10 -3 and 10 8 with nested leave-one-out cross-validation for each dimension separately (as implemented in [START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF]). Brain predictions Ŷ = W X were evaluated with a Pearson correlation on the test set:

R = Corr(Y test , Ŷtest ) (2)
For the MEG source noise estimate, the correlation was also performed after source projection:

R = Corr(KY test , K Ŷtest ) (3)
with K ∈ R n×m the inverse operator projecting the n MEG sensors onto m sources. Correlation scores were finally averaged across cross-validation splits for each subject, resulting in one correlation score ('brain score') per voxel (or per MEG sensor/time sample) per subject.

Brain score and similarity: Network → Brain mapping

To estimate the functional similarity between each artificial neural network and each brain, we followed the same analytical pipeline used for noise ceiling, but replaced X with the activations of the deep learning models. Specifically, using the same cross-validation, and for each subject separately, we trained a linear mapping W ∈ R o,n with o the number of activations, to predict brain responses Y from the network activations X. X was normalized across words (mean=0, std=1).

To account for the hemodynamic delay between word onset and the BOLD response recorded in fMRI, we used a finite impulse response (FIR) model with five delays (from 2 to 10 seconds) to build X * from X. W was found using the same ridge regression described above, and evaluated with the same correlation scoring procedure. The resulting brain correlation scores measure the linear relationship between the brain signals of one subject (measured either by MEG or fMRI) and the activations of one artificial neural network (e.g a word embedding). For MEG, we simply fit and evaluated the model activations X at each time sample independently.

In principle, one may orthogonalize low-level representations (e.g. visual features) from high-level network models (e.g. language model), to separate the specific contribution of each type of model. This is because middle layers have access to the word-embedding layer, and can, in principle, simply copy some of its activations. Similarly, word embedding can implicitly contain visual information: e.g. frequent words tend to be visually smaller than rare ones. In our case, however, the middle layers of transformers were much better than word embeddings, which were much better than visual embeddings. To quantify the gain ∆R achieved by a higher-level model M 1 (e.g. the middle layers of a transformer) and a lower level model M 2 (e.g. a word embedding) we thus simply compared the difference of their encoding scores:

∆R M1 = R M1 -R M2 (4

Permutation feature importance

To systematically quantify how the architecture, language accuracy, and training of the language transformers impacted their ability to linearly map onto brain activity, we fitted, for each subject separately, a Random Forest across the models' properties to predict their brain scores, using scikit-learn's RandomForest (35; 36). Specifically, we input the following features to the random forest: the training task (causal language modeling "CLM" vs. masked language modeling "MLM"), the number of attention heads ∈ [START_REF] Turing | Computing machinery and intelligence[END_REF][START_REF] Devlin | BERT: pre-training of deep bidirectional transformers for language understanding[END_REF], the total number of layers ∈ [START_REF] Turing | Computing machinery and intelligence[END_REF][START_REF] Devlin | BERT: pre-training of deep bidirectional transformers for language understanding[END_REF][START_REF] Mikolov | Efficient estimation of word representations in vector space[END_REF], dimensionality ∈ [128, 256, 512], training step (number of gradient updates, ∈ [0, 4.5M ]), language modeling accuracy (top-1 accuracy at predicting a masked word) and the relative position of the representation (a.k.a 'layer position', between 0 for the word-embedding layer, and 1 for the last layer). The performance of the Random Forest was evaluated for each subject separately with a Pearson correlation R using five-split cross-validation across models.

"Permutation feature importance" summarizes how each of the covarying properties of the models (their task, architecture, etc.) specifically impacts the brain scores (35). Permutation feature importance was implemented with scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF] and is summarized with ∆R: the decrease in R when shuffling one feature (using 50 repetitions). For each subject, we reported the average decrease across the cross-validation splits (Figure 2). The resulting scores (∆R) are expected to be centered around 0 if the corresponding feature does not impact the brain scores , and positive otherwise.

Population statistics

To estimate the robustness of our results, we systematically performed second-level analyses across subjects. Specifically, we applied Wilcoxon signed-rank tests across subjects' estimates to evaluate whether the effect under consideration was systematically different from the chance level. The p-values of individual voxel/source/time samples were corrected for multiple comparisons, using a False Discovery Rate (Benjamini/Hochberg) as implemented in MNE-Python ((73)).

Error bars and ± refer to the standard error of the mean (SEM) interval across subjects.

Brain parcellation

In Section 2.2, Section 2.5 and Figure 3, we focus on particular regions of interest using the Brodmann's areas from the PALS parcellation of freesurfer 4 . The superior temporal gyrus (BA22) is split into its anterior, middle and posterior parts to increase granularity. For clarity, we rename certain areas as specified in the table below.

Anatomical and temporal characteristics of average brain responses to reading.

When and where do textual sentences elicit brain activity? As expected (54; 51; 74; 52), average fMRI and MEG responses to written words peak in a distributed and bilateral cortical network, including the primary visual cortex, the left fusiform gyrus, the supra-marginal, and the superior temporal cortices, as well as the motor, premotor and infero-frontal areas (Figure S4A). MEG source reconstruction, based on structural MRI and minimum norm estimates, further clarifies the dynamics of this cortical network: on average, word onset elicits a series of brain responses originating in V1 around ≈100 ms and continuing within the left posterior fusiform gyrus around 200 ms, the superior and middle temporal gyri, as well as the pre-motor and infero-frontal cortices between 150 and 500 ms after word onset (Figure S4A, Video 1).

Noise ceilings.

To compare the brain responses to the activations of deep language models, it is important to estimate the level of signal-to-noise ratio (Figure S4). For this, we fit, for each subject separately, a "noise-ceiling" model across subjects:

for each recording of each subject and each sentence Y train , we fit a linear model W from the recordings of all other subjects who read the same sentence X train to predict each voxel and each MEG sensor at each time sample, separately.

Using a cross-validation scheme across sentences, we then evaluate the Pearson correlation R between ( 76)). More precisely, we fit and evaluate an 2 -penalized linear model to predict each of these features given the transformer's embedding and plot this decoding performance as a function of the language performance of the model (Figure S6). While the word embedding and middle layers similarly predict word-level features (word length and part-of-speech of the word), the two high-level syntactic features (number of open and pending nodes) are better predicted by the middle layers of transformers. Finally, the decoding performance of the two syntactic features varies with the layer and the performance, in a manner strikingly similar to the brain score. These analyses suggest that middle layers are more "brain-like" than extremity layers because they learn to encode abstract linguistic properties like syntax.

Video materials

Below the captions for the two videos provided in supplementary. Figure S6: What linguistic information drives the brain score? A. From the stimulus, we compute three linguistic features: the part-of-speech of the words (i) (as given by Spacy), and two higher-level syntactic features: the number of pending nodes (ii) and open nodes (iii). These two syntactic features are derived from the constituency trees of the sentences, following [START_REF] Nelson | Neurophysiological dynamics of phrase-structure building during sentence processing[END_REF]. B-D. A 2-penalized linear regression is fit to predict the three linguistic features from the word embeddings (green), and middle layers (red) of the causal models studied in Figure 2B. The decoding performance is reported on the y-axis (accuracy at predicting the part-of-speech for B, r-squared for C, D and E). E. MEG scores (averaged across sensors and time) of the embeddings given their language modeling performance (top-1 accuracy at predicting the next word, Figure 2B).

Figure 2 :

 2 Figure 2: Language transformers tend to converge towards brain-like representations. A. Bar plots display the average MEG score (across time and channels) of six representative transformers varying in tasks (causal vs masked language modeling) and depth (4-12 layers). The green and red bars correspond to the word-embedding and middle layers, respectively. The star indicates the layer with the highest MEG score. B.Average MEG scores (across subjects, time, and channels) of each of the embeddings (dots) extracted from 18 causal architectures, separately for the input layer (word embedding, green) and the middle layers (red). C. Zoom of (B), focusing on the best neural networks (i.e. accuracy >35 %), revealing a slight plateau and divergence of the middle and input layers, respectively. D. Permutation importance reveals how each property of the language transformers specifically contribute to the brain scores (∆R). All properties significantly contribute to the brain scores (∆R > 0, all p < 0.0001 across subjects). Ordered pairwise comparisons of the permutation scores are marked with a star ('*': p < .05, '**': p < .01, '***': p < .001). E-H. Same as A-D, but evaluated on fMRI recordings. All error bars are the 95% confidence intervals across subjects.

  Figure3: The similarity between language transformers and the brain reveals the spatio-temporal hierarchy of language representations. Lexical and compositional representations can be isolated from the word embedding (green) and one middle layer (red) of a typical language transformer (here, a 12-layer causal transformer). To account for low-level visual representations, we also compute the brain scores of a convolutional neural network trained on character recognition (blue). A. Mean (across subjects) fMRI scores obtained with the convolutional neural network (blue), the word embedding layer (green), and the ninth layer of a 12-layer transformer (red). All colored regions display significant fMRI scores across subjects after FDR correction for multiple comparisons. B. The temporal resolution of MEG allows to precisely track the unfolding of the three types of representation over time. In B, the mean MEG scores averaged across all time samples and subjects. C. Left: mean MEG scores averaged across all sensors. Right: mean MEG gains averaged across all sensors: i.e. green: [word embedding] -[visual embedding]; red: [compositional embedding] -[word embedding]). D. Mean MEG gains in four regions of interest. For a whole-brain depiction of the MEG gains, see Video 2.

) 4 . 4 . 1

 441 Convergence analysisAll neural networks but the visual CNN were trained from scratch on the same corpus (cf. section 4.1.1). We systematically computed the brain scores of their activations on each subject, sensor (and time sample in the case of MEG) independently. For computational reasons, we restricted model comparison on MEG encoding scores to ten time samples regularly distributed between [0, 2]s. Brain scores were then averaged across spatial dimensions (i.e. MEG channels or fMRI surface voxels), time samples, and subjects to obtain the results in Figure2. To evaluate the convergence of a model, we computed, for each subject separately, the correlation between (1) the average brain score of each network and (2) its performance or its training step (Figure2and S5). Positive and negative correlations indicate convergence and divergence, respectively. Brain scores above 0 before training indicate a fortuitous relationship between the activations of the brain and those of the networks.

Video 1 .

 1 Anatomical and temporal hierarchy of reading. Average brain responses elicited by the onset of visual words (≈2,700 words were presented to each of the 95 subjects), as estimated with minimum source estimates (MNE) of the single-trial responses constrained by the individual subjects' anatomy (cortical surface extracted from T1 scans).

Figure S4 :

 S4 Figure S4: Noise ceilings. A. Grand average MEG source estimates to word onset (t=0) for 7 regions typically associated with reading (V1: purple, M1: green, fusiform gyrus: dark blue, supramarginal gyrus: light blue, superior temporal gyrus: orange, infero-frontal gyrus: yellow and fronto-polar gyrus: red), normalized to their peak response. Vertical bars indicate the peak time of each region. The full (not normalized) data is displayed in Video 1. B. MEG noise ceilings, approximated by predicting brain responses of a given subject from those of all other subjects. Colored lines depict the mean noise ceiling in each region of interest. The grey line depicts the best noise ceiling across sources. C. Same as (D) in sensor space. D. Noise ceiling estimates of fMRI recordings.

Figure S5 :

 S5 Figure S5: Correlation between the network's performance and brain score. A-B.Standardized beta coefficients between the language modeling performance of the network and its MEG (A) or fMRI (B) scores. For each subject, the brain scores are first scaled (0-mean, 1-std). Then, a linear regression is fit to predict the brain score (averaged across channels and time for MEG, across voxels for fMRI) of each layer of 100 networks (all 512-dimensional, with 12 layers and 8 heads) given their language performance (top-1 accuracy). The beta coefficients of the language performance are reported (y-axis). Results are consistent across 4-, 8-, and 12-layer transformers, trained on a causal (top) or masked (bottom) language modeling task. Error bars are the standard error of the mean beta coefficients across subjects. C. Pearson correlation between the performance of the 100 transformers (all 512-dimensional, with 12 layers and 8 heads) and the brain score of their word embedding (top) and ninth layer (bottom), for each voxel. Correlation scores are computed for each (subject, voxel) pair, then averaged across subjects. Only significant voxels are displayed, as assessed with a two-sided Wilcoxon test across subjects and corrected for multiple comparison using false discovery rate across voxels (threshold: .001).

  FigureS6: What linguistic information drives the brain score? A. From the stimulus, we compute three linguistic features: the part-of-speech of the words (i) (as given by Spacy), and two higher-level syntactic features: the number of pending nodes (ii) and open nodes (iii). These two syntactic features are derived from the constituency trees of the sentences, following[START_REF] Nelson | Neurophysiological dynamics of phrase-structure building during sentence processing[END_REF]. B-D. A 2-penalized linear regression is fit to predict the three linguistic features from the word embeddings (green), and middle layers (red) of the causal models studied in Figure2B. The decoding performance is reported on the y-axis (accuracy at predicting the part-of-speech for B, r-squared for C, D and E). E. MEG scores (averaged across sensors and time) of the embeddings given their language modeling performance (top-1 accuracy at predicting the next word, Figure2B). F. MEG scores of the embeddings given their ability to predict the number of open nodes.

Figure S7 :

 S7 Figure S7: Performance of the 32 transformer architectures. Best perplexity (the lower the better) and top-1 accuracy (the higher the better) of 32 transformer architectures, evaluated on a test test of 180K words from Wikipedia.Transformers are trained with a masked ('mlm') or causal ('clm') language modeling objective. They vary in their dimensionality ("Dim"), number of layers ('Layers') and number of attention heads ('Heads'). The models are trained on a set of 280K words from Wikipedia (in Dutch). The training is stopped when the perplexity on a validation set does not decrease for 5 epochs.

  1) the true brain responses of subject Y test and (2) the predicted brain responses Ŷtest = W • X test for each voxel and each MEG sensor separately. This procedure can be thought of as approximating an optimal black box: i.e. evaluating a one-hot encoder of brain responses is trained and evaluated on each element of a unique sentence. Noise ceiling peaks within the Middle layers better map onto brain responses than input and output layers. Why is there such a difference between layers? To tackle the question, we measure the level to which the 32,400 transformer embeddings linearly predict two types of linguistic features: part-of-speech (i.e a lexical feature), and the number of open and pending nodes (i.e compositional syntactic features (

	Fusiform:	0.120 ± 0.004 p < 10 -8
	Infero-frontal:	0.139 ± 0.005 p < 10 -8
	M1:	0.042 ± 0.003 p < 10 -8
	STG:	0.129 ± 0.004 p < 10 -8
	Supramarginal:	0.078 ± 0.003 p < 10 -8
	V1:	0.150 ± 0.006 p < 10 -8
	Supplementary Table 1. Average noise ceiling within each region-of-interest. Mean, standard error of the mean and
	p-values across subjects.
	6.3 Probe analysis of the language transformer.

expected language network (75) (Figure

1 F-H

). These estimates are relatively low: for example, fMRI noise ceilings reach, on average, R = 0.129 (±0.004 SEM across subjects) in the superior temporal gyrus, whereas MEG noise ceilings peak at R = 0.069 ± 0.001.

Fronto-polar cortex: 0.054 ± 0.003 p < 10 -8

For simplicity, we refer to 'middle layers' as the layer l ∈ [n layers /2, 3n layers /4], Figure

2A and E

Each algorithm was trained each on 8 GPUs using early stopping with training perplexity criteria, 16 streams per batch, 128 words per stream, epoch size of 200 000 streams, 0.1 dropout, 0.1 attention dropout, gelu activation, inverse (sqrt) adam optimizer with learning rate 0.0001, 0.01 weight decay.

https://github.com/attardi/wikiextractor

https://surfer.nmr.mgh.harvard.edu/fswiki/PALS_B12

Label

Corresponding Brodmann's areas V1 BA17 Fusiform BA37 Angular BA39 aSTG BA22-anterior mSTG BA22-middle pSTG BA22-posterior Supramarginal BA40 Infero-frontal BA44 / BA45 / BA47 Fronto-polar BA10 Temporo-polar BA38

Ethics

This study was conducted in compliance with the Helsinki Declaration. No experiments on living beings were performed for this study. These data were provided (in part) by the Donders Institute for Brain, Cognition, and Behaviour after having been approved by the local ethics committee (CMO -the local "Committee on Research Involving Human Subjects" in the Arnhem-Nijmegen region).

Acknowledgement

This work was supported by ANR-17-EURE-0017, the Fyssen Foundation, and the Bettencourt Foundation to JRK for his work at PSL. sequential recruitment of the cortical hierarchy of reading, starting from early visual cortex, moving up through the expected location of the visual word form area, and then igniting the temporal, prefrontal and parietal areas typically associated with language processing. Although these effects are bilateral, the typical left-lateralization associated with language processing can be observed.

Video 2. The main levels of the hierarchy of language revealed by deep neural networks. Single-trial encoding scores obtained for three representative embeddings reveal the types of representations that are generated within each region and at each time instant. Blue, green and red colors indicates when and where brain responses to words are specifically predicted by visual, word and compositional embeddings, respectively (a.k.a gain in brain scores).

The animated legend illustrates the same data without the anatomy: each dot corresponds to a brain source, radius corresponds to effect size (center: no effect, circle: maximum effect), and angle corresponds to the type of representation (visual, lexical or compositional). Overall, these results show when and where the brain transforms visual representations into lexical and compositional representations.