Across-speaker articulatory reconstruction from sensorimotor cortex for generalizable brain-computer interfaces

Objec&ve: Speech brain-computer interface s (BCIs) can restore speech features like arHculatory movements from brain acHvity. However, for individuals with vocal tract paralysis, lack of arHculatory movements can pose a challenge for speech BCI development. To address this challenge, our study aims at extracHng generalizable arHculatory features from a group of naHve Dutch speakers and reconstrucHng these features from brain data of a separate group of able-bodied individuals. Approach: We applied a tensor component analysis (TCA) model to extract generalisable arHculatory features from a publicly available arHculatory movement dataset. To reconstruct arHculatory features from the brain, we analyzed data of three able-bodied parHcipants P1, P2 and P3 with high-density electrocorHcography (HD-ECoG) electrode arrays implanted over the sensorimotor cortex. For each parHcipant, a separate TCA model was applied to extract neural features. A gradient boosHng regression model was used to reconstruct arHculatory features from neural features. ReconstrucHon performance was measured as Pearson’s correlaHon coefficient (PCC) between the reconstructed and the generalizable arHculatory features. Results: The extracted arHculatory features showed even contribuHons across parHcipants, indicaHng that these features captured generalizable arHculatory kinemaHc pa\erns . By using these features, we were able to reconstruct arHculatory features from brain data. PCC between the reconstructed and original arHculatory features were significantly above chance for all three parHcipants, with mean PCCs of 0.80, 0.75 and 0.76 for P1, P2 and P3 respecHvely. Significance: With the rapid development of speech BCI, our research demonstrates that speech - related arHculatory features can be restored from HD-ECoG signal using generalizable arHculatory features derived from able -bodied individuals. With the potenHal to reconstruct audio or speech - related facial movements from the re constructed arHculatory features, our framework may provide a new way for developing speech BCIs for people unable to produce mouth movements . .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint Introduc)on To produce speech, our brain coordinates the movement of vocal tract muscles , referred to as speech arHculators: lips, jaw, tongue, and larynx (Browman & Goldstein, 1992; CharHer et al., 2018; Silva et al., 2024) . Due to motor impairments caused by medical condiHons such as amyotrophic lateral sclerosis (ALS) or brain stem stroke, some people may lose the ability to control their vocal tract muscles to produce speech. The loss of speech can lead to barriers in social li fe and substanHally impact life quality (Walshe & Miller, 2011). Brain-computer interface s (BCIs) may be able to help these people by restoring their communicaHon ability. This can be done, for example, by converHng neural signals from brain areas involved in speech producHon to text on computer display or synthesized speech as has been shown in several individuals (Silva et al., 2024). However, it is unclear how generalizable these BCIs are across people. Also, these BCIs require a lot of data to train the decoding model. Most speech BCIs use signal from the ventral sensorimotor cortex (vSMC), where arHculatory kinemaHcs informaHon is strongly encoded (CharHer et al., 2018; Mugler et al., 2018). Therefore, reconstrucHng the intermediate arHculatory movement features is a logical approach for speech BCIs. Recent BCI applicaHons have demonstrated successful speech synthesis from arHculatory features reconstructed from limited neural signals (Anumanchipalli et al., 2019; Metzger et al., 2023). Accurate reconstrucHon of vocal tract and facial movement has also led to a successful implementaHon of a digital facial avatar that a person with motor impairment c ould use for enhancing communicaHon and self -expression (Metzger et al., 2023). So far, research has been focusing on connecHng brain signals to arHculaHon in the same individuals who retain their ability to control their vocal tract muscles to some extent. However, people who suffer from severe vocal tract paralysis may completely lose the ability to control vocal tract muscles and cannot produce any speech- movement or sound (Duffy, 2012). For these people, no arHculatory or acousHc data are available to train BCI decoding models . For individuals unable to produce any speech movements, transfer learning could provide a soluHon for developing arHculaHon-based speech-BCIs. Transfer learning applies knowledge learned from one dataset to a different but related dataset (Jayaram et al., 2016; Weiss et al., 2016). In the context of speech-BCIs, transfer learning can be used to extract arHculatory features shared across healthy individuals and then map them onto the brain acHvity of individuals with vocal tract paralysis. Two key research quesHons underlie the integraHon of speech-BCIs with transfer learning: (1) given individual variaHon in vocal tract anatomy, is it possible to extract generalizable arHculatory features across healthy parHcipants? (2) is it possible to reconstruct generalizable arHculatory features from brain acHvity of a different group of individuals? To address these quesHons, we developed a novel BCI framework that extracts generalizable arHculatory features across healthy parHcipants and reconstructs them from brain data collected from a separate group of parHcipants without using their own arHculatory movement data. We use d a tensor component analysis (TCA) model to extract generalizable arHculatory features (Kolda & Bader, 2009; Williams et al., 2018). We then used a gradient boosHng regression model to reconstruct the arHculatory features from high-density ECoG (HD-ECoG) data of three different parHcipants, unseen by the TCA model. We show that the arHculatory features can be reconstructed from HD-ECoG with a significant Pearson’s correlaHon with the original arHculatory features , which demonstrates the possibility of cross-parHcipants arHculatory-brain mappings . Our result highlights the potenHal of the novel framework for developing generalizable BCI models that could facilitate speech reconstrucHon in people with vocal tract paralysis. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint

Electromagne,c ar,culography (EMA) dataset We used arHculatory movement s recorded during word producHon from a publicly available Electromagne-c Ar-culography (EMA) dataset (Wieling et al., 2016). The EMA data was collected by a portable 16-channel device (WAVE, Northern Digital InC.) at a sampling rate of 100 Hz during a word producHon experiment. A total of 21 speakers of the Low-Saxon Dutch dialect and 19 speakers of the Central Dutch dialect parHcipated in the word producHon experiment. During the experiment, parHcipants performed two tasks. Task I was a picture naming task, in which parHcipants were presented with 70 pictures of objects sequenHally and were instructed to nam e the objects out loud in their dialect. Task II was a consonant-vowel-consonant ( CVC) sequences reading task, in which parHcipants were presented with 27 CVC sequences on screen sequenHally and were instructed to read them out loud in the standard Dutch dialect (Figure 1a). Therefore, the sHmuli included 97 Dutch words in total (70 object names and 27 CVC sequences). SHmuli in both tasks were repeated twice in random order. For both tasks seven EMA sensors were a\ached to jaw, lower lip (LL), right corner of the lip (right lip, RL), upper lip (UL), tongue Hp (TT), tongue body (TB) and tongue dorsal (TD) posiHons, to track movements of the corresponding arHculator along three dimensions: up-down, anterior- posterior and lem-right direcHons (Figure 2). PosiHons in the lem-right direcHon were excluded from all analyses because speech movement are typically symmetrical and provide li\le extra informaHon. We confirmed this with an analysis of variance a long each direcHon, where scores for the lem -right direcHon for tongue was very low (Table S1). The Hmestamps with pronunciaHon onset and offset of each word were included in the dataset. Figure 1. Experimental paradigms of EMA and HD -ECoG data collec>on. a. During task I of the EMA data collec8on, par8cipants were presented with pictures of objects (70 in total) and were instructed to name them. During task II, par8cipants were asked to read consonant -vowel-consonant (CVC) sequences (27 in total) presented on screen. Every picture/ CVC sequence was repeated twice in random order. This figure was adapted from images provided by Prof. Dr. Wieling, used with permission from the author. b. During the HD-ECoG data collec8on, par8cipants were presented with words from the same set of 70 object names and 27 CVC sequences (97 Dutch words in total) used in the EMA experiment. Unlike the EMA experiment, all words were presented in textual form. The task trial dura8on was 2 seconds for P1 and P2, and 1.8 seconds for P3. The dura8on of inter- trial-interval was 1.2 sec onds for P1 and P2, and 1.8 sec onds for P3. Dutch words shown in the figure: bal (ball), bijl (axe), taak (task), paak (CVC leXer sequence, not a real Dutch word), tor(beetle). .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint Figure 2. SagiAal view of loca>ons of EMA sensors. Each sensor loca8on is marked with a sphere, and the corresponding label is shown. There was one EMA sensor placed per loca8on: jaw, lower lip, right lip, upper lip, tongue 8p, tongue body and tongue dorsal loca8ons to track movements of each corresponding ar8culator in three direc8ons anterior-posterior (X-axis), up-down (Y-axis), and le[-right (Z-axis) direc8ons. In our study, we only included data collected from parHcipants speaking the Central Dutch dialect , given its greater similarity with the standard Dutch compared to the Low -Saxon Dutch dialect , as evidenced by a shorter Levenshtein distance (Heeringa, 2004). To ensure a consistent number of trials and comparable arHculator movement amplitude , we addiHonally excluded data from two Central Dutch parHcipants with missing trials , eight parHcipants with interrupted trials and one parHcipant with excepHonally large mouth movements (large amplitude of EMA signal s beyond the range of mean±2.5*std ). EMA data from eight remaining parHcipants (U1 to U8) were included in our analyses. HD-ECoG dataset Par$cipants Three male parHcipants (P1, P2, P3, age 31, 26, and 40, respecHvely) were admi\ed to the University Medical Centre Utrecht for diagnosHcs and treatment of medicaHon-resistant epilepsy or brain tumour. Two parHcipants (P1 and P2) underwent clinical ECoG implantaHon for epilepsy diagnosHc procedures and one parHcipant (P3) underwent electrode implantaHon during awake tumour resecHon surgery. All parHcipants provided wri\en informed consent for implantaHon of a high -density ECoG grid for research purposes. The protocols were approved by the Medical Ethical Commi\ee of the University Medical Cent re Utrecht in accordance with the DeclaraHon of Helsinki (20 24). All parHcipants consented to use of an high-density (HD) ECoG grid over the somatosensory cortex, which for P1 and P2 was implanted in addiHon to standard clinical grids with low density configuraHon. Here, we only focused on the HD-ECoG grids. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint Data acquisi$on All three parHcipants were implanted with a HD-ECoG grid over the sensorimotor cortex on the lem hemisphere (Figure 3) , with different configuraHons (Table 1). HD-ECoG signals were recorded at a sampling rate of 2000 Hz for P1 and P3 (Blackrock Microsystems LLC, Salt Lake City, USA), and at 2048 Hz for P2 (Micromed, Treviso, Italy). For P1 and P3 , audio recordings were recorded by a microphone connected to the neural recording system during the task. For P2, audio recordings were obtained with a separate microphone installed in the parHcipant’s room. All audio recordings were recorded at a sampling rate or 30000 Hz. For P1 and P3, the audio was recorded by Blackrock and was automaHcally synchronised with brain data. For P2, the audio was recorded by the PresentaHon® somware (Version 23.0, Neurobehavioral Systems, Inc., Berkeley, CA, www.neurobs.com) on the sHmulus presentaHon computer and was manually synchronised with brain data based on the PresentaHon log. Figure 3. Loca>ons of high-density (HD)-ECoG grid electrodes. High-density ECoG grids were implanted over the sensorimotor cortex on the le[ hemisphere. In P3, the HD grid was partly on top of another small grid blocking signals. These electrodes were removed from all analyses (marked by red crosses). The yellow lines indicate the loca8on of the central sulcus. Table 1. HD-ECoG configura>ons across par>cipants. ParHcipant Number of electrodes Inter-electrode distance (IED) Exposed diameter P1 128 3 mm 1 mm P2 32 4 mm 1 mm P3 96 3 mm 1 mm Experimental procedures PresentaHon® somware (Version 23.0, Neurobehavioral Systems, Inc., Berkeley, CA, www.neurobs.com) was used to run the experiment s for P1 and P2, and PyQT (h\ps://wiki.python.org/moin/PyQt) was used to run the experiment for P3. As in the EMA experiment, all HD-ECoG parHcipants were instructed to read out the same 97 Dutch words, including 70 object names from the picture naming task and 27 CVC sequences from the CVC sequence reading task. However, unlike the picture naming task in the EMA experimen t, all words in the HD-ECoG experiment were presented to parHcipants in textual form to avoid visual responses evoked by pictures (Figure 1b). All words were presented in random order and each word appeared on the screen twice. The duraHon of word trials and the inter-trial interval varied across parHcipants as P3 performed the task in the operaHng room , where Hme for research was limited . Thus, in P1 and P2 , the duraHon of a word trial was 2 seconds, and the duraHon of an inter-trial interval was 1.2 seconds. In P3, both the word trial and the inter-trial interval had a duraHon of 1.8 seconds. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint Data Analysis The data analysis framework (Figure 4) consists of two pipelines: the EMA analysis pipeline and the HD-ECoG analysis pipeline. Each pipeline includes three steps: data preprocessing, tensor construcHon, and feature extracHon. From these steps , we obtained arHculatory features extracted from the EMA tensor and neural features extracted from the HD -ECoG tensor. We then used these features for the arHculatory feature reconstrucHon. Figure 4. The framework for ar>culatory fea ture reconstruc>on from HD-ECoG signals. The le[ panel shows the analysis pipeline of the EMA and the right panel shows the analysis pipeline of the HD-ECoG dataset. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint EMA data preprocessing and tensor construc$on First, we idenHfied missing value s in the EMA data for each parHcipant (U1: 1.35%, U2: 1.14%, U3: 0.16%, U4: 0.41%, U5: 0.71%, U6: 0.35%, U7: 0.42%, U8: 0.51%). Missing values were imputed using linear interpolaHon available in the Pandas library of Python (h\ps://pandas.pydata.org). Then, EMA data were smoothed with a 5th -order, 20 Hz low-pass Bu\erworth filter. Amer that, EMA data in the lem-right direcHon (z-axis) were excluded. We then removed baseline from the remaining EMA data by subtracHng the mean posiHons for each combinaHon of arHculator and direcHon during a separate 10- second rest period, in which parHcipants were asked to refrain from speaking or swallowing (Wieling et al., 2016). Amer preprocessing, we segmented EMA data by using the Hme window from 0 to 1 second around the speech onset for each word trial and obtained a three-dimensional tensor : 100 x 14 x 194, which corresponds to 100 Hme points, 14 combinaHons of arHculators (7) and direcHons (2), and 194 word trials (97 unique words repeated twice) for each parHcipant. We then concatenated these tensors across parHcipants and obtained a four-dimensional tensor (8 x 100 x 14 x194). We experimented with shiming the Hme window around the speech onset (prior to speech onset) and different window sizes, but the current Hme window contained the most informaHon . HD-ECoG data preprocessing, channel selec$on and tensor construc$on The HD-ECoG data analysis pipeline was applied to each HD-ECoG parHcipant’s dataset. First, we visually inspected if there were channels with noisy or flat signals in the data. No channel was removed during this step for any of the parHcipants. We then removed 14 channels in P3 because of they lay on top of another grid and therefore were not recording signals directly from the brain Hssue (Figure 3). Next, for each remaining channel, line noise (50Hz and its harmonics) was removed with a notch filter. Common average referencing was then applied to remove common noise and trends , and HD-ECoG signals were downsampled to 512 Hz. Amer that, for each channel high frequency band (HFB; 60-170 Hz) signals were extracted in 1 Hz bins using Morlet wavelet decomposiHon and were averaged across frequency bins. HFB signals were then log-transformed and downsampled to 100 Hz to match the sampling rate of EMA . Finally, for each channel, the resulHng HFB signals were baseline-corrected by subtracHng the mean value s of a 10-second rest period before the first word trial. The speech onset moments were detected by Azure speech -to-text service (h\ps://learn.microsom.com/en -us/azure/ai-services/speech-service/fast-transcripHon- create?tabs=locale-specified) and manually corrected using Praat (Boersma, 2001) . We used the Hme window from 0.25 sec before to 0.75 sec amer speech onset to segment conHnuous recordings into word trials. The window starHng at 0.25 s before speech onset captured all informaHon of the first phoneme in the sensorimotor cortex and has been shown to produce opHmal decoding results in previous studies (Jiang et al., 2016; Ramsey et al., 2018). The 1 second window length captured the arHculatory movements of most words (95 percenHle of the distribuHon of word duraHons: 0.69 sec). Rest trials were composed of 1 second of data starHng amer the end of the preceding word trial. We idenHfied channels with significant responses to the task by comparing the mean HFB acHvity between speech and silence using a two-tailed paired t-test. Since most channels showed very significant p- values (p0.8, corresponding to a t- value of 12.44) to ensure only the most robust task-specific acHvaHon were included. The selected channels are referred to as acHve channels in subsequent analysis. Amer channel selecHon, we constructed a three-dimensional tensor (100 x the number of selected channels x 194) for each HD-ECoG parHcipant. Unlike the EMA analysis, we did not concatenate the HD-ECoG tensors across parHcipants to construct a four-dimensional tensor, as the channel dimension differed across parHcipants in number, locaHon and electrode spacing (Figure 3 and Table 1). .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint Feature extrac$on Tensor component analysis Tensor decomposiHon was applied to extract a set of matrices . Vectors within these matrices are referred to as components. Individual component captures the variability along each dimension. We analyzed the characterisHcs of these components quanHtaHvely and qualitaHvely. The tensor decomposiHon and characterisHcs analysis are referred to as tensor component analysis (TCA) (Williams et al., 2018). CANDECOMP/PARAFAC decomposiHon is one of the most commonly used tensor decomposiHon

s (Kolda & Bader, 2009). An n-dimensional tensor 𝑋 ∈ ℝ!!× !" × ⋯ × !# is decomposed into the combinaHon of R components and can be approximated by, 𝑋 ≈ ∑ 𝑎% ∘ 𝑏% ∘ 𝑐% ∘ ⋯ ∘ 𝑛% & %'( , where ∘ is the outer product , and 𝑎% ∈ 𝑖( , 𝑏% ∈ 𝑖) , … , 𝑛% (𝑟 = 1, ⋯ , 𝑅) are component s, which are the rth columns of the factor m atrices 𝐴, 𝐵, ⋯ 𝑁, respecHvely. Elements in each factor matrix are factor loadings. The CANDECOMP/PARAFAC decomposiHon procedure is omen opHmized by the alternaHng least squares method . This method has two steps. The first is to iniHalize factor matrices randomly or by certain criteria. The second step is t o opHmize one factor matrix at a Hme while keeping the others fixed. The opHmizaHon is achieved by minimizing the objecHve funcHon between the original tensor and the approximated tensor. For example, when opHmizing factor matri x 𝐴, we have the objecHve funcHon, 𝐴 ← ∑(𝑋 − ∑ 𝑎%; ∘ 𝑏% ∘ 𝑐% ∘ ⋯ ∘ 𝑛% & %'( )*+ ,%-.!/ , in which the minimizaHon is a linear -square matrix problem . All factor matrices are opHmized iteraHvely. When opHmizing each factor matrix, this procedure is iterated unHl the objecHve funcHon converges or the number of iteraHons reaches a pre-set threshold. EMA tensor decomposi&on and ar&culatory feature extrac&on We used the nonnegaHve CANDECOMP/PARAFAC decomposiHon opHmized by hierarchical alternaHng least squares to decompose the EMA tensor. The decomposiHon pipeline was implemented using the Python library tensortools (Williams et al., 2018). Amer applying the nonnegaHve CANDECOMP/PARAFAC decomposiHon, t he EMA tensor was decomposed into 29 components and four factor matrices: the parHcipant factor, Hme factor, arHculator-by-direcHon factor and word trial factor. The decomposiHon process was repeated four Hmes and resulted in four TCA runs with consistent results across the runs (Figure S3). We then checked the goodness-of-fit by inspecHng the cumulaHve explained variance of all 29 components, which was 0.89 averaged across four runs (Figure 5a). Adding more components did not substanHally increase the explained variance. To reduce components , we first selected components by the run with the highest explained variance. Since the explained variance of the first component was over 75%, we used the mean explained variance from the 2 nd to the 29 th component as the threshold to determine the number of components to select. The 10 th component was the last component with an explained variance above this threshold. Therefore, we only retained the first 10 components (Figure 5b). In the next step, our goal was to select components that generalized across all parHcipants. To measure the .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint generalizability of each component, we calculate the commonality value of each component in the parHcipant factor as following (Hurley & Rickard, 2009; Williams et al., 2018), com(𝑎%) = (∑ 𝑎%! ))0 !'( ∑ 𝑎% ! )0 !'( In which 𝑎% is the r th component in the parHcipant factor matrix and 𝑎%! is the i th element in that component. The commonality value reflects how many parHcipants have similar factor loadings, i.e., how many parHcipants contribute similarly to a component. If all parHcipants have similar factor loadings, the commonality of this component equals the number of parHcipants. In contrast, if only one parHcipant has a non-zero factor loading, the commonality of this component equals 1, indicaHng that only that parHcipant contributes to this component . We permuted elements across components 1000 Hmes to simulate the distribuHon of commonality (Figure 5c). Figure 5d shows examples of a component with high commonality and a component with low commonality. We used the mean commonality (4.57) as the threshold, as this value indicates that more than half of eight parHcipants contributed equally to the component . Using a h igher threshold yields too few features and could compromise reconstrucHon performance. By using this threshold, six components with commonality above it were selected for further analyses. We refer to selected components in these factor matrices as features. Specifically, in our reconstrucHon analysis, we refer to features from the word trial factor as the arHculatory features. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint Figure 5. Tensor component selec>on in the EMA dataset . a. Cumula8ve explained variance of all 29 tensor components averaged across 4 runs. b. Explained variance per component (star8ng from component 2 since component 1 explained approximately 80% explained variance) averaged across 4 runs. The mean explained variance was indicated by the red line and used as the threshold for component pre-selec8on. Component 10th is the last component with an explained variance above the threshold . c. The simulated distribu8on of commonality obtained by permu8ng factor loadings across the 29 components. The red doXed line indicates the mean commonality of the simulated distribu8on. d. Examples of a component with high commonality and a component with low commonality, which was 7.93 and 3.03, respec8vely HD-ECoG tensor decomposi&on and neural feature extrac&on As in the arHculatory feature extracHon, we extracted features from the HD-ECoG tensors by using the nonnegaHve CANDECOMP/PARAFAC decomposiHon opHmized by hierarchical alternaHng least squares. The HD-ECoG tensor s were decomposed into 10, 14, and 17 components for P1, P2, and P3, respecHvely, corresponding to 80% explained variance. We chose 80% explained variance as the cut- off point for model fiyng because adding more components did not increase the explained variance substanHally. To select reproducible features across repeHHons, we used Pearson’s correlaHon coefficient (PCC) between word trial factor loadings to select components . Since HD-ECoG tensor decomposiHon yield less components than EMA tensor, only one-step component selecHon is needed . First, we shuffled word trial loadings across components 1000 Hmes and calculated PCC between repeHHons for all components . Figure 6 shows the distribuHon of PCC of word trial factors between repeHHons for each HD-ECoG parHcipant. We calculated the 95th percenHle of the distribuHon and used it to select components with PCC above this threshold from the word trial factor as neural features. This resulted in 7, 6, and 6 components for parHcipants P1, P2, and P3, respecHvely. The mean PCC of extracted components between repeHHons was 0.59 for P1, 0.30 for P2 and 0.49 for P3. As in the EMA analysis, we combined vectors from the selected components to create a factor matrix for each dimension of the HD-ECoG tensor. We obtained three factor matrices: the channel factor, the Hme factor, and the word factor, which captured spaHal, temporal, and kinemaHc variaHons across words in neural acHvity, respecHvely. We also refer to components from these factor matrices as features. In our reconstrucHon analysis, features from the word trial factor matrix are referred to as the neural features. Figure 6. Selec-on of neural features in the HD-ECoG data. The panels show the distribuHons of PCC of word factors between word repeHHons across all components over 1000 shuffles for HD-ECoG parHcipants P1, P2, and P3, respecHvely. The significance threshold (α=0.05) (indicated by the green do\ed line, which represents the 95th percenHle of the distribuHon) of each distribuHon was 0.33, 0.17, and 0.31 for parHcipant P1, P2 and P3, respecHvely. We selected components with a PCC above these thresholds. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint Generalisable ar,culatory feature reconstruc,on from neural features To capture the complex relaHonships between neural data and arHculatory data, we used the gradient boosHng regression model to reconstruct the generalisable arHculatory features from the neural features (Friedman, 2001; 2002; HasHe, 2009) . This model was chosen for its robust performance on small datasets so that overfiyng can be avoided . ReconstrucHon performance was evaluated by the PCC between the reconstructed and original arHculatory features and was cross-validated across all words using the leave-one-word-out scheme . We used the permutaHon test to determine the significance of reported PCC values by shuffling word labels of all trials 1000 Hmes. PermutaHon was performed on data of repeHHon 1 and repeHHon 2 separately. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint

HD-ECoG channel selec,on For P1, P2, and P3, respecHvely, 19.5% (25/128), 75% (24/32) and 25.6% (21/82) of HD-ECoG channels showed very strong HFB increases during word trials. These channels were referred to as acHve channels (Figure 7). Figure 7. T-value maps comparing high-frequency band (HFB) responses during speech periods and silence periods. The channel colour represents t-values, as shown in the colormap. The threshold for t-values is 12.44, corresponding to a large effect size (Cohen’s d>0.8). Only channels with a t-value above the threshold were color-coded to highlight the most robust task-specific ac8va8on. In P3, the greyed-out channels were lay on top of another grid. The grey line indicates the loca8on of the central sulcus. Label “A” represents the anterior direc8on and label “P” represents the posterior direc8on. Features extracted from the EMA tensor Figure 8 shows the Hme, arHculator-by-direcHon, parHcipant and word trial factors extracted from the EMA tensor . The Hme factor shows temporal profiles of kinemaHc pa\erns shared across all arHculators. The arHculator-by-direcHon factor illustrates the strength of movements per arHculator along different direcHons. Considering both factors together reveals how the kinemaHc pa\erns evolve across Hme for each arHculator. For example, features 2 and 3 capture the tongue movements along the anterior-posterior direcHon within the first 0.5 seconds amer the speech onset. Overall, most movements occurred within 0.5 sec amer speech onset, which is the mean duraHon of word producHon. In the parHcipant factor, the commonality values of features 1-6 were 7.92, 7.88, 7.39, 7.32, 6.37, and 4.72, respecHvely, indicaHng that over half of parHcipants contribute similarly to the selected features. This result means that the extracted features generalize well across parHcipants. The word trial factor reflects the spaHotemporal arHculatory features for each word trial. In the word trial factor, the PCC between repeHHon 1 and 2 were 0.97, 0.93, 0.83, 0.95, 0.62, and 0.49 . PCC values were significant (P<0.05) for the first four features. Features 5 and 6 show a lower PCC, as well as lower commonality. Taken together, features 5 and 6 may reflect trial-specific arHculatory kinemaHcs pa\erns. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint Figure 8. Ar>culatory features extracted from EMA data using TCA. The par>cipant factor of the EMA tensor. Each bar represents the contribu8on of each par8cipant to the feature . The ar>culator-by-direc>on factor of the EMA tensor. Each bar shows movements recorded by a sensor along anterior -posterior (x) or up-down (y) direc8on. We grouped sensors into three groups by ar8culators (lip, jaw, and tongue) and colored them in blue, green and orange according to the group . The dark colors (dark blue, dark green and dark orange) represent movements along the anterior-posterior (x) direc8on, and the light colors (light blue, light green, and light orange) represent movements along the up-down (y) direc8on. The >me factor of EMA tensor. The black line indicates the speech onset, and the red line indicates the mean dura8on of spoken wo rds averaged across all trials and across all par8cipants. The word factor of the EMA tensor. We showed the Pearson’s correla8on (r) between factor loadings of repe88ons for all words, as well as the linear regression fiXed line between repe88ons (shown as the red line). .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint Features extracted from HD-ECoG tensor As in the EMA tensor, each dimension in the HD -ECoG tensor corresponds to one factor. Supplementary Figure S1 show the Hme, channel factor and word factor for P1, P2, and P3. Features in the Hme factor show the Hme profiles of basis funcHons that reflect the HFB acHvity shared across all acHve channels. In all parHcipants, we idenHfied three types of brain responses by inspecHng the Hme profiles: early response, mid response and late response. For each parHcipant, the basis funcHon reaching its peak amplitude prior to speech onset is marked as the early response; the basis funcHon showing rising acHvity amer speech onset and reaching its peak amplitude before mean speech offse t is marked as the mid response; and the basis funcHon reaching its peak amplitude around or amer the mean speech offset is marked as the late response (Figure 9). Features in the channel factor illustrate how the basis funcHons in the Hme factor are weighted for each channel , suggesHng the anatomical localizaHon of basis funcHons. In P1 and P3, channels with large weights are most located posteriorly, while in P2, channels with large weights are most located anteriorly. When examining the anatomical localizaHon of different types of brain responses, we did not observe a clear disHncHon in P1 and P3. However, in P2, the early response is evenly located across all acHve channels, while the mid and late responses are more anteriorly located. Features in the word trial factor reflect the spaHotemporal neural features for each word trial, i.e., the combinaHon of anatomical distribuHon of neural acHvaHon and Hme profiles of brain response . The PCC corresponding to early responses appeared to be lower than those corresponding with mid and late responses. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint Figure 9. Three types of brain responses during the task: early, mid response and late responses. All temporal profiles were min-max normalized to enable direct comparison across par8cipants. The early, mid, late response are shown by the blue, orange and green line, respec8vely. We iden8fied the complete set of 8me features for each response type . For the early and late response, only one 8me feature was iden8fied per par8cipant , whereas for the middle response, two 8me features were iden8fied in P2. The speech onset is shown by the black dashed line, and the mean speech offset by the red dashed line. Generalisable ar,culatory feature reconstruc,on We reconstructed arHculatory features from neural features by using the gradient boosHng regression model . Figure 10 demonstrates the reconstrucHon performance averaged across word repeHHons. Among three HD -ECoG parHcipants, P1 achieved the best reconstrucHon performance, with a mean PCC of 0.8 0 (p<0.05), while P2 and P3 achieved a mean PCC of 0.75 (p<0.05) and 0.76 (pculatory feature reconstruc>on from neural features . For each par8cipant, we calculated the real distribu8on of PCC (shown in blue). We also calculated the simulated distribu8on of PCC by shuffling word labels of all features 1000 8mes (shown in yellow). The red dashed line denotes the significance threshold (α = 0.05), and the blue dashed line denotes the mean PCC of the real distribu8on. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint

In our study, by recording brain acHvity from the vSMC using HD-ECoG grids and leveraging a publicly available EMA dataset , we demonstrated a generali zable BCI framework that reconstructed arHculatory features from three parHcipants without relying on their arHculatory or acousHc data. For individuals with severe vocal tract paralysis, this generalizable BCI framework , which is independent of individual-specific arHculatory or acousHc data, could be a promising way to restore speech from brain acHvity. Because arHculator movements were not measured in the HD-ECoG parHcipants, we used a TCA model to extract arHculatory features from a public EMA dataset . We then selected features with high commonality values, indicaHng that these arHculatory features capture arHculatory movement pa\erns shared across speakers. Earlier studies also a\empted to align arHculatory features across speakers. For example, Bocquelet and colleagues used a linear model to map the arHculatory features of individual speakers to the feature space of a reference speaker (Bocquelet et al., 2016). However, their approach can be biased by the arHculatory characterisHcs of the chosen reference, such as vocal tract shape. In contrast, our approach employed commonality to ensure that the selected features showed generalizable arHculatory kinemaHcs across healthy speakers, which is parHcularly advantageous when both acousHc and arHculatory data cannot be obtained from target users . A previous study not having arHculatory movement data available reported good results using an acousHc-to-arHculatory inversion model relying on parHcipants’ own audio data to generate arHculatory data for speech BCI development (Anumanchipalli et al., 2019). However, their approach requires parHcipants to produce intelligible speech, which is not feasible for individuals with severe vocal tract paralysis. Having extracted a set of generalizable arHculatory features across speakers, we next examined the interpretability of the extracted neural features. In all three parHcipant we idenHfied three types of brain responses: early, mid and late. Compared with previous studies, the early and late responses we observed resemble the transient response, which shows increased acHvity around speech onset/offset, and the mid response resembles the sustained response , which shows increased acHvity during the u\erance (Conant et al., 2018; Salari et al., 2018). However, we did not observe a clear disHncHon in the anatomical organizaHon among these responses, while Salari and colleagues found that transient responses were more anteriorly located and the sustained responses were more posteriorly located in their parHcipants with HD grids (Salari et al., 2018). One possible explanaHon is the relaHvely small number of acHve channels in our analysis . Beyond the anatomical organizaHon of these brain responses, we also considered their possible funcHonal roles: early and late responses may reflect the arHculator movements at the speech onset/offset, such as mouth opening or closing . AlternaHvely, the early response may also reflect movement planning (Salari et al., 2018). For the mid response, it may represent a nonspecific signal for holding the vocal tract configuraHon (Conant et al., 2018). To further understand the neural representaHon during word producHon, more HD -ECoG data will be needed. To reconstruct arHculatory features from brain acHvity, we train ed a gradient boosHng regression model using word trial factors from both the EMA and HD -ECoG TCA model s. In all three parHcipants, the reconstructed and the original arHculatory features significantly correlated with each other . However, the correlaHon coefficients were only slightly above the significance threshold (α = 0.05) , indicaHng that the arHculatory features may be similar across words . One possibility is that the arHculatory features of different words cluster closely in the low-dimensional space . Further studies should include more repeHHons per word and train a word -level classifier to evaluate the discriminability of arHculatory features. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint To the best of our knowledge, this is the first study to extract generalizable arHculatory features from healthy individuals and reconstruct them from brain acHvity of different individuals. Previous work showed that speakers share a similar state -space of arHculatory features, but only used this to reconstruct voice, demonstraHng the cross-subject arHculatory-acousHc transfer (Anumanchipalli et al., 2019) . Our study extends their work by showing the cross-subject arHculatory-neural transfer, demonstraHng that the generalizable arHculatory features can be reconstructed from brain acHvity, even if the two modaliHes are from different individuals. This finding is also consistent with prior research, where authors report structural similarity between the arHculatory features from healthy parHcipants and neural features from an individual with paralysis (Wille\ et al., 2023). Such cross- subject arHculatory-neural transfer is beneficial for speech BCI design: even for individuals with severe vocal tract paralysis, it may sHll be possible to reconstruct arHculatory movements by developing models using arHculatory data from healthy individuals. Recent advances in speech BCIs have shown the possibility of reconstrucHng speech directly from brain signal without relying on intermediate arHculatory features (Wairagkar et al., 2024; Li\lejohn et al., 2025). However, previous research suggests that arHculatory features are more robustly encoded in the vSMC than acousHcs features, and therefore can be learned faster with limited neural data (CharHer et al. , 2018; Conant et al. , 2018; Anumanchipalli et al. , 2019) . Thus, incorporaHng intermediate arHculatory features may enhance the speech decoding performance (Anumanchipalli et al., 2019) . Apart from audio synthesis, the intermediate arHculatory features can also be used to restore speech-related orofacial movements (Metzger et al., 2023). . Taken together, our results indicate that our generalizable BCI framework could be a promising approach for restoring the communicaHon for individuals with severe vocal tract paralysis, for whom arHculatory or acousHc data are not available. However, our current use of EMA data, which only contains 7 sensors a\ached to the upper vocal tract , restricts the ability to capture complete vocal tract movements. Using more complete arHculatory features that span the enHre vocal tract space may enable speech synthesis with human-like fidelity (Wu et al., 2023). Another limitaHon of our study is the vocabulary size. Our restricted vocabulary of 97 Dutch words is likely not enough to capture the full range of arHculatory gestures in natural conversaHon. Expanding to larger vocabularies and to sentences will be important for assessing reconstrucHon performance in more realisHc contexts. Finally, the proposed pipeline has only been validated offline. In the future, adapHng this framework for online use could enable real-Hme arHculatory-based speech BCIs.

Overall, we demonstrated that our proposed framework could reconstruct generalizable arHculatory features from brain acHvity from a separate group of able-bodied speakers, even when these speakers’ arHculatory or audio data were not available . Using these generalizable arHculatory features has potenHal for developing speech BCIs that can restore full communicaHon for individuals with severe vocal tract paralysis, while reducing the need for a large amount of training data.

This work is supported by Dutch Brain Interface IniHaHve (DBI2), project number 024.005.022 of the research programmed GravitaHon, which is financed by the Dutch Ministry of EducaHon, Culture, and Science (OCW) via the Dutch Research Council (NWO). .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint

Anumanchipalli, G.K., Char-er, J. & Chang, E.F. (2019) Speech synthesis from neural decoding of spoken sentences. Nature, 568, 493-498. Bocquelet, F., Hueber, T., Girin, L., Savariaux, C. & Yvert, B. (2016) Real--me control of an ar-culatory -based speech synthesizer for brain computer interfaces. PLoS computa/onal biology, 12, e1005119. Boersma, P . (2001) Praat, a system for doing phone-cs by computer. Glot. Int., 5, 341-345. Browman, C.P . & Goldstein, L. (1992) Ar-culatory phonology: An overview. Phone/ca, 49, 155-180. Char-er, J., Anumanchipalli, G.K., Johnson, K. & Chang, E.F. (2018) Encoding of ar-culatory kinema-c trajectories in human speech sensorimotor cortex. Neuron, 98, 1042-1054. e1044. Conant, D.F., Bouchard, K.E., Leonard, M.K. & Chang, E.F. (2018) Human sensorimotor cortex control of directly measured vocal tract movements during vowel produc-on. Journal of Neuroscience, 38, 2955-2966. Duffy, J.R. (2012) Motor speech disorders: Substrates, differen/al diagnosis, and management. Elsevier Health Sciences. Friedman, J.H. (2001) Greedy func-on approxima-on: a gradient boos-ng machine. Annals of sta/s/cs, 1189-1232. Friedman, J.H. (2002) Stochas-c gradient boos-ng. Computa/onal sta/s/cs & data analysis, 38, 367-378. Has-e, T. (2009) The elements of sta-s-cal learning: data mining, inference, and predic-on. Springer. Heeringa, W.J. (2004) Measuring dialect pronuncia-on differences using Levenshtein distance. Hurley, N. & Rickard, S. (2009) Comparing measures of sparsity. IEEE Transac/ons on Informa/on Theory, 55, 4723-4741. Jayaram, V., Alamgir, M., Altun, Y ., Scholkopf, B. & Grosse-Wentrup, M. (2016) Transfer learning in brain-computer interfaces. IEEE Computa/onal Intelligence Magazine, 11, 20-31. Jiang, W., Pailla, T., Dichter, B., Chang, E.F. & Gilja, V. (Year) Decoding speech using the -ming of neural signal modula-on. 2016 38th Annual Interna-onal Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). IEEE, City. p. 1532-1535. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint Kolda, T.G. & Bader, B.W. (2009) Tensor decomposi-ons and applica-ons. SIAM review, 51, 455-500. Lidlejohn, K.T., Cho, C.J., Liu, J.R., Silva, A.B., Yu, B., Anderson, V.R., Kurtz-Miod, C.M., Brosler, S., Kashyap, A.P . & Hallinan, I.P . (2025) A streaming brain-to-voice neuroprosthesis to restore naturalis-c communica-on. Nature neuroscience, 1-11. Metzger, S.L., Lidlejohn, K.T., Silva, A.B., Moses, D.A., Seaton, M.P ., Wang, R., Dougherty, M.E., Liu, J.R., Wu, P . & Berger, M.A. (2023) A high-performance neuroprosthesis for speech decoding and avatar control. Nature, 620, 1037-1046. Mugler, E.M., Tate, M.C., Livescu, K., Templer, J.W., Goldrick, M.A. & Slutzky, M.W. (2018) Differen-al representa-on of ar-culatory gestures and phonemes in precentral and inferior frontal gyri. Journal of Neuroscience, 38, 9803-9813. Ramsey, N.F., Salari, E., Aarnoutse, E.J., Vansteensel, M.J., Bleichner, M.G. & Freudenburg, Z. (2018) Decoding spoken phonemes from sensorimotor cortex with high-density ECoG grids. Neuroimage, 180, 301-311. Salari, E., Freudenburg, Z., Vansteensel, M. & Ramsey, N. (2018) Spa-al-temporal dynamics of the sensorimotor cortex: sustained and transient ac-vity. IEEE Transac/ons on Neural Systems and Rehabilita/on Engineering, 26, 1084-1092. Silva, A.B., Lidlejohn, K.T., Liu, J.R., Moses, D.A. & Chang, E.F. (2024) The speech neuroprosthesis. Nature Reviews Neuroscience, 25, 473-492. Tomasi, G. & Bro, R. (2006) A comparison of algorithms for figng the PARAFAC model. Computa/onal Sta/s/cs & Data Analysis, 50, 1700-1734. Wairagkar, M., Card, N.S., Singer-Clark, T., Hou, X., Iacobacci, C., Hochberg, L.R., Brandman, D.M. & Stavisky, S.D. (2024) An instantaneous voice synthesis neuroprosthesis. bioRxiv. Walshe, M. & Miller, N. (2011) Living with acquired dysarthria: the speaker's perspec-ve. Disability and rehabilita/on, 33, 195-203. Weiss, K., Khoshgojaar, T.M. & Wang, D. (2016) A survey of transfer learning. Journal of Big data, 3, 1-40. Wieling, M., Tomaschek, F., Arnold, D., Tiede, M., Bröker, F., Thiele, S., Wood, S.N. & Baayen, R.H. (2016) Inves-ga-ng dialectal differences using ar-culography. Journal of Phone/cs, 59, 122-143. Willed, F.R., Kunz, E.M., Fan, C., Avansino, D.T., Wilson, G.H., Choi, E.Y ., Kamdar, F., Glasser, M.F., Hochberg, L.R. & Druckmann, S. (2023) A high-performance speech neuroprosthesis. Nature, 620, 1031-1036. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint Williams, A.H., Kim, T.H., Wang, F., Vyas, S., Ryu, S.I., Shenoy, K.V., Schnitzer, M., Kolda, T.G. & Ganguli, S. (2018) Unsupervised discovery of demixed, low-dimensional neural dynamics across mul-ple -mescales through tensor component analysis. Neuron, 98, 1099-1115. e1098. Wu, P ., Li, T., Lu, Y ., Zhang, Y ., Lian, J., Black, A.W., Goldstein, L., Watanabe, S. & Anumanchipalli, G.K. (2023) Deep speech synthesis from MRI-based ar-culatory representa-ons. arXiv preprint arXiv:2307.02471. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint Supplementary materials .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint Figure S1 ECoG Components of P1, P2, and P3. From le[ to right in each panel, we showed the channel factor, the 8me factor, and the word factor for par8cipant P1, P2 and P3, respec8vely. In the channel factor, the grey lines indicate the loca8on of central sulcus. In the 8me factor, the black dash lines indicate the speech onset and the red dash lines indicate the speech offset averaged across word trials. For components with the flat line, no brain ac8vity was reflected during the corresponding period . In the word factor, the red line s are the linear regression fiXed line between feature loadings of repe88on 1 and 2. We also showed the Pearson’s correla8on (r) between feature loadings of repe88on 1 and 2 for each selected feature. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint Figure S2 Distribu>on of word produc>on during HD-ECoG data collec>on. We accumulate all word produc8on dura8on, which is the interval between speech onset and offset across all word trials for three HD -ECoG par8cipants and calculate the distribu8on. Figure S3. Similarity scores between TCA runs with different number of components. For each TCA run (blue dot), a similarity score was calculated across all TCA runs and the run with the lowest reconstrucHon error .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint Table S1. Average variance per ar-culator-by-direc-on across par-cipants . For each arHculator, the variance of movements along each direcHon: anterior -posterior (X), up-down (Y), and lem-right (Z) was calculated and average across parHcipants. Ar-culator Direc-on Anterior-posterior (x) Up-down (Y) Lej-right (Z) UL 1.41 1.62 1.52 RL 4.27 1.20 4.09 LL 6.18 4.38 3.19 JAW 2.23 2.55 1.69 TD 16.57 8.77 3.96 TB 17.22 9.46 3.93 TT 18.81 9.52 5.24 Similarity score The similarity of two fi\ed TCA models were calculated based on the angles between latent factors (Tomasi & Bro, 2006; Williams et al., 2018). For two fi\ed four-factor TCA models with R components , {𝐴, 𝐵, 𝐶, 𝐷} and {𝐴1, 𝐵1 , 𝐶1, 𝐷1}, the similarity score is, 𝑚𝑎𝑥 𝑤 ∈ Ω 1 𝑅 G H(1 − I𝜆% − 𝜆2(%)K maxN𝜆%, 𝜆2(%)OP (𝑎5𝑎1 ∙ 𝑏5𝑏1 ∙ 𝑐5𝑐1 ∙ 𝑑5𝑑1)] & %'( , where Ω denotes all possible permutaHons of all factors , and 𝑤 is a parHcular permutaHon. When calculaHng the similarity score, every Hme one component is selected from each of the four factors, each scaled to unit length and 𝜆% denotes the product of these scalings. The similarity score is averaged across all components. Amer enumeraHng all possible permutaHons, the maximal similarity score is taken as the final similarity score between these two TCA models. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted October 8, 2025. ; https://doi.org/10.1101/2025.10.08.680888doi: bioRxiv preprint