The Hierarchical Cortical Organization of Human Speech Processing

MRI data collection and preprocessing.

Structural MRI data and BOLD fMRI responses from each subject were obtained while they listened to ∼2 h and 20 min of natural stories. For five of the subjects, these data were collected during two separate scanning sessions that lasted no more than 2 h each. For two of the subjects (S1, author A.G.H., and S5, author W.A.d.H.) the validation data (two repetitions of a single story) were collected in a third separate session. MRI data were collected on a 3 T Siemens TIM Trio scanner at the University of California, Berkeley, Brain Imaging Center, using a 32-channel Siemens volume coil. Functional scans were collected using a gradient echo EPI sequence with TR = 2.0045 s, echo time = 31 ms, flip angle = 70°, voxel size = 2.24 × 2.24 × 4.1 mm, matrix size = 100 × 100, and field of view = 224 × 224 mm. Thirty-two axial slices were prescribed to cover the entire cortex. A custom-modified bipolar water excitation radiofrequency pulse was used to avoid signals from fat tissue. Anatomical data were collected using a T1-weighted MP-RAGE (Brant-Zawadzki et al., 1992) sequence on the same 3 T scanner.

Each functional run was motion corrected using the fMRIB Linear Image Registration Tool (FLIRT) from FSL 4.2 (Jenkinson and Smith, 2001). All volumes in the run were then averaged to obtain a high-quality template volume. FLIRT was also used to automatically align the template volume for each run to the overall template, which was chosen to be the template for the first functional run for each subject. These automatic alignments were manually checked and adjusted for accuracy. The cross-run transformation matrix was then concatenated to the motion-correction transformation matrices obtained using MCFLIRT, and the concatenated transformation was used to resample the original data directly into the overall template space. Low-frequency voxel response drift was identified using a second-order Savitsky–Golay filter with a 120 s window, and this was subtracted from the signal. After removing this time-varying mean, the response was scaled to have unit variance (i.e., z scored).

Structural imaging and fMRI were combined to generate functional anatomical maps that included localizers for known regions of interests (ROIs). These maps were displayed either as 3D structures or as flatmaps using custom software (pycortex). For this purpose, cortical surface meshes were first generated from the T1-weighted anatomical scans using Freesurfer software (Dale et al., 1999). Five relaxation cuts were made into the surface of each hemisphere, and the surface crossing the corpus callosum was removed. The calcarine sulcus cut was made at the horizontal meridian in V1 using retinotopic mapping data as a guide. Known auditory and motor ROIs were then localized separately in each subject using standard techniques. To determine whether a voxel was responsive to auditory or motor stimuli, repeatability of the voxel response was calculated as an F statistic given by the ratio of the total variance in the responses over the residual variance. The residual variance was obtained by comparing responses in individual trials to the mean responses obtained from multiple repeats of the stimulus played back or of the motor action. AC localizer data were collected in one 10 min scan. The subject listened to 10 repeats of a 1 min auditory stimulus, which consisted of 20 s segments of music, speech, and natural sound. Motor localizer data were collected during one 10 min scan. The subject was cued to perform six different motor tasks in a random order in 20 s blocks (10 blocks per motor task for a total to 60 blocks). For the hand, mouth, foot, speech, and rest blocks, the stimulus was simply a word located at the center of the screen (e.g., “Hand”). For the Hand cue, the subject was instructed to make small finger-drumming movements with both hands for as long as the cue remained on the screen. Similarly, for the “Foot” cue the subject was instructed to make small toe movements for the duration of the cue. For the “Mouth” cue, the subject was instructed to make small mouth movements approximating the nonsense syllables balabalabala for the duration of the cue—this requires movement of the lips, tongue, and jaw. For the “Speak” cue, the subject was instructed to continuously subvocalize self-generated sentences for the duration of the cue. For the saccade condition, the written cue was replaced with a fixed pattern of 12 saccade targets and the subject was instructed to make frequent saccades between the targets. After preprocessing, a linear model was used to find the change in BOLD response of each voxel in each condition relative to the mean BOLD response. Weight maps for the foot, hand, and mouth responses were used to define primary motor area (M1) and somatosensory area (S1) for the feet (M1F, S1F), hands (M1H, S1H), and mouth (M1M, S1M); supplementary motor areas for the feet and hands; secondary somatosensory area for the feet (S2F) and, in some subjects, the hands (S2H); and, in some subjects, the ventral premotor hand area. The weight map for saccade responses was used to define the frontal eye field, frontal operculum eye movement area, intraparietal sulcus visual areas, and, in some subjects, the supplementary eye field. The weight map for speech production responses was used to define Broca's area and the superior ventral premotor speech area (sPMv).

Linear model fitting.

The relationship between the speech stimulus represented in various feature spaces and the BOLD response was fitted with a linear filter estimated for every single voxel (i.e., VM) (Fig. 2). This filter is equivalent to a finite impulse response (FIR) model (Nishimoto et al., 2011; Huth et al., 2012) and to spatiotemporal receptive fields where the spatial dimensions correspond to the vector dimensions of particular feature spaces. The filters modeled here were estimated at the following four delays: t-2 s, t-4 s, t-6 s, t-8 s. Therefore, the equation for the FIR can be written as follows: Because auditory cortex has shorter hemodynamic delays than does visual cortex (Belin et al., 1999), this model incorporates a 2 s time delay that was not used in earlier vision publications from our group (Nishimoto et al., 2011; Huth et al., 2012). Models were fitted for each feature space independently, for all pairwise combinations of feature spaces, and for the combination of all three feature spaces (for a total of seven linearized models for each voxel).

The linear model weights h⃗_i were estimated with regularized regression methods and cross-validation to avoid overfitting. Before fitting, both the BOLD response and the stimulus values in the feature space were z scored. Because each story had significantly different spectral, articulatory, and semantic content, and because the BOLD signals adapted to these average levels, a different mean and variance was used in the z-scoring operation for each story.

Regularization was performed using ridge regression. A separate ridge regularization parameter was estimated for each voxel, in each subject, and for each model (i.e., each combination of feature spaces). Ridge parameter estimation was performed by repeating a cross-validation regression procedure 50 times, in each subject, for each model. On each cross-validation iteration, 800 time points (consisting of 20 random blocks of 40 consecutive time points each) were selected at random and removed from the training dataset. Then the model parameters were estimated on the remaining 2937 time points, and each of 20 possible regularization hyperparameters were log-spaced between 10 and 10,000. These weights were used to predict responses for the 800 reserved time points, and R² was computed from these data (R² gives the fraction of variance that the predictions explain in the responses). After the 50 cross-validation iterations were complete, a regularization–performance curve was obtained by averaging the sample R² values across the 50 iterations. This curve was used to select the best regularization hyperparameter for the current model in each voxel and in each subject. Finally, the selected hyperparameter and all 3737 training time points were used to estimate the final model parameters.

For voxelwise models that combined multiple feature spaces, two different regularization approaches were investigated. First, the subspaces obtained from performing ridge regression on the individual features were combined, and a single regression analysis was performed in that joint subspace (see Joint ridge regression section below). Second, the features were concatenated without performing any dimensionality reduction and a new optimal ridge parameter was estimated for this joint space. Although results were similar in both cases, the second approach tended to yield a smaller error in the variance partitioning scheme. Only the results obtained with a single ridge parameter are shown in the Results section. However, we also describe the method for the joint ridge regression approach in the Materials and Methods section, because it is a principled approach for comparing nested models, and it can be used to verify the validity of other approaches for combining feature spaces in the context of regularized regression.

All model fitting and analysis was performed using custom software written in Python, which made heavy use of the NumPy (Oliphant, 2006) and SciPy (Jones et al., 2007) libraries.

Signal detection and model validation.

Before voxelwise modeling was performed, the validation dataset was used to determine which voxels were significantly active in response to the stories (story-responsive voxels). The validation set consisted of BOLD responses to two repeats of the same story, as described above (Fig. 2). The correlation between these paired BOLD responses was calculated, and then the significance of this correlation was computed using an exact test. This test gives the probability of finding the observed correlation assuming that the two response vectors are bivariate Gaussians distributed with zero covariance. These p values were then corrected for multiple comparisons across voxels within each subject using the false discovery rate (FDR; Benjamini and Hochberg, 1995). The results of this analysis were used to show which voxels responded to speech but were not used to constrain or bias the voxelwise modeling analyses in any way.

After voxelwise model estimation was performed using the model estimation dataset, we estimated which voxels were significantly predicted by each model. First, the correlation between the actual BOLD response in the validation dataset (averaged across the two repetitions) and the model predictions were calculated, and then the significance of the observed correlation was computed as above. While the correlation between predicted response and actual mean response is an appropriate metric for assessing significance, it is biased downward due to noise in the validation data (Sahani and Linden, 2003; Hsu et al., 2004; David and Gallant, 2005). This is because the actual mean response is calculated using a finite number of repetitions (in this case, 2), and so it contains both signal and residual noise. This noise level is likely to vary across voxels due to vascularization and magnetic field inhomogeneity. The noise in the validation dataset was accounted using the method developed in the study by Hsu et al. (2004). In this method, the raw correlation is divided by the expected maximum possible model correlation (called the noise ceiling) for each voxel.

Joint ridge regression.

The following section describes how to perform a joint ridge regression operation that preserves the individual shrinkages obtained in the single-feature space ridge regressions. This preservation of individual shrinkages in the joint model is important when comparing nested models.

Ridge regression as well as other forms of regularization shrink the parameter space of input features to prevent overfitting. These dimensionality reduction operations rely on feature spaces that have uniform physical dimensions and statistical properties in a Euclidian space. The joint models fitted here combined feature spaces with different units, different numbers of parameters, p, and different degrees of correlation across those parameters. When regularized linear models are fitted separately for each stimulus feature spaces, the solution that best generalizes to a novel dataset will be obtained with a different optimal value for the ridge parameter λ: stimulus feature spaces that have more dimensions and/or are less predictive of the response will require more shrinkage to prevent overfitting. When a model is based on a union of two feature spaces, such as one with a large number of parameters (e.g., semantic) and one with a low number of parameters (e.g., articulation), using a single large shrinkage value to prevent overfitting with this large combined feature space might significantly reduce the contribution to the prediction from the features of the smaller feature space. More importantly, using the same shrinkage in the joint model (i.e., the same projection to the same subspaces) is required to estimate the total variance (union in set theory), and from there both the shared (intersection in set theory) and unique [relative complement (RC) in set theory] variances explained (the variance partitioning in set theoretical terms is further developed below). The total variance explained of two or more regularized (shrunken) feature spaces is the variance explained by the union of these shrunken feature subspaces.

The voxelwise models predict responses, r⃗, from stimulus parameters S. Here r⃗ is a column vector (n × 1) corresponding to the BOLD response of a single voxel as a function of n discrete sampling points in time. S is a (n × p) matrix where each row corresponds to the stimulus at the same n time points as r⃗, and the columns correspond to the values describing the stimulus in its feature space for a number of time slices: the number of columns (p) is equal to the dimension of the feature space (k) times the number of delays used in the FIR model (here p = 4 × k since four time slices are used). The maximum likelihood solution for the multiple linear regression is given by the following normal equation: where h⃗ is the column vector of regression weights (p × 1) also called the model parameters or the filter. The angle brackets (〈 and 〉) stand for averaging the cross-products across time (across rows). Note, more specifically, that the correct unbiased estimate of the stimulus–response cross-covariance (the numerator) and the stimulus auto-covariance (the denominator) are as follows: The prediction can then be obtained by: The inverse of the symmetric and positive definite stimulus auto-covariance matrix can easily be obtained from its eigenvalue decomposition or equivalently from the singular value decomposition (SVD) of S. The SVD of S can be written as follows: where V is a (p × p) matrix of orthonormal input vectors in columns (or left singular vectors), W is a diagonal (p × n) matrix of positive single values, and U is a (n × n) matrix of orthonormal output vectors in columns (or right singular vectors).

The eigenvalue decomposition of S^TS is then given by the following: where W² is the (p × p) diagonal matrix of eigenvalues.

To prevent overfitting (when p is large relative to n), a regularized solution for w⃗ can be obtained by ridge regression. The ridge regression is the maximum a posteriori solution (MAP) with a Gaussian prior on w⃗ with zero mean and covariance matrix given by Iλ.

Under these assumptions, the MAP solution is as follows: If the normal equation can be interpreted as solving for w⃗ in the whitened-stimulus space (uncorrelated by the rotation given by V and normalized W), then the ridge regression decorrelates the stimulus space and provides a weighted normalization where the uncorrelated stimulus parameters with small variance (or small eigenvalues) are shrunken more than those with higher variance (or higher eigenvalues). The level of this relative shrinkage is controlled by the hyper-parameter λ, and its optimal value is found by cross-validation (see Linear model fitting above).

In the following equations, S₁ is the (n × p₁) matrix of features for the first stimulus feature space, and S₂ is the (n × p₂) matrix of features for the second stimulus feature space. S₁₂ is the (n × (p₁ + p₂)) matrix that combines the features from both spaces simply by column concatenation. In the joint ridge approach, the regression is performed in the rotated and scaled basis obtained for each of the models. The stimulus space in that new basis set is noted with a prime in the equations below. The decorrelation step is then performed in the joint stimuli but without performing any additional normalization After whitening the stimuli, a correlation coefficient matrix is created from the covariance matrix to decorrelate the stimuli without further normalization. The stimulus covariance matrix in this new stimulus space (denoted with the prime) can be obtained with S₁₂^′TS₁₂^′ divided by n − 1, or the following: where σ² is variance and c is the covariance between individual parameters in each of the two feature spaces. The first index is for the feature space corresponding to model 1 or 2, and the second index runs over the parameters in that feature space. As one can notice from the form of this covariance matrix, the stimulus parameters are uncorrelated within each subset, but because of the relative shrinkage performed by the ridge they are not perfectly white. Therefore, the variance in the diagonals is not exactly equal to 1 but is slightly smaller and with decreasing values along each block diagonal. If at this stage the weights of the linear regression were to be obtained using a normal equation, the shrinkage performed in the ridge solution would be inverted. To prevent this unwanted normalization, the covariance matrix can be replaced with the correlation matrix. Dividing by the correlation matrix will decorrelate the stimulus features across the two component models while preserving the exact shrinkage from the separate ridge regressions. The correlation matrix obtained from the covariance matrix is given by the following: The combined ridge filter is then calculated as follows: and is used to obtain predictions from the combined model with the following equation: For clarity, the derivation used here involved joint ridge regression on two models, but it can be extended to the joint ridge regression so that any number of feature spaces can be combined.

Partitioning of variance.

To quantify the unique contribution of different stimulus features to the BOLD responses, we estimated the variance explained (R²) uniquely by each individual feature space and the variance explained by the intersections of various combinations of these feature spaces (Lescroart et al., 2015). For this purpose, the results obtained from fitting models from individual feature spaces and combinations of two and three feature spaces were used to estimate R² for all of the nested models. Set theory was then used to calculate the common (as a set intersection) and unique (as a set difference) variances explained. (See Fig. 5 for a graphical representation of this process.) To be succinct, in the remainder of this section, the variance explained by the three feature spaces will be written as sets A–C. First, R² values for the following nested models were directly obtained using the following linear model fitting and cross-validation procedures described above: The shared variances explained by the intersections of two sets was then obtained from the following: Similarly, the variance explained by the intersection of all three sets was obtained from the following: The variance explained by the intersections of two models that did not include the variance explained by the intersection of all three models was then calculated from the following: Finally, the variance solely explained by one model, with no overlap of variance explained by any of the other models, was calculated. This is known as the RC for each pair of models. The relative complement of BC, or BC^RC, is the portion of the variance explained exclusively by model A: Set notation is used here because of its simplicity and its intuitive graphical representation of the results. However, one can easily rewrite these quantities in terms of R² and the sum of errors. For example, if SS₀ is used to represent the total sum of square errors (or the SS of a 0th order model predicting the mean response), then we have the following: and:

Correction of variance partition estimates.

Because empirical estimates of the variance explained by single and joint models contain sampling noise, the set theoretical approach detailed above sometimes produced results that were not theoretically possible. These sampling errors occurred both using the joint regression algorithm or when using a single ridge shrinkage parameter for each model. For example, the estimated variance explained by A ∪ B in the held-out validation dataset was sometimes smaller than the variance explained by A or B alone, due to overfitting of the larger A ∪ B model and sampling error. This happened most often by combining the semantic model, which has a large number of parameters and good predictive power, with the spectral or articulatory model, either of which has a small number of parameters and little additional predictive power in many regions of the brain. This situation produced nonsensical results, such as variance partitions with negative values. To mitigate this problem, a post hoc correction was applied to the estimated variance explained by each model in each voxel. This correction moved the estimates to the nearest values that produced no nonsensical results. Mathematically, this involved estimating a bias term for the variance explained by each model in each voxel. We began by assuming that the estimated variance explained by some model (R²), X̂, is a biased estimate of the true variance explained, X*: X̂ = X* + b_X.

Because there were seven models (each feature space alone, each pair, and all three together), this formulation yields seven bias parameters (b_x). Furthermore, because we know that the size of each variance partition must be at least equal to zero, the set theory equations that give the size of each partition can be used to define seven inequality constraints on the bias terms. Assuming that we want to find the smallest set of bias parameters (in an L2 sense) that produce no nonsensical results, this allowed us to set up a constrained function minimization problem, as follows: where h are our seven inequality constraints.

This procedure was applied separately to the estimated values of the variance explained for each voxel. Applying this correction to simulated data verified that this scheme significantly decreases error, variance, and bias in the estimated variance partition sizes.

Auditory cortex axes and centers of mass.

The center of mass was calculated along two axes of the auditory area for our seven variance partitions. First, all of the voxels within the auditory cortex (as defined by our localizer) were projected onto the following two different axes: an anterior–posterior axis and a medial–lateral axis. The medial–lateral axis was defined as the geodesic distance along the cortical surface from the crown of the superior temporal gyrus (STG), in millimeters. The anterior–posterior axis was defined as the geodesic anterior–posterior distance from the intersection of Heschl's gyrus and STG. The center of mass of each model and axis was then calculated for each subject, as follows: where a_i is the location of the ith voxel projected on the chosen axis (in millimeters), and r is the voxel model performance, expressed as partial correlation, or the positive square root of the partial R².

Bootstrapping was used to calculate the SEs of these calculated centers of mass. For each model, subject, and axis, 1000 points were sampled (with replacement) along the chosen auditory axis, 1000 times. The center of mass was then calculated for each sample, and the SEM was computed from these data.

Linear mixed-effects modeling.

Linear mixed-effects models (lme) were used to compare average responses in the left versus right hemispheres for all the cortical voxels as well as for specific ROIs. In these statistical tests, the subject ID was the random effect. The lme tests were run in R using the lme4 library. For post hoc tests, p values were corrected for multiple comparisons using the FDR (Benjamini and Hochberg, 1995).

Mapping semantic selectivity within auditory cortex.

Variance partitioning shows that some of the response explained by the semantic features in AC cannot be explained away by spectral or articulatory features. However, the semantic model could still be picking up on other features that are correlated with semantics but were not included in the variance partitioning. One possibility is global modulatory effects such as attention or arousal, which one could easily imagine as being correlated with semantic content. If semantic models were capturing global modulatory signals, we would expect to find homogeneous semantic tuning across AC. To test for this possibility, we examined the variability of semantic selectivity within AC. The weights for the semantic model were projected onto the top 10 principal components (PCs) of the semantic models from all seven subjects (for details, see Huth et al., 2016; the PCs used here were taken from the analysis described in that article). This was done to reduce the semantic models to a dimensionality that could be easily examined while preserving as much structure as possible. We then selected the subset of the voxels where the unique contribution of the semantic features was the largest variance partition according to the variance partitioning analysis and computed the variance of the PC projections across that set of voxels for each of the 10 PCs. This analysis was performed for each subject and each hemisphere.

To obtain a baseline, we repeated this analysis for 200 random regions of the same size as AC in each subject and hemisphere. Random regions were formed by randomly selecting a point on the cortex and then selecting nearby vertices (according to geodesic distance) until the new region had the same number of vertices as AC for that subject and hemisphere.