We modeled the high gamma (70-150 Hz) amplitude recorded on 331 speech-responsive electrodes located over the left superior temporal gyrus (STG) in 11 participants while they passively listened to 438 naturally spoken sentences from the Texas Instruments and Massachusetts Institute of Technology (TIMIT) acoustic-phonetic corpus (Garofolo et al., 1993) . High gamma amplitudes in neural voltage recordings are known to correlate with the firing rates (Dubey and Ray, 2020; Manning et al., 2009; Ray et al., 2008; Ray and Maunsell, 2011; Scheffer-Teixeira et al., 2013) and dendritic processes ( Bédard et al., 2006 ; Leszczyński et al., 2020 ; Miller et al., 2009 ; Suzuki and Larkum, 2017) of neurons near the electrode ( Buzsáki et al., 2012 ), and we use them here as a proxy for the level of population activity under the ECoG electrodes. Using our model, we show that high gamma responses to speech stimuli across hundreds of electrodes can be parsimoniously represented as a combination of a few lowdimensional latent state responses to specific feature events in the stimulus. Two latent states in particular, corresponding to the sentence onset and peak rate features, reflect a large proportion of the explained variance in the model, and their dynamic properties suggest specific computational roles in the speech perception network.

Successful previous models of high gamma activity over STG have taken two different approaches: using supervised regression to model single-electrode responses as a function of spectral or linguistic characteristics in the audio speech signal (Aertsen and Johannesma, 1981; Holdgraf et al., 2017; Mesgarani et al., 2014; Oganian and Chang, 2019; Theunissen et al., 2001) , and using unsupervised dimensionality reduction to infer latent states without reference to the characteristics of the audio stimulus (Hamilton et al., 2018) .

The advantage of the single-electrode regression models is that they characterize the relationship between the neural responses and acoustic features in the speech signal. In the models, the high gamma responses on individual electrodes are considered to be the result of a convolution of time-dependent receptive fields with corresponding time series of acoustic features. The classic spectrotemporal receptive field (STRF) model, for example, uses a mel spectrogram of the stimulus as the acoustic feature representation, resulting in a framework where the neural receptive fields act as a linear filter on the speech spectrogram (Theunissen et al., 2001). Based on the observation that electrode activity over STG reflects information at the level of phonetic features rather than individual phonemes (Mesgarani et al., 2014) , Oganian and Chang (Oganian and Chang, 2019) used an event-based feature representation to capture these effects and to show that some electrodes additionally have responses triggered by sentence onsets and sharp transients in the acoustic envelope of the speech signal, called peak rate events. While these models have been instrumental in describing the response patterns on individual electrodes, they fail to capture latent dynamics that are shared across multiple electrodes, which could uncover computational principles at work at a larger spatial scale.

An alternative approach uses unsupervised dimensionality reduction to investigate latent structure in neural responses to speech (e.g. (Hamilton et al., 2018) ). Using convex nonnegative matrix factorization, they showed that electrodes can be naturally classified into two groups, <onset= electrodes that have a short increase in high gamma activity at the onset of a sentence, and <sustained= electrodes that show increased high gamma activity throughout the stimulus. This observation is also apparent using principal component analysis, in which the first component has characteristic sustained profile, and the second component has the onset profile (See Supplementary Figure S1). Note that the high gamma signals are not intrinsically low-dimensional: 2 dimensions capture only 24% of the variance in speech responsive electrodes (comparable to 16.9% of the variance in all electrodes captured in the first two clusters of (Hamilton et al., 2018) ) and 189 dimensions are necessary to capture 80% of the variance. This could be related to the high-dimensional nature of the task: in an unsupervised framework in which the system responds to stimulus features, the response dimensionality needs to be at least as high-dimensional as the task itself (Gao et al., 2017; Stringer et al., 2019) . Furthermore, both of these components are timelocked to sentence onset, and it is difficult to connect them or higher components to other speech features, possibly because the dynamics related to other features are not orthogonal to the sentenceonset subspace or to each other. In particular, the dependence of the neural responses on the peak rate events is not apparent from this analysis, and a model that could capture latent dynamics related to peak rate would be valuable for describing population encoding of shorter timescales.

We chose to use a model that combines the advantages of the regression and dimensionality reduction approaches, using a multivariate integrative reduced rank regression model (iRRR) (Li et al., 2019) to estimate the latent dynamics attributed to each speech feature separately. This group-reduced-rank model partitions the expected neural activity into a separate latent state for each feature, choosing the best latent dimensionality for each feature while penalizing the total dimensionality across all features. The model uses a multivariate adaptation of the event-based regression framework of Oganian and Chang (Oganian and Chang, 2019) . In matrix form, the model has the following structure: Ā = ∑ÿ=1 ÿÿ ÿ + Ā (1) Where Ā is the Ā × þ matrix of high gamma amplitude values across electrodes and timepoints, each ÿÿ (Ā × ÿ ) represents the delayed feature events for feature Ā , and Ā (Ā × þ ) is Gaussian noise, assumed to be uncorrelated across electrodes (Ā : number of timepoints; þ : number of electrodes, ÿ : number of delays, ā : number of features). In Ā , ÿ , and Ā , the timepoints corresponding to subsequent sentence stimuli are stacked together. The coefficient matrices ÿ (D ×  N) are the multivariate temporal response functions (MTRFs), representing the responses of each electrode to the given feature across electrodes and delays (up to 750ms).

Only speech responsive electrodes over STG were used for this analysis, defined using single-electrode fits to a linear spectrotemporal model (see electrode selection in Methods). Figure 1A shows the electrodes that were used, colored by the testing r2 value of the fitted spectrotemporal model: STG electrodes with r2>0.05 were used for subsequent analyses (þ = 331).

We fit the regression model using integrative reduced-rank regression (iRRR) (Li et al., 2019) , which applies a penalty based on a weighted sum of the nuclear norms of the feature matrices (see Methods for more detail):

þ Ā ∈ý × 21ÿ ||Ā 2 ∑ÿ=1 ÿÿ ÿ ||ℱ2 + λ ∑ÿ=1 ÿ || ÿ ||∗ { ÿ,ÿýýý}ÿ=1 = argmin (2) where || ⋅ ||ℱ represents the Frobenius (L2) norm, the ÿ s are chosen to balance the regularization across features, and λ is a regularization parameter. The notation || ⋅ ||∗ represents the nuclear norm, or the sum of the singular values of the bracketed matrix. The nuclear norm penalty acts as an L1 penalty on the singular values of each feature matrix, so the regression tends to find solutions where the feature matrices are low-rank (i.e. sparse in the singular values). Because many of the singular values will be zero, the fitted feature matrices can be represented using a low-dimensional singular value decomposition:

ÿ = ā ÿ ÿÿ Ăÿÿ (3) where ā ÿ is ÿ × Ā , ÿÿ is Ā × Ā , and Ăÿÿ is Ā × þ , for some Ā < þ . In other words, the full multivariate feature receptive fields can be represented with a small number of patterns across time (columns of ā ÿ ), patterns across electrodes (rows of Ăÿÿ ), and corresponding weights (values on the diagonal of ÿÿ ). The number of dimensions Ā can be different for each feature, and it comes from balancing the contribution 196 197 198 199 200 201 202 203 204 205 206 207 208 of the feature to the first term (the mean squared error) with the contribution of the feature to the second term (the nuclear norm penalty), relative to other features. Increasing the tuning parameter λ will tend to increase the total number of dimensions used across all features.

For comparison, we also fit the same model using ordinary least squares (OLS) and ridge regression where a separate regularization parameter was chosen for each electrode. All models were fit using 10fold cross validation. For the iRRR and ridge models, the regularization parameters were fit with an additional level of 5-fold cross-validation nested within the outer cross-validation. regression (Ridge). 95% confidence intervals were estimated using the standard error of the mean across cross-validation folds (see Methods). Significance was assessed for comparisons using two-sided paired ttests across cross-validation folds, *** p<0.0005. C: Total explained variance, computed as the testing r2 computed over all speech-responsive electrodes. D: Group nuclear norm, meaning the penalty term from the iRRR model (see Equation 2). E: The effective number of parameters for the fitted models. F: Unique explained variance for each feature (over all speech-responsive electrodes), expressed as a percentage of the variance captured by the full model. Comparing individual features, both timing features have significantly more unique explained variance than all phonetic features, after Bonferroni correction over pairs (left). Also shown is the unique explained variance for the combined timing features (sentence onset and peak rate) and the combined phonetic features (right). When the features are grouped, the phonetic features capture more unique explained variance than the timing features. iRRR outperforms models that treat each electrode individually, and sentence onset and peak rate capture more of the variance than phonetic features Figure 1C-E compare the three different fitting frameworks: OLS, ridge regression, and iRRR. Because the regression framework is the same for all three, the fitted models have very similar total explained variance (r2 computed over all electrodes, Figure 1C), but iRRR by design achieves a much smaller nuclear norm (Figure 1D), which results in solutions that can be described with 94% fewer parameters than OLS and ridge regression (Figure 1E). The fact that the iRRR model captures as much information as the singleelectrode models using far fewer parameters suggests that substantial feature-related information is shared across electrodes.

In order to compare the contribution to the model of the different features, we fit reduced versions of the iRRR model with each feature left out. From there, we could compute the percent of explained variance by comparing the r2 of the full model (ÿ2Ă ) to the r2 of the model without feature Ā (ÿ22ÿ): 100 × (ÿ2Ă 2 ÿ22ÿ)/ÿ2Ă (4) Figure 1F shows the result of this analysis: sentence onset and peak rate explain a larger percentage of the full model variance than each of the phonetic features (p<0.0005 for all comparisons using a twosided paired t-test after Bonferroni correction). This suggests that these two timing features reflect a substantial amount of the speech-induced response across STG.

When the features are grouped into timing (sentence onset and peak rate) and phonetic (all other features) groups, both groups explain a large proportion of the variance (15% and 22%, respectively). Comparing the groups, however, the phonetic features explain more of the unique variance than the timing features (p<0.0005, two-sided paired t-test). This could be surprising in light of the individual feature comparisons: while timing features capture more explained variance than phonetic features when compared individually, when combined they capture less explained variance. This is likely due to (1) correlations between individual phonetic features that lead to lower individual unique explained variance and (2) the fact that more electrodes respond to sentence onset and peak rate than individual phonetic features (Oganian and Chang, 2019) , meaning that sentence onset and peak rate have more widespread spatial support than the more spatially localized phonetic features. This more widespread 249 250 251 spatial support means that the iRRR model is better able to consolidate the activity patterns across multiple electrodes, i.e. capture the latent dynamics, for the sentence onset and peak rate features than for the phonetic features. Accordingly, the following two sections describe the latent state representations for the sentence onset and peak rate features in more detail.

The model fit captures known response differences between pSTG and mSTG In Hamilton and colleagues9 (Hamilton et al., 2018) unsupervised model, the <onset= cluster of electrodes was found to occur primarily over the posterior portion of STG (pSTG). This observation led them to propose that pSTG may play a role in detecting temporal landmarks at the sentence and phrase level, because the short-latency, short-duration responses to sentence onsets in pSTG would be able to encode the event time with high temporal resolution. This idea fits well within a long history of evidence that stimulus responses in mSTG have longer latencies and longer durations than those in pSTG (Hamilton et al., 2020; Jasmin et al., 2019; Yi et al., 2019) . Here, the model fits recapitulate these known differences between mSTG and pSTG.

As discussed above (Equation 3), the feature response matrices that are fitted by the iRRR model can be decomposed into a small number of components across time (<time components=, columns of ā ÿ ), components across electrodes (<spatial components=, rows of Ăÿÿ ), and corresponding weights (values on the diagonal of ÿÿ ). Figure 2 shows the Sentence Onset and Peak Rate fitted feature matrices decomposed in this way (Since ā ÿ and Ăÿ are orthonormal, their columns are unit vectors: as a result, their units are arbitrary and can be best interpreted in relative terms).

Figures 2A and B show the time components scaled by their corresponding weights, and Figures 2C and D show the first two spatial components. To illustrate how the low dimensional components map back to the response functions for individual neurons, Figures 2E and F show the individual electrode response functions (rows of f), colored by the spatial component from Figures 2C and D.

Looking at the left panel of Figures 2C and 2E, we can see that electrodes that have large values in the first spatial component (red circles in Figure 2C, left) have relatively larger overall responses to sentence onset events (red lines in Figure 2E, left). These electrodes occur primarily over pSTG, which is in line with previous findings (Hamilton et al., 2018) .

For peak rate, the first component plays the same role: electrodes that have larger values in the first spatial component (Figure 2D, left) have relatively larger overall responses to peak rate events (Figure 2F, left). Electrodes with large peak rate responses are not limited to pSTG like sentence onset electrodes: rather, they are distributed over all of STG. In other words, the encoding of peak rate in STG is not focal but is distributed over centimeters of cortex, suggesting a representation on a large spatial scale. Interestingly, the second component does appear to have a spatial distinction between pSTG and mSTG: electrodes with positive values for the second component tend to occur over pSTG, while electrodes with negative values for the second component tend to occur over mSTG (Figure 2D, right). The negative and positive values distinguish response functions by their temporal response profile: positive values correspond to electrodes that have an early peak rate response, while negative values correspond to electrodes that have a late peak rate response (Figure 2F, right). This suggests that peak rate responses over pSTG are faster than peak rate responses over mSTG. 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328

Feature latent states have rotational dynamics that capture continuous relative timing information To show how the latent states behave during the presentation of a stimulus, we used the fitted model to predict the dynamics in each latent state during the presentation of the sentence <They9ve never met, you know= (Figure 3). Predictions from the model can be computed in latent space using the decomposition defined in Equation 3:

(5) Āÿ; āþÿā = ÿÿā ÿ ÿÿ The sentence onset latent space has 5 dimensions and the peak rate latent space has 6 dimensions. While the sentence onset feature only occurs once at the beginning of the stimulus, evoking a single response across the sentence onset dimensions, the peak rate feature occurs several times, and the dynamics of the peak rate latent state do not go back to baseline in between peak rate events (Figure 3B and C). Plotting the top three dimensions, which capture more than 98% of the variance in the coefficient matrices ( ÿ), shows cyclical dynamics for both sentence onset and peak rate (Figure 3D and

E): the sentence onset state rotates once at the beginning of the sentence, and the peak rate latent state rotates 3-4 times, once after each peak rate event.

To quantify this effect, we used jPCA (Churchland et al., 2012) to identify the most rotational 2 dimensional subspace within the top three components of ÿ. These planes capture 31.8% and 20.3% of the variance in the sentence onset and peak rate coefficient matrices, respectively, and they highlight the cyclical dynamics that were visible in the top 3 dimensions (Figure 3F and G).

Note that seeing cyclical dynamics in the latent states is not necessarily surprising: the coefficient matrices ÿ describe smooth multivariate evoked responses that will tend to start and end at the same baseline. We highlight them here to motivate a geometrical argument for the computational role of the peak rate responses (see Discussion) and to make the case that the structure of the peak rate responses enables them to act as a temporal context signal against which other features are organized. In order for the peak rate latent state to play this role, the trajectories should be sufficiently spread out in latent space to enable downstream areas to decode the time relative to the most recent peak rate event using just the instantaneous latent state. We investigate whether this is true in the next section.

Latent states from the model can be used to decode time relative to feature events So far, we have described how the model is fit using known feature event times, and how the fitted model can be used to predict responses given new feature events. We also wanted to know whether the model fit could be used to decode the timing of events, which would indicate that sufficient information is contained in the feature responses for downstream areas to use them as temporal context signals.

The set of spatial components for each feature defines a feature-specific subspace of the overall electrode space. The projection of the observed high gamma time series onto this subspace is an approximation of the feature latent state (note that it is not exact, because the different feature subspaces are not orthogonal to each other): Ā̃ÿ,āÿĀĀ = ĀĂÿ (6) We asked whether this latent projection time series could be used to decode the time since the most recent feature event.

Figure 4 shows the result of this analysis: a perceptron model was trained to decode the time since the most recent feature event up to 750 ms, given either the activity on the full set of electrodes or the projection of the electrode activity onto the corresponding feature subspace (see Methods). The decoder for sentence onset performs slightly better when using all electrodes, which may be due to the large proportion of the overall activity that is time-locked to sentence onsets (see Supplementary Figure S1). For all other features, however, decoder performance using the reduced-dimensional latent subspaces performs even better than decoding using the full dimensional activity across electrodes (paired t-test over 10 cross validation folds, p<0.05 with Bonferroni correction across features). Because no information is gained in the projection operation, this is an indication that projecting onto the latent subspaces increases the signal to noise ratio, i.e. removes activity that is irrelevant to decoding relative time.

Participants included 11 patients (6M/5F, age 31 +/- 12 years) undergoing treatment for intractable epilepsy. As a part of their clinical evaluation for epilepsy surgery, high-density intracranial electrode grids (AdTech 256 channels, 4mm center-to-center spacing and 1.17mm diameter) were implanted subdurally over the left peri-Sylvian cortex. All procedures were approved by the University of California, San Francisco Institutional Review Board, and all patients provided informed written consent to participate. Data used in this study was previously reported in (Hamilton et al., 2018) .

Stimuli consisted of 499 English sentences from the TIMIT acoustic-phonetic corpus (Garofolo et al., 1993) , spoken by male and female speakers with a variety of North American accents. Stimuli were presented through free-field Logitech speakers at comfortable ambient loudness (~70 dB), controlled by a custom MATLAB script. Participants passively listened to the sentences in 4 blocks, each lasting about 4 minutes. A subset of 438 sentences were selected for analysis that were heard once by all 11 subjects. The sentences had durations between 0.9 and 2.6s, with a 400ms intertrial interval. Neural recordings and electrode localization Neural recordings were acquired at a sampling rate of 3051.8 Hz using a 256-channel PZ2 amplifier or 512-channel PZ5 amplifier connected to an RZ2 digital acquisition system (Tucker-Davis Technologies, Alachua, FL, USA).

Electrodes were localized by coregistering a preoperative T1 MRI scan of the individual subject9s brain with a postoperative CT scan of the electrodes in place. Freesurfer was used to create a 3d model of the individual subjects9 pial surfaces, run automatic parcellation to getindividual anatomical labels, and warp the individual subject surfaces into the cvs_avg35_inMNI152 average template (Desikan et al., 2006; Fischl et al., 2004) . More detailed procedures are described in (Hamilton et al., 2017) .

For each electrode, the high gamma amplitude time series were extracted from the broadband neural recordings as follows (Hamilton et al., 2018; Oganian and Chang, 2019) . First, the signals were downsampled to 400 Hz, rereferenced to the common average in blocks of 16 channels (blocks shared the same connector to the preamplifier), and notch filtered at 60, 120, and 180 Hz to remove line noise and its harmonics. These LFP signals were then filtered using a bank of 8 Gaussian filters with center frequencies logarithmically spaced between 70 and 150 Hz. Using the Hilbert transform, the amplitude 21 of the analytic signal was computed for each of these frequency bands, and for each electrode the high gamma amplitude was defined as the first principal component across these 8 frequency bands. Finally, the high gamma amplitude was further downsampled to 100Hz and z-scored based on the mean and standard deviation across each experimental block.

In order select speech-responsive electrodes over STG, electrodes were included (1) if they were located over the STG, as identified in the Freesurfer anatomical parcellation of the individual subject cortical surface, and (2) if their high gamma activity was well-predicted by a linear spectrotemporal model (Hamilton et al., 2018) .

For this single electrode analysis, the model had the form of a spectrotemporal receptive field (STRF): (ā) = ∑ÿ ∑τ ÿ(Ā, ā 2 τ)´(τ, Ā) + ÿ(ā) (7) where is the high gamma amplitude on a single electrode, ÿ is the mel spectrogram of the speech audio signal over frequencies between 75Hz and 8kHz, coefficients ´ vary across frequencies and delays between 0 and 500ms, and ÿ is the zero-mean Gaussian error term. Ridge regression was used to fit the models (see Model fitting below for details of the ridge regression framework): the data were split into 80% training and 20% testing data sets, the training data was used to choose the alpha parameter according to a 5-fold cross-validation, the full training data was fit using the chosen ³ parameter, and the r2 was assessed on the testing data (see Explained Variance Calculation below for computation of r2). Electrodes with r2>0.05 were included in subsequent analyses. The selected electrodes and their corresponding r2 values are shown in Figure 1A.

Ā = ∑ÿ=1 ÿÿ ÿ + Ā The multivariate temporal receptive field model has the following structure:

(8) Where: ● Ā is the Ā × þ matrix of z-scored high gamma amplitude values across electrodes and timepoints.

The time dimension represents a concatenation of all 438 sentence stimuli that were heard by every subject, from 500 ms before sentence onset until 500 ms after sentence offset (132,402 timepoints, later split for cross validation, see Model Fitting below). The electrode dimension includes speech-responsive electrodes from all subjects (331 electrodes). ● Each ÿÿ (Ā × ÿ ) represents the delayed feature events for feature Ā. The first column contains the feature events across time (1 representing an event occuring, 0 otherwise. For peak rate, events were coded by a real-valued magnitude, see Figure 1B). Following columns contain the same time series, offset by time-delays between 10 ms and 750 ms (76 delays). There were 12 features: sentence onset, peak rate, dorsal, coronal, labial, high, front, low, back, plosive, fricative, and nasal (described below). ● Ā (Ā × þ ) is Gaussian noise, assumed to be uncorrelated across electrodes ● ÿ (ÿ × þ ) are the coefficient matrices, i.e. the multivariate temporal response functions (MTRFs), representing the responses of each electrode to the given feature across electrodes and delays ● Ā : number of timepoints; þ : number of electrodes, ÿ : number of delays, ā : number of features. Sentence onset was defined as the sound onset time for the sentence stimulus. Peak rate was extracted by taking the derivative of the analytic envelope of the speech signal: the peak rate event times were the times when the derivative reached a maximum, and the peak rate magnitude was the value of the derivative at that time point (Oganian and Chang, 2019) . Phonetic feature event times (dorsal, coronal, labial, high, front, low, back, plosive, fricative, nasal) were extracted from time-aligned phonetic transcriptions of the TIMIT corpus, which were timed to the onset of the respective phonemes in the speech signal (Garofolo et al., 1993) .

The model was fit using ordinary least squares (OLS), ridge regression, and iRRR. The difference between the three is the objective function that is minimized to choose the fitted coefficient matrices: 1 { ÿ,ÿ þ }ÿ=1 = þaĀr∈gým×in 2ÿ ‖Y 2 ∑ÿ=1 ÿÿ ÿ‖ℱ2

þ Ā ∈ý × 21ÿ||Ā 2 ∑ÿ=1 ÿÿ ÿ ||ℱ2 + ³ ∑ÿ=1 || ÿ ||ℱ2 { ÿ,ÿÿýĀþ}ÿ=1 = argmin

þ Ā ∈ý × 21ÿ ||Ā 2 ∑ÿ=1 ÿÿ ÿ ||ℱ2 + λ ∑ÿ=1 ÿ || ÿ ||∗ { ÿ,ÿýýý}ÿ=1 = argmin (9a) (9b) (9c) The weights used for the iRRR model were chosen to balance the different features (Li et al 2019) : ÿ = σ(ÿÿ, 1) {√þ + √ÿ(ÿÿ)} /Ā

(10) where σ(ÿÿ, 1) is the first singular value of the matrix ÿÿ and ÿ(ÿÿ) = ÿ is the rank of matrix ÿÿ. All predictors ÿÿ and responses Ā were column-centered before fitting the models.

In order to compute confidence intervals for model performance metrics, models were fit using 10-fold cross validation, using group cross validation to keep time points corresponding to the same sentence stimulus in the same fold. For ridge regression and iRRR, an additional nested 5-fold cross validation was used to choose the ³ and λ parameters within each fold of the outer cross-validation. 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828

Note that the approach of using a regression framework to fit a group-reduced rank model of neural activity has been used before (Aoi et al., 2020; Aoi and Pillow, 2019) : the iRRR framework differs in that it uses an L1 relaxation, resulting in a convex optimization formulation that can be fit efficiently using alternating direction method of multipliers.

ÿ2 = 1 2 þþþþāÿÿāĀ where the ÿÿÿþĀ is the residual sum of squares computed on the testing dataset: ÿÿÿþĀ = ||Ā 2 ∑ÿ=1 ÿÿ ÿ ||ℱ2 (12) and ÿÿāĀā is the total sum of squares computed on the testing dataset:

(13) ∑ÿ=1 ÿ || ÿ ||∗ The group nuclear norm (Figure 1D) was computed as the penalty term in the iRRR model: (14) Because OLS and ridge regression yield full-rank coefficient matrices, the number of parameters (Figure 1E) used for both is ÿþ . For iRRR, the number of parameters is Ā (ÿ + þ + 1 ), based on the singular value decomposition described in Equation 3, reproduced here:

ÿ = ā ÿ ÿÿ Ăÿÿ (15) Unique explained variance for each feature (Figure 1F) was computed by fitting a reduced iRRR model without the feature Ā , and then comparing the total explained variance of the full model ÿ2Ă to the total explained variance of the reduced model ÿ22ÿ. The reduced iRRR model was fit using the same λ value as the full model, chosen using nested cross validation on the full model as described above. For the <all timing= category, the reduced model was fit without sentence onset and peak rate, and for the <all phonetic= category, the reduced model was fit without the phonetic features. The unique explained variance was expressed as a percentage of the full model: 100 × ÿ2Ăÿ2Ă2ÿ−2Ā (16) All metrics are reported in terms of the mean across the 10 folds of the cross validation, and 95% confidence intervals are ±ā9,0.975Ā/√10, where Ā is the sample standard deviation across the 10 cross validation folds. Note that these confidence intervals do not account for the dependence between cross-validation folds due to reuse of samples in training and testing sets, and may therefore be smaller than the true intervals (Austern and Zhou, 2020; Bates et al., 2021; Bengio and Grandvalet, 2004) .

Significant differences between conditions were assessed using paired two-tailed t-tests across crossvalidation folds (Dietterich, 1998) for the following comparisons (with the resulting p-value ranges): 1. Total explained variance for OLS vs Ridge (p>0.05), OLS vs iRRR (p<0.0005), and Ridge vs iRRR (p<0.0005). 2. Unique explained variance of sentence onset vs each acoustic-phonetic feature and peak rate vs each acoustic-phonetic feature. Here the p-values were Bonferroni corrected across the (2 timing features times 10 acoustic-phonetic features) 20 comparisons. After correction, all comparisons were significant with p<0.0005. 3. Unique explained variance of the combined timing features vs the combined acoustic-phonetic features (p<0.0005).

Similar to the confidence intervals described above, the significance tests did not account for the dependence between cross-validation folds and may therefore have an inflated type II error (Austern and Zhou, 2020; Bates et al., 2021; Bengio and Grandvalet, 2004) .

Given a fitted model, the predicted latent response to a stimulus matrix ÿÿ is (reproduced from Equation 5): Āÿ; āþÿā = ÿÿā ÿ ÿÿ (17) where, as before, ÿÿ (Ā × ÿ ) represents the delayed feature events for feature Ā , ā ÿ is the ÿ × Ā time components for feature Ā , and ÿÿ is a diagonal matrix containing the weights for each component (Ā × Ā ). Āÿ; āþÿā is a Ā × Ā matrix representing the predicted response within the Ā -dimensional latent space of the feature. Figure 3 shows the predicted sentence onset and peak rate responses to the sentence <They9ve never met, you know=. jPCA The plane of fastest rotation for the sentence onset and peak rate latent states (Figure 3C) was identified by applying jPCA (Churchland et al., 2012) to the feature coefficient matrices ÿ. Using jPCA, we modeled the temporal receptive fields in the coefficient matrix as a linear dynamical system evolving over delays: ýþ Ā (ā) = ý ÿ(ā) (18)

ýā where ā indexes the delay dimension of ÿ, so the dynamical system describes the evolution of an þ dimensional dynamical system over ÿ timepoints. By approximating the derivative on the left hand side using first differences, the transition matrix ý can be fit using regression. Furthermore, the purely rotational component of the transition matrix can be isolated by constraining the matrix ý to be skewsymmetric, having purely imaginary eigenvalues that come in complex conjugate pairs. The pair of eigenvectors with the largest magnitude eigenvalues describes the plane with the fastest rotations.

It is important to note that jPCA identifies planes with fast rotational dynamics, regardless of whether they capture a large proportion of the variance of the dynamics in the original dynamical system. Classic jPCA uses PCA in preprocessing in order to confine the analysis to six dimensions of largest variance. Here, the iRRR model chooses Ā dimensions for each feature that are most valuable to the overall fit of the model. Hence there was no need to perform additional PCA to reduce the dimensionality. However, because the coefficient matrices had dimensions capturing very little variance, we did subselect components to capture 98% of the variance of the coefficient matrices. For both sentence onset and peak rate, this corresponded to the top 3 components. Hence the jPCA plane represents the plane of maximal rotation within a 3-dimensional subspace capturing 98% of the variance in the 5-dimensional (or 6-dimensional) coefficient matrix for sentence onset (or peak rate). If we had used more components for the jPCA computation, the rotational dynamics would be stronger but they would capture much less of the variance (using Ā dimensions vs using 3 dimensions: 2.8% vs 31.8% for sentence onset and 4.8% vs 20.3% peak rate), making them less informative about the overall population dynamics.

Once the jPCs were computed using the coefficient matrices, the predicted trajectory for a given stimulus (Figure 3F and G) is calculated as:

(19) Āÿ;ĀĀÿý = ÿÿ ÿ ÿ = [Ā1 + Ā2, ÿ(Ā1 2 Ā2)] where Ā1 and Ā2 are the eigenvectors with largest eigenvalues of the skew-symmetric matrix ý defined above. ÿ is therefore the þ × 2 projection matrix from electrode space onto the plane of highest rotation from jPCA.

For the decoding analysis (Figure 4), a perceptron model was trained to predict the time relative to the most recent feature event (up to 750 ms). The model was designed using the MLPRegressor class of the sklearn package, with one hidden layer with 20 hidden units using a logistic activation function. We used a simple perceptron model in order to account for possible nonlinearities in the mapping from electrode space / feature latent space to relative times.

Using the same cross-validation framework that was used for iRRR model fitting, the perceptron model was trained using the training data (high gamma amplitudes) either across all electrodes Ā or using the projected data onto the latent state subspace (reproduced from Equation 6): Ā̃ÿ,āÿĀĀ = ĀĂÿ (20) where Ăÿ is the þ × Ā matrix of electrode components for feature Ā , as above. The Ā × Ā matrix Ā̃ÿ,āÿĀĀ is an approximation of the latent state across time, but it may be contaminated by activity from 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 other features because the Ăÿ matrices do not describe orthogonal subspaces. It also contains activity from noise.

Performance of the models was assessed using r2 (Equation 11) on the held-out testing data for the cross-validation fold. The 95% confidence intervals were computed using the t distribution as described above, and the performance of the models trained on all-electrodes was compared to the performance of the models trained on the latent projections using a two-sided paired t test, as described above (Model Performance Metrics), Bonferroni corrected across the 12 features. The sentence onset model performed better using all electrodes than the latent projection (corrected p<0.05), while the models for all other features performed better using the latent projection than using all electrodes (corrected p<0.05).

Custom Python code to perform the iRRR fits is available online (https://github.com/emilyps14/iRRR_python), which is a port of the Matlab implementation by the original authors (https://github.com/reagan0323/iRRR, Li et al 2019) . Python code for the analysis pipeline described above is also available (https://github.com/emilyps14/mtrf_python). We thank Antin and colleagues (Antin et al., 2021) for their implementation of jPCA in the Python programming language (https://github.com/bantin/jPCA), which we used to perform the jPCA. Table S1. Clinical and demographic details for subjects. Hem = hemisphere of implantation. Subject ID SL01 SL02 SL03 SL04 SL05 SL06 SL07 SL08 SL09 SL10 SL11

L L L L L L L L L L L Hem

Age

Sex