All subjects provided informed written consent to participate in the study, which was approved by the Institutional Review Boards of Columbia University (Columbia University Human Research Protection Office), New York University (NYU) Langone Medical Center (Human Research Protections) and the University of Rochester (Research Subjects Review Board).

Measuring durational variability of speech phonemes

To illustrate the variability of speech structures, we measured the durational variability of phonemes in the LibriSpeech corpus59 (Fig. 1a,b and Extended Data Fig. 1). Phoneme alignments were computed using the Montreal Forced Aligner60,61. The duration estimates for word-initial and word-final phonemes were occasionally contaminated by periods of silence that were included in the phoneme’s segmentation; therefore, we discarded word-initial and word-final phonemes from our estimates to ensure they did not inflate our variability estimates. For each phoneme, we then measured the distribution of durations across all speakers and utterances in the corpus. We then calculated the central 95% interval of this distribution and measured the ratio between the upper and lower boundaries of this interval as a measure of durational variability.

Cross-context correlation analysis

We used the ‘cross-context correlation’ to estimate context invariance from both computational models and neural data. In this analysis, we first compile the response timecourses to all segments of a given duration in a segment-by-time matrix (Fig. 2a). Each row contains the response timecourse surrounding a single segment, aligned to segment onset. Different rows thus correspond to different segments, and different columns correspond to different lags relative to segment onset. We compute a separate matrix for each of the two contexts being compared. The central segment is the same between contexts, but the surrounding segments differ.

Our goal is to determine whether there is a lag when the response is the same across contexts. We instantiate this idea by correlating corresponding columns across segment-aligned response matrices from different contexts (schematized by the linked columnar boxes in Fig. 2a). At segment onset (Fig. 2a, first box pair), the cross-context correlation should be near zero because the integration window must overlap the preceding segments, which are random across contexts. As time progresses, the integration window will start to overlap the shared segment, and the cross-context correlation should increase. If the integration window is less than the segment duration, there will be a lag in which the integration window is fully contained within the shared segment, and the response should thus be the same across contexts, yielding a correlation of 1 (Fig. 2a, second box pair).

Our stimuli enable us to investigate and compare two types of contexts: cases in which a segment is a subset of a longer segment and thus surrounded by its natural context (Fig. 1c, left) and cases in which a segment is surrounded by randomly selected segments of the same duration (Fig. 1c, right). We computed the cross-context correlation by comparing natural and randomly selected contexts as well as two different randomly selected contexts (see our previous publication7 for details), but the results were very similar when only comparing natural and random contexts. We note that any response that is selective for naturalistic structure (for example, a word) will, by definition, show a difference between natural and random contexts, and thus our paradigm will be sensitive to this change. In practice, we found that the cross-context correlation was close to the noise ceiling for segment durations of at least 333 ms (Fig. 3a,e and Extended Data Fig. 5), and this was true even when only comparing natural and random contexts. This fact demonstrates that segment durations of 333 ms produce similar responses to those of longer segments (1 s or longer) in the auditory cortex.

Estimating integration windows from computational model responses

We tested whether our approach could distinguish time-yoked versus structure-yoked integration by applying our analyses to the outputs of computational models. Our methods were similar to our neural analyses, except that in the case of computational models, the response is noise-free and we are not constrained by experiment time, and therefore can measure responses to many segments. We measured the cross-context correlation using 30 segment durations with 100, 27 s sequences per duration (segment durations of 20, 40, 60, 80, 120, 160, 200, 240, 280, 320, 400, 480, 560, 640, 720, 800, 880, 960, 1,040, 1,120, 1,200, 1,280, 1,440, 1,600, 1,760, 1,920, 2,080, 2,240, 2,400, 2,560 ms). Segments were excerpted from LibriSpeech following time compression and stretching by a factor of \(\surd 3\) as in the neural experiments (time stretching and compression were implemented using waveform synchronous overlap and add, as implemented by SOX in Python62).

To estimate the integration window, we applied our cross-context correlation analysis, which is described both in the main text and in our previous publication7 (Fig. 2a) (comparing natural and random contexts). We then calculated the peak cross-context correlation value for each segment duration as a measure of context invariance. Finally, we interpolated the correlation versus segment duration curve to determine the smallest segment duration needed to achieve a threshold cross-context correlation value (threshold set to 0.75).

STRF model

Following standard practice, our STRF model was defined by applying a linear transformation to a spectrogram representation of sound (mel spectrogram with torchaudio; sample_rate = 16,000, n_fft = 1,280, hop_length = 320, f_min = 40, f_max = 8,000, n_mels = 128, power = 2.0, norm = ‘slaney’). To make the STRF model more realistic, the STRFs were fit to approximate human intracranial responses to natural speech. The STRFs were fit using regularized regression (ridge regression) against lagged spectrogram features (tenfold cross-validation was used to select the regularization parameter63). The weights from the regression analysis define the STRF window. We used a window size of 1 s sampled at 100 Hz. A 50 ms half-Hanning window was applied to the beginning and end of each STRF to suppress edge artifacts64. The data and stimuli used to fit the STRF models have been described previously46 and consisted of 566 electrode responses from 15 patients. Models were fit to all electrodes from that prior study46. Each patient listened to 30 min of speech excerpted from a children’s storybook (Hank the Cowdog) and an instructional audio guide (four voice actors; two male, two female).

Phoneme integration model

As a simple model of structure-yoked integration, we instantiated a model that integrated phonemic features within a window whose temporal extent varied inversely with the speech rate. The features of the phoneme model were defined by 22 binary phonetic features, each indicating the presence (1) or absence (0) of a single feature (for example, manner of articulation)46. As is common, we used onset features in which the presence or absence of a feature is represented by a ‘1’ at the onset of each phoneme. As with the STRF, the window was fit to neural data to make it more similar to that from a neural experiment. The window was fit in the same way by regressing time-lagged features against the neural response, and the weights from the regression analysis defined the window. We then stretched or compressed this window by interpolation when predicting responses to stretched and compressed speech, respectively, to simulate a structure-yoked response.

DANN model

We computed integration windows from a popular DANN speech recognition model (DeepSpeech2)35,65. The DANN consists of two convolutional layers, five recurrent layers (long-short-term memory cells) and one linear readout layer, and it was trained to transcribe text from a mel spectrogram using a connection-temporal-classification loss (applied to graphemes). The model was trained using 960 h of speech from the LibriSpeech corpus using standard data augmentation techniques35 (background noise, reverberation and frequency masking) to make recognition more challenging (25 epochs; optimization was implemented in PyTorch using the Adam optimizer). Each layer is defined by a set of unit response timecourses, and we measured the integration of each unit timecourse by applying the TCI analysis to its response (in the same manner as all other models). Like most DANN models, DeepSpeech2 is acausal, but this is not problematic for measuring its integration window35 because acausality shifts the cross-context correlation to earlier lags, and our measure of context invariance (the peak correlation across lags) is invariant to shifts.

Intracranial recordings from human patientsParticipants and data collection

Data were obtained from 15 patients undergoing treatment for intractable epilepsy at the NYU Langone Hospital (six patients) and the Columbia University Medical Center (CUMC) (nine patients) (seven male, eight female; mean age, 36 years, standard deviation, 13 years). Three additional patients were tested in a follow-up experiment (all female; ages 29, 30 and 50 years) in which we measured responses to speech that was naturally faster or slower. Two of these patients were tested at the University of Rochester Medical Center (URMC), and one patient was tested at NYU. Electrodes were implanted to localize epileptogenic zones and delineate these zones from eloquent cortical areas before brain resection. NYU patients were implanted with subdural grids, strips and depth electrodes depending on the clinical needs of the patient. CUMC patients were implanted with depth electrodes. NYU patients were compensated $20 per hour. URMC patients were compensated $35 per hour. CUMC patients were not compensated owing to Institutional Review Board prohibition.

Stimuli

Segments of speech were excerpted from a recording of a spoken story from the Moth Radio Hour (Tina Zimmerman, Go In Peace). The recording was converted to mono and resampled to 20 kHz; pauses longer than 500 ms were excised. The stimuli were then compressed and stretched while preserving pitch using a high-quality speech vocoder STRAIGHT66. There were five logarithmically spaced segment durations (37, 111, 333, 1,000, 3,000 ms). Each stimulus was 27 s and was composed of a sequence of segments of a single duration. There were two stimuli per segment duration, each with a different ordering of segments. We also tested a 27 s stimulus composed of a single, undivided excerpt. Shorter segments were created by subdividing longer segments. Each stimulus was repeated twice (in one subject, the stimuli were repeated four times). In two participants, we tested an additional segment duration (9 s) and used a longer stimulus duration (45 s) with only one ordering of segments rather than two.

We conducted a subsequent experiment to test whether similar results would be observed for speech that is naturally faster or slower. Specifically, we recorded the same talker (author S.V.N.-H.) producing a common set of sentences at either a fast or slow pace. There were 28 fast sentences and 11 slow sentences, such that the total duration of the material was approximately the same (40.41 s, 41.24 s). There were 11 shared sentences that were present for both the fast and slow conditions, plus 17 additional sentences that were only tested for the fast condition (because more sentences were needed). For the 11 shared sentences, the slow sentences were 2.53 times longer than the fast sentences (3.75 s versus 1.48 s). Each sentence was root mean square-normalized after removing very low frequencies below 50 Hz (fourth-order Butterworth filter). We then generated segment sequences in the same manner as that described above (segment durations of 62.5, 125, 250, 500 ms).

Preprocessing

Our preprocessing pipeline was similar to prior studies7,67. Electrode responses were common-average referenced to the grand mean across electrodes from each subject. We excluded noisy electrodes from the common-average reference by detecting anomalies in the 60 Hz power band (measured using an IIR resonance filter with a 3 dB down bandwidth of 0.6 Hz; implemented using MATLAB’s iirpeak.m). Specifically, we excluded electrodes whose 60 Hz power exceeded five standard deviations of the median across electrodes. Given that the standard deviation is itself sensitive to outliers, we estimated the standard deviation using the central 20% of samples, which are unlikely to be influenced by outliers (we divided the range of the central 20% of samples by that which would be expected from a Gaussian distribution of unit variance). After common-average referencing, we used a notch filter to remove harmonics and fractional multiples of the 60 Hz noise (60, 90, 120, 180; using an IIR notch filter with a 3 dB down bandwidth of 1 Hz; the filter was applied forward and backward; implemented using MATLAB’s iirnotch.m).

We measured integration windows from the broadband gamma power response timecourse of each electrode. We computed broadband gamma power by measuring the envelope of the preprocessed signal filtered between 70 and 140 Hz (implemented using a sixth-order Butterworth filter with 3 dB down cutoffs of 70 Hz and 140 Hz; the filter was applied forward and backward; envelopes were measured using the absolute value of the analytic signal, computed using the Hilbert transform; implemented using fdesign.bandpass in MATLAB). We have previously shown that the filter does not strongly bias the measured integration time because the integration window of the filter (~19 ms) is small relative to the integration window of the measured cortical responses7. Envelopes were downsampled to 100 Hz. We detected occasional artifactual time points as time points that exceeded five times the 90th percentile value for each electrode (across all time points for that electrode), and we interpolated these outlier time points from nearby non-outlier time points (using ‘piecewise cubic Hermite interpolation’ as implemented by MATLAB’s interp1.m function).

As is standard, we time-locked the intracranial EEG recordings to the stimuli by either cross-correlating the audio with a recording of the audio collected synchronously with the intracranial EEG data or by detecting a series of pulses at the start of each stimulus that were recorded synchronously with the intracranial EEG data. We used the stereo jack on the experimental laptop to either send two copies of the audio or to send audio and pulses on separate channels. The audio on one channel was used to play sounds to subjects, and the audio/pulses on the other were sent to the EEG recording system. Sounds were played through a Bose Soundlink Mini II speaker (at CUMC), an Anker Soundcore speaker (at NYU) or a Genelec 8010a speaker (at CUMC). Responses were converted to units of percent signal change relative to silence by subtracting and then dividing the response of each electrode by the average response during the 500 ms before each stimulus.

We selected electrodes with a significant test–retest correlation (Pearson correlation) across the two presentations of each stimulus. Significance was measured with a permutation test, in which we randomized the mapping between stimuli across repeated presentations and recomputed the correlation (using 1,000 permutations). We used a Gaussian fit to the distribution of permuted correlation coefficients to compute small P values. Only electrodes with a highly significant correlation relative to the null were retained (P < 10−5). We also excluded electrodes for which the test–retest correlation fell below 0.05, which resulted in 132 total electrodes.

Following standard practice, we localized electrodes as bright spots on a post-operative computer tomography image or dark spots on an MRI, depending on which was available. The post-operative computer tomography or MRI was aligned to a high-resolution, pre-operative MRI that was undistorted by electrodes. Each electrode was projected onto the cortical surface computed by Freesurfer from the pre-operative MRI, excluding electrodes greater than 10 mm from the surface. We used the same correction procedure as in our prior studies7,67 to correct gross-scale errors by encouraging points that are nearby in 3D space but far apart on the 2D cortical surface (for example, two abutting gyri) to preferentially be localized to regions where sound responses are common.

Estimating integration windows from noisy neural responses

Neural responses are noisy and therefore never produce the same response, even if the response is context-invariant. Moreover, we are limited in the number of segment durations and segments that we can test owing to limited experimental time with patients. To address these challenges, we measure a noise ceiling for our cross-context correlation when the context is identical, using repeated presentations of the same segment sequence. The noise ceiling is a stimulus-dependent measure because it reflects the relative strength of the stimulus-driven and noise variance, which necessarily varies with the stimuli. Therefore, we compute a separate noise ceiling for each time lag, segment duration and speech rate. We then estimate the integration window that best predicts the cross-context correlation by finding a parametric window that best predicts the cross-context correlation pooling across all lags and segment durations. The predictions are computed by multiplying a noise-free prediction from the parametric model window by the measured noise ceiling. We have previously described and justified our parametric model in detail and have extensively tested the method, showing that it can correctly estimate integration windows from a variety of ground-truth models without substantial bias using noisy, broadband gamma responses with similar signal-to-noise ratios as those in actual neural data7. We report results from 110 electrodes (out of 132) for which the model predictions were highly significant (P < 0.001; measured using a significance test with phase-randomized predictions7).

In our original formulation, the window was parametrized using a gamma probability density function, with three parameters that control the width, delay and shape of the window (note that the window is not treated as a probability distribution; the gamma window just provides a convenient parametric form). We found previously that the best-fitting delay is highly correlated with the width and is close to the minimum possible value for a given width and shape, and therefore constrained the delay to take this minimum value to reduce the number of free parameters. For a structure-yoked response, we predict that the window will scale with the temporal scaling of the stimulus structures, which will change the width but not the shape. Therefore, we constrained the shape of the window to be the same for compressed and stretched speech for our main analysis. When we did not constrain the shape to be the same, there was more variance between the estimates for stretched and compressed speech, and overall we observed slightly higher structure yoking (median structure-yoking index of 0.15 for untied shapes versus 0.04 for tied shapes; Extended Data Fig. 4). To determine whether these differences were primarily a result of lower reliability or genuine differences in the integration window shape, we estimated integration windows using two different splits of data (non-overlapping segments). We found that tying the shapes increased the reliability of the estimates, measured as the Spearman correlation between splits (untied correlation, 0.614; tied correlation, 0.716). We also found that tying improved predictions for untied estimates. Specifically, we found that correlating the untied estimates from one split with the tied estimates from another split increased the correlation (0.668) compared with just correlating two untied estimates from different splits (0.614). These results suggest that the primary effect of tying is to enhance the reliability of the estimates. We also measured the average cross-context correlation within each annular ROI as a simple, model-free way of assessing the average integration window for stretched and compressed speech (Fig. 3e). The results of this model-free analysis support our model-dependent results by showing that integration windows increase in non-primary regions but change little with structure duration. Although these results do not rule out the possibility that there might be small changes in the integration window between stretched and compressed speech, they point to the same conclusion: that time-yoked integration predominates throughout the auditory cortex, even in non-primary regions.

We also applied our parametric window analysis to our STRF and phoneme integration models, which revealed the expected result with time-yoked integration for our STRF model and structure-yoked integration for our phoneme model (Extended Data Fig. 2). It was not possible to apply our parametric window analysis to our DANN model because, unlike the brain, the model’s responses are acausal, and a gamma-distributed window is therefore inappropriate.

Timecourse rescaling

We measured the degree to which the neural response timecourse rescales by correlating the response timecourse to stretched and compressed speech after rescaling the timecourse for the compressed speech (accomplished by resampling; that is, upsampling by a factor of three, and then discarding the additional time points). After rescaling, we measured the average correlation between all pairs of stimulus repetitions, and we measured a ceiling for this across-condition correlation (stretched versus rescaled compressed) by measuring the average test–retest correlation across repetitions for the same condition (stretched versus stretched and rescaled compressed versus rescaled compressed). We only used responses to the intact 27 s stimuli for this analysis. We applied the same analysis to the trained and untrained DANN models.

Statistics and reproducibility

No formal tests were used to determine the sample size, but the number of participants (15 in our primary experiment) was larger than in most intracranial studies, which often test fewer than ten participants. The only data inclusion or exclusion criterion was the presence of a reliable response to sound (described above). Unlike most intracranial studies, we obtained responses from a sufficient number of participants to perform across-subject statistics and used a LME model to account for subjects as a random effect (using lmefit.m in MATLAB). To evaluate whether the integration windows differed between stretched and compressed speech, we used the following model:

$$\Delta i \sim 1+\left(1|}\right)$$

which models the logarithmic difference in integration windows between stretched (is) and compressed (ic) speech using a fixed effects intercept plus a subject-specific random intercept.

To examine the effect of distance (d) on overall integration windows (\(\bar\)), we modeled the overall integration window across stretched and compressed speech as a function of distance to the primary auditory cortex:

$$\bar \sim 1+d+\left(1+d|}\right).$$

To evaluate the effect of distance on structure yoking, we modeled the difference in integration windows (which is proportional to the structure-yoking index) as a function of distance:

$$\Delta i \sim 1+d+\left(1+d|}\right).$$

To evaluate whether the structure yoking increased at longer timescales, we modeled the difference in integration windows as a function of the overall integration window:

$$\Delta i \sim 1+\bar+\left(1+\bar|}\right).$$

We evaluated whether there was a significant effect of timecourse rescaling by fitting a model on the difference in correlations between rescaled and non-rescaled responses:

$$_}}-_} \;}} \sim 1+\left(1|}\right).$$

We used a Bayesian implementation of LME models because we found that maximum likelihood estimates often resulted in zero-variance random effects terms for subjects. Bayesian LME models were implemented in STAN (CmdStan v.2.35). For a model with both intercepts and slopes, the Bayesian model took the form:

$$y \sim N(_+_x+_}+_}x,\sigma )$$

where β0 and β1 are the fixed intercepts and slopes, respectively, and S0,s and S1,s are the random intercepts and slopes for subjects. We used weakly informative priors for all parameters:

$$\sigma \sim \rm(0,2.5)$$

We report the standard deviation and 90% credible intervals of the posterior distribution for the fixed effect parameters of interest. Data (y) were standardized to unit variance, and predictors (x) were standardized and demeaned. The standardization factors were accounted for when reporting effect sizes (for example, multiplying by the standard deviation of the data when reporting intercept terms). We verified convergence by checking \(\hat\), which was always close to 1 (±0.001), and by examining trace plots, which showed clear evidence of mixing (ten chains, 10,000 samples per chain, 1,000 sample burn-in period). Some chains occasionally had divergent transitions, which we addressed using a non-centered parametrization and a high adapt_delta parameter (0.99). Parameters were initialized to the following values: β0, β1: 0, τ0,s, τ1,s: 0.1, σ: 1. Results were robust to all of these choices (that is, similar for centered parametrization, robust to initialization). Posterior distributions showed clear evidence of unimodality.

Bootstrapping was used to compute error bars (Fig. 3c,d), resampling both subjects and electrodes with replacement (resampling subjects, and for each sampled subject, resampling electrodes). Error bars plot the central 68% interval (equivalent to one standard deviation) of the bootstrapped distribution.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.