Structural connectivity predicts functional activation during lexical and sublexical reading

A critical question in neuroscience is the extent to which structural connectivity of the brain predicts localization of brain function. Recent research has suggested that anatomical connectivity can predict functional magnetic resonance imaging (fMRI) responses in several cognitive domains, including face, object, scene, and body processing, and development of word recognition skills (Osher et al., 2016; Saygin et al., 2016). However, this technique has not yet been extended to skilled word reading. Thus, we developed a computational model that relates anatomical connectivity (measured using probabilistic tractography) of individual cortical voxels to fMRI responses of the same voxels during lexical and sublexical reading tasks. Our results showed that the model built from structural connectivity was able to accurately predict functional responses of individual subjects based on their structural connectivity alone. This finding was apparent across the cortex, as well as to specific regions of interest associated with reading, language, and spatial attention. Further, we identified the structural connectivity networks associated with different aspects of skilled reading using connectivity analyses, and showed that interconnectivity between left hemisphere language areas and right hemisphere attentional areas underlies both lexical and sublexical reading. This work has important implications for understanding how structural connectivity contributes to reading and suggests that there is a relationship between skilled reading and neuroanatomical brain connectivity that future research should continue to explore.


Structural connectivity predicts cortical activation during lexical and sublexical reading
A fundamental question in cognitive neuroscience is the extent to which structural brain connectivity contributes to cognitive processing. Of particular interest to the current work, reading is a relatively recent human invention that, unlike other cognitive skills (such as speech), requires effortful, explicit instruction in order to be successful. Reading ability has been shown to rely on dissociable, but overlapping, functional and structural networks of brain regions that subserve lexical and sublexical processing. Interestingly, however, while there might be great variability in reading instruction, writing systems, and even processing modality (e.g., alphabetic scripts versus Braille), the cognitive and neural architecture of the reading network appears to develop approximately the same across individuals (e.g., Perfetti, 2011;Rueckl et al., 2015). This suggests that there is a consistent underlying functional and structural neural architecture that potentiates the development of skilled reading. However it is also clear that, even amongst skilled readers, there is heterogeneity of reading processes (Andrews, 2012) that are associated with consistent differences in neural activation (Welcome and Joanisse, 2012). This suggests that there is not a simple, unitary definition of skilled reading, and, instead, that skilled reading may come in many different forms.
Functionally, orthographic processing is associated with a ventral occipitotemporal circuit of brain regions consisting of the left inferior temporal gyrus, lateral extrastriate regions, and left fusiform gyrus that encompass the anatomical 'visual word form area' (VWFA; e.g., Dehaene and Cohen, 2011;Glezer et al., 2009;Stevens et al., 2017;Taylor et al., 2013). Lexical reading has been shown to recruit this ventral system, whereby whole-word reading processes can be optimally promoted through reading of exception words (EWs; words with irregular spelling to sound correspondences, e.g., 'bowl'; Borowsky et al., 2006;Borowsky et al., 2007;Cummine et al., 2015Cummine et al., , 2013. Further, sublexical reading has also been shown to activate the VWFA, as well as recruit a dorsal temporoparietal circuit consisting of the angular and supramarginal gyri and posterior superior temporal gyrus (Taylor et al., 2013). Sublexical reading can be promoted through reading of pseudohomophones (PHs; non-words that when decoded phonetically sound like real words, e.g., 'bohl'; Borowsky et al., 2006Borowsky et al., , 2007, which have similar phonology and semantics as real-words. Phonological representation is associated with the posterior superior temporal gyrus, angular gyrus, and supramarginal gyrus as well as the left inferior frontal gyrus (IFG, i.e., Broca's area; Taylor et al., 2013).
The importance of visuospatial attention in reading has recently been stressed. For example, during reading development, spatial attention plays an important role in remediating reading impairments, whereby spatial attentional training (via action videogames) leads to significant improvements in reading ability (Franceschini et al., 2012;Franceschini et al., 2017). In skilled readers, Ekstrand et al. (2019aEkstrand et al. ( , 2019b found evidence that lexical and sublexical reading strategies are differentially associated with the attentional orienting regions outlined by Corbetta and Shulman (2002). Specifically, lexical reading relies more strongly on ventral, reflexive attentional orienting areas (i.e., the right temporoparietal junction, TPJ), whereas sublexical reading relies more strongly on dorsal, voluntary orienting areas (i.e., the right superior parietal lobule, SPL, and intraparietal sulcus, IPS). Thus, spatial attentional regions of the brain appear to play an integral role in both reading development and skilled reading.
Skilled reading relies on adequate communication between these regions, as well as visual encoding and motor output regions via a network of white matter pathways, which can be reconstructed using diffusion tensor imaging (DTI). Previous research has shown that anatomical connectivity can be used to predict fMRI activation for several cognitive processes, including face, object, scene, and body processing. Seminal work by Saygin et al. (2012) examined the ability of voxel-wise DTI connectivity to predict face selectivity in the right fusiform gyrus. Results from this study were robust and suggested that fMRI activation to faces could be accurately predicted from DTI connectivity. DTI connectivity also predicted fMRI activation better than group fMRI average models, suggesting increased sensitivity of this technique to identify individual differences in task-based fMRI activation. Osher et al. (2016) extended these findings to four visual categories (faces, objects, scenes, and bodies) across the entire brain. Their results indicated that models built from DTI connectivity outperformed group fMRI average models (whereby the researchers argue that group average data is currently the only alternative means of predicting voxel-wise neural responses in a new participant) and were able to successfully predict functional responses across the four visual processing categories. This suggests that individual DTI connectivity can be used to predict brain responses to cognitive functions.
Of particular interest to the current study, previous research has suggested that voxel-wise anatomical connectivity is strongly associated with reading ability, particularly during reading development. Saygin et al. (2016) examined white matter connectivity to the VWFA in children pre-literacy (age 5) and after reading acquisition (age 8) to see if early VWFA connectivity could predict subsequent reading acquisition. To achieve this, the researchers identified the location of the VWFA at age 8 and created a model that utilized the DTI connectivity of the child at age 5 to predict activation in the VWFA at age 8. Their results showed that even prior to functional selectivity in the VWFA for words (i.e., at pre-literacy), there is a distinctive pattern of structural connectivity that is able to predict subsequent reading regions. However, the relationship between voxel-wise anatomical connectivity and skilled reading has yet to be explored using a computational modeling approach similar to Osher et al. (2016). Examination of the relationship between structure and function of skilled reading will be essential for identifying biomarkers of skilled reading, which in turn may inform the assessment of literacy skills and interventions.
Thus, the current study seeks to investigate the extent to which underlying white matter connectivity (measured via DTI tractography) is able to predict fMRI activation during both lexical reading and phonetic decoding in skilled readers. To do this, we will use a similar technique to Osher et al. (2016) that models the relationship between whole-brain structural DTI connectivity and task-based fMRI activation during both lexical and sublexical reading. In line with Saygin et al. (2016), we hypothesize that there will be a strong relationship between structural connectivity profiles and brain function that will generalize to reading in skilled readers. Therefore, using computational modeling techniques, individual structural connectivity should predict fMRI activation in reading tasks, particularly in areas such as the left fusiform gyrus (including the anatomical VWFA), IFG (i.e., Broca's area), and supramarginal and angular gyri, as well as spatial attentional areas that may contribute to reading from Corbetta and Shulman's (2002) model, including the right TPJ, SPL/IPS, IFG, and frontal eye field (FEF). Further, we will be able to examine which connections are significant for predicting task-based activation, thus uncovering structural networks associated with lexical and sublexical reading.

Participants
Thirty participants (mean age 27.1, 15 males) performed DTI scans and lexical (EW) and sublexical (PH) reading tasks during fMRI. All participants spoke English as their first language and reported normal or corrected-to-normal vision. The participants gave written informed consent to participate in the study and all testing procedures were approved by the University of Saskatchewan Research Ethics Board and have therefore been performed in accordance with the ethical standards laid down in the 1964 Declaration of Helsinki and its later amendments. Raw data were generated at the University of Saskatchewan. Data and relevant code from this study are available upon direct request by contacting the corresponding author, R.B.

DWI acquisition parameters and tractography
All imaging was conducted using a 3T Siemens Skyra scanner. Wholebrain anatomical scans were acquired using a high resolution magnetization-prepared rapid acquisition gradient echo (MPRAGE) sequence consisting of 192 T1-weighted echo-planar imaging (EPI) slices of 1-mm thickness (no gap) with an in-plane resolution of 1 Â 1 mm (field of view (FOV) ¼ 256; repetition time (TR) ¼ 1900 ms; echo time (TE) ¼ 2.08 ms).
DTI data were acquired using 195 EPI slices of 4-mm thickness (no gap) with an in-plane resolution of 1.72 Â 1.72, (FOV ¼ 220; TR ¼ 3700 ms; TE ¼ 95 ms; diffusion weighting isotropically distributed along 60 directions; b-value 1000 s mm À2 , with a b 0 volume interspersed every 10 diffusion directions). The top two coil sets (16 channels) of a 20-channel Siemens head-coil were used, with the bottom set for neck imaging (four channels) turned off. Preprocessing included alignment to the b 0 images using FSLs eddy-correct tool (http://fmrib. ox.ac.uk/fsl) to correct for head motion and eddy current distortions, removal of non-brain tissue using the Brain Extraction Tool (BET) from FSL (Smith, 2002), and registration to the high resolution anatomical (T1-weighted structural) scans using FSLs flirt (FMRIB's Linear Image Registration Tool; Jenkinson, Bannister, Brady and Smith, 2002;Jenkinson and Smith, 2001). Next, the GPU version of FSLs bedpostx (Bayesian Estimation of Diffusion Parameters Obtained using Sampling Techniques; Hern andez et al., 2013), ran on a NVIDIA GTX 1070 GPU with 8 GB of RAM, was used to build sampling distributions of the diffusion parameters at each voxel necessary for probabilistic tractography.
Tractography proceeded as follows. First, each of the 268 regions from the Shen et al. (2013) atlas were transformed into diffusion space using FSLs flirt and were checked and corrected for registration errors (if necessary). The DTI-registered parcels were then used as seed and target regions for fiber tracking. Fiber tracking was performed using the GPU version of FSLs probtrackx tool (Hernandez-Fernandez et al., 2016), which uses probabilistic tractography to create a connectivity distribution at each voxel in the seed region (5000 streamline samples per voxel) to each of the target regions, with the distance correction option. This procedure results in a vector of connection probabilities from each voxel in the seed region to all other brain regions.

FMRI protocol
For each of the functional tasks, T2*-weighted single shot gradientecho EPI scans were acquired using an interleaved ascending EPI sequence, consisting of 65 vol of 25 axial slices of 4-mm thickness (1-mm gap) with an in-plane resolution of 2.65-mm Â 2.65-mm (FOV ¼ 250) using a flip angle of 90 . The top two coil sets (16 channels) of a 20-channel Siemens head-coil were used, with the bottom set for neck imaging (four channels) turned off. Acquisition slices were positioned to prioritize complete coverage of the cortex. Additional foam padding was used to reduce head motion. In order to acquire verbal behavioral responses, we used a sparse-sampling (gap paradigm) fMRI method that allows the participant to respond during a gap in image acquisition (TR ¼ 3300 ms, with a 1650 ms gap of no image acquisition; TE ¼ 30 ms; flip angle ¼ 90; e.g., Borowsky et al., 2007Borowsky et al., , 2013. Participants responded vocally (and were instructed to minimize jaw and mouth movements) during the regular, periodic 1650 ms gap in the image acquisition that followed the offset of each volume of image acquisition, which allowed the participants to respond with no noise interference from the MRI. Lexical (EW) and sublexical (PH) reading tasks occurred in two separate runs in the same session.

Stimuli and procedure
Stimuli were presented using a PC running EPrime software (Psychology Software Tools, Inc., http://www.pstnet.com) via MRI compatible goggles (Cinemavision Inc., http://www.cinemavision.biz). Continuous synchronization between the MRI and the experimental paradigm was maintained by detection of the leading edge of the fiberoptic signal emitted by the MRI by a Siemens fMRI trigger converter at the beginning of each acquisition volume that was then passed to the EPrime PC via the serial port. The order of the EW of PH reading conditions was counterbalanced between participants.
Reading tasks. The trial progression for each of the reading tasks was as follows. Participants were presented with 30 target stimuli (either EWs or PHs depending on the task) in a random order with five stimuli presented per block (for a total of six blocks), interspaced with 6 blocks of relaxation. The order of the EW and PH tasks was counterbalanced between participants. The list of stimuli can be found on OSF (https://osf. io/grn9c/). A black central fixation cross (0.6 in height) on a white background was presented. Following this, a jitter of 100, 200, 300, 400, or 500 ms (presented randomly) occurred before presentation of the EW or PH stimulus. This jitter was included to provide more accurate estimates of activation across conditions by staggering the temporal relationships between trial types, thus sampling different components of the hemodynamic response function (Amaro and Barker, 2006).
Participants were asked to read the stimulus aloud as quickly and accurately as possible during the gap in acquisition when the stimulus was presented (1650 ms). EW and PH stimuli were matched on several of the characteristics available from the E-Lexicon Database (http:// elexicon.wustl.edu/), specifically length (t(48) ¼ 0.436, p ¼ 0.665) and log10 base word frequency (t(48) ¼ -.176, p ¼ 0.861). In line with Ekstrand et al. (2019aEkstrand et al. ( , 2019b, phonetic decoding of PHs was used to examine sublexical reading as opposed to pseudowords based on PHs identical phonological representation to their word counterpart and identical meaning. Thus, PHs offer the greatest experimental control for examining differences between lexical and sublexical reading by ensuring that differences in activation are due solely to differences in decoding strategy, not phonology or semantics. During relaxation, a central fixation cross was presented on the screen.

FMRI analyses
Prior to analysis, the functional scans for the EW and PH reading tasks were merged (i.e., concatenated) across time to create a single functional volume, whereby the PH trials were appended to the EW trials, in order to allow for comparison of EW versus PH reading at the individual subject level. FMRI preprocessing and analysis was performed using FSL's FEAT (FMRI Expert Analysis Tool) protocol Version 6.0 (FMRIB, Oxford, UK, http://www.fmrib.ox.ac.uk/fsl/). Preprocessing included MCFLIRT linear slice-time/motion correction (Jenkinson et al., 2002), BET brain extraction (Smith, 2002), spatial smoothing using a Gaussian kernel of full-width half maximum 5-mm, grand-mean intensity normalization of the entire 4D dataset by a single multiplicative factor, high-pass temporal filtering (0.01 Hz; Gaussian-weighted least-squares straight line fitting, with sigma ¼ 16.0s), and normalization to the Montreal Neurological Institute (MNI) 152 T1 2-mm template (however all modeling was performed in the participant's native space). For more accurate registration, the fMRI images were registered first to the high-resolution MPRAGE scan for each participant (6-df linear registration) before registration to the MNI 152 template (12-df linear registration). Functional images were then resampled using 2-mm isotropic voxels and smoothed with a Gaussian kernel of 5 mm full-width at half-maximum into standard space.
Individual subject level comparisons of EW and PH reading, as well as EW greater than PH and PH greater than EW contrasts were performed for each participant, resulting in four t-statistic images. Specifically, firstlevel analyses for each participant were performed using a sinusoidal double-gamma hemodynamic response function convolution that modeled the 5 stimuli over each block versus rest for the EW and PH runs, as well as an additional contrast that modeled the EW blocks versus the PH blocks. Time-series statistical analysis was carried out using FILM with local autocorrelation correction (Woolrich et al., 2001). The t-statistic images were standardized using the same technique as Osher et al. (2016). Specifically, we subtracted the mean functional value of the whole brain from the functional response at each voxel and divided it by the standard deviation of the whole brain. All modeling was then performed on these standardized t-statistic images. We then transformed the whole-brain standardized t-statistic images into DTI space using FSLs flirt by first registering the original functional images to the DTI image and then applying this transformation matrix to the t-statistic images (Jenkinson and Jenkinson et al., 2002). Next, the standardized t-statistic images were masked with the same Shen et al. region masks into 268 anatomical parcels that were the same size as the DTI connectivity images.

Overall group fMRI analysis
We also performed an analysis of the overall group activation patterns for all 30 participants using the first-level analyses for each participant. Group analyses were performed using FSL's FLAME 1 (FMRIB's Local Analysis of Mixed Effects). Results from the whole-brain analyses were thresholded by Z > 3.1 and a corrected cluster significance threshold of p < 0.05 (Worsley, 2001).

Modeling methods
Our modeling approach was comparable to the approach used by Osher et al. (2016) and was implemented using in-house Python code. Participants were divided into two groups whereby modeling for Group 1 (N ¼ 15) was validated using leave-one-out cross-validation (LOOCV) and modeling for Group 2 (N ¼ 15) was performed by applying the final model from all of the participants in Group 1 to each of the participants in Group 2 to evaluate how well the model can generalize to new data. Each participant's anatomical brain was divided into 268 cortical regions using the Shen et al. (2013) atlas in their native space, allowing for individual anatomical variations during modeling.

Group 1
For Group 1, to predict function from connectivity, we employed the LOOCV approach, whereby the connectivity and functional data of a single participant was excluded and a model was trained on the remaining participants before being applied to the left-out subject. We repeated the LOOCV for all participants to create independent predictions for each subject. The modeling proceeded as follows (see Fig. 1): each of the 268 regions from the Shen parcellation for each subject was used as a seed parcel, whereby every voxel of the seed parcel had a functional response to the fMRI contrast (a 1 Â N (number of voxels in the seed regions) vector), as well as DTI connectivity to 267 target parcels (N x 267 matrix, where rows are voxels and columns are connectivity to each of the target regions). Neural responses and DTI connectivity were concatenated (i.e., combined into one matrix, where rows represent voxels across all participants and columns represent connectivity of the voxel to target regions) for all but the left-out participant. Next, the relationship between the fMRI response of a voxel and its DTI connectivity was modeled using a linear regression implemented by the Stats-Models linear regression library of Python language (https://www .statsmodels.org/). This resulted in a 1 Â 267 vector of coefficients relating the relevance of the DTI connectivity from the seed parcel to each of the target parcels for predicting the fMRI response in the seed Fig. 1. Overview of the modeling procedure. Each participant's brain was first divided into 268 regions from the Shen et al. (2013) atlas, as shown by the colored brain in c). Each region was then modeled separately using the following procedure: a) Voxel-wise DTI connectivity from the modeled region to the remaining 267 regions was concatenated for all but one participant (i.e., the left-out participant). b) FMRI t-statistic values corresponding to each of the voxels in a) are concatenated for all but the left-out participant. c) A linear regression (represented by Ä) models the relationship between DTI connectivity in a) and fMRI activity in b). This results in a vector of coefficients (depicted as a greyscale vector) of length 267 (i.e., the number of columns in a) representing each target region) reflecting the contribution of each target region to predicting the fMRI response. d) The left-out participant also has a DTI connectivity matrix with 267 columns that e) the function, f(x), from c) is applied to each voxel in the left-out participants connectivity matrix resulting in f) a vector of predicted fMRI activation for each voxel. Predicted responses are then compared with the actual fMRI responses for each voxel. This procedure is then repeated for each of the other 267 seed regions in c) for each participant, with every participant in Group 1 left out iteratively in order to generate independent predictions for each participant. To predict fMRI activation for Group 2, a final model, f(x), is generated from all of Group 1's voxel-wise connectivity and fMRI data, which is then applied to each participant in Group 2. This entire procedure is repeated for each contrast (i.e., lexical (EW), sublexical (PH), lexical vs. sublexical (EW > PH/PH > EW)). parcel. We then applied these coefficients to the 267 Â N DTI connectivity of the left-out subject, resulting in a predicted fMRI value for each voxel of the left-out participant's seed parcel.
This procedure was repeated for each of the 268 seed regions of the Shen et al. (2013) atlas and concatenated in order to produce predictions for the entire brain of the left-out participant. We then compared the activation predicted by the model to the participant's actual fMRI activation images for each contrast and calculated the absolute error (AE; i.e., the absolute value of actual minus predicted activation) for each voxel. Finally, we created a model using all fifteen participants from Group 1's DTI connectivity and fMRI data that we then applied to the other 15 independent participants in Group 2.

Group 2
The overall model coefficients from all fifteen participants in Group 1 were then applied to an independent group of subjects' (N ¼ 15) individual DTI connectivity data to produce predicted fMRI maps in a similar way to Group 1. We then calculated prediction accuracies by examining the AEs between actual and predicted activation (in the same way as Group 1).

Model validation
In order to assess the validity of our generated models, we compared the performance of our connectivity models to group activation models both across the cortex as well as in regions of interest (ROIs). The SciPy (http://www.scipy.org) Stats module was used to compare mean AEs (MAEs; i.e., the average of all AEs across the brain) in MNI space between the connectivity and group activation models, as discussed below.

Comparison to group activation models
Group activation models were also created using LOOCV using a similar technique to Osher et al. (2016). First, each participant's functional data was transformed to their high-resolution anatomical using FSL's flirt. Next, we performed a nonlinear registration of each participant's anatomical image to MNI space using the symmetric diffeomorphic normalization method for non-linear transformation from Advanced Normalization Tools (ANTs; Avants et al., 2008Avants et al., , 2011. This transformation was then applied to the functional data in anatomical space. All participant's (excluding the left-out participant) fMRI images in standardized space were then superimposed to create composite maps (i.e., the predicted activation for the left-out participant was the average activation from all other participants in the group). We then used this group averaged fMRI image as the input to FSLs FEAT using the same contrasts as those used on the individual participant data. This resulted in t-statistic images for each of the contrasts (i.e., EW, PH, EW > PH, PH > EW) that underwent the same standardization as the t-statistic images used in the DTI connectivity model (i.e., mean functional value of the whole brain was subtracted from the functional response at each voxel and divided by the whole-brain standard deviation). This predicted image was then transformed back into the left-out participant's native space and AEs were calculated. This was repeated for each of the participants in Group 1 to create 15 independent predictive models based on group activation. To calculate MAEs for each participant, whole-brain AE images were then transformed back into MNI using the reverse transforms of those previously applied from FSL's flirt and ANTs.
For Group 2, the group activation model was created from the average activation for all of the subjects in Group 1, the resulting model was transformed to each participant's native space, and AEs and MAEs were calculated in a similar way as for Group 1.

Regions of interest
Each participant's whole-brain AE images were transformed to the MNI 152 T1 2-mm template using FSL's flirt and diffeomorphic normalization method for non-linear transformation from ANTs to ensure all ROIs were comparable across participants. ROIs were derived from parcels in the Shen et al. (2013) 268 region parcellation in the participant's native space and corresponded to the left anatomical VWFA (four regions: Shen 198,199,200,201), IFG (two regions: Shen 151, 156), and temporoparietal regions (two regions: Shen 182,183). Further, based on the proposed integral role of spatial attention in reading and the importance of right hemisphere white matter tracts in reading, we also examined the ability of our model to predict activation in primary spatial attentional regions. Specifically, in the ventral stream we examined the right IFG (two regions: Shen 16,19) and TPJ (two regions: Shen 47,48), and in the dorsal stream the right FEF (Shen 32) and SPL/IPS (Shen 43), resulting in a total of 14 ROIs.
AEs for each region were calculated in a similar way to the wholebrain connectivity versus group average fMRI comparison. A voxelwise paired-samples t-test was performed per participant across all grey-matter voxels in each ROI. A Bonferroni-corrected threshold of p < 0.05/(total number of subjects in both groups times the number of ROIs) (30 Â 14) ¼ 1.19 Â 10 À4 was used to determine the number of participants whose activation pattern was better predicted by the connectivity than the group average activation models. We also calculated MAEs over all of the voxels in each ROI (i.e., the average activation in each ROI) for each participant and performed a paired-sample t-test for each ROI.

Permutation tests
In order to ensure that specific voxel-wise correspondence between activation and connectivity is what is important for accurate prediction and not other factors such as the number of parameters of the model, we also performed permutation testing for each of our ROIs from the final model from Group 1. The procedure for permutation testing was identical to our modeling methods, with the following exception. Prior to modeling, fMRI activation for each seed voxel was shuffled randomly across participants, without altering the DTI connectivity matrix, the linear model was fit and this procedure was repeated 5000 times for each ROI. We then created a distribution of R 2 scores (i.e., the proportion of variance in fMRI activity accounted for by connectivity) from the permutations and evaluated the real DTI model R 2 scores against this distribution.

Functionally relevant DTI networks for reading
In order to uncover the DTI structural connectivity networks associated with each reading contrast, we examined the coefficients from each of the 268 regression models derived from the participants in Group 1. For each anatomical parcel, we calculated the DTI connections that were significant predictors of functional activation for each contrast. Specifically, we assessed significance of the coefficients at a conservative, Bonferroni corrected p-value for multiple comparisons of all beta coefficients across the whole brain of p crit ¼ 7.00 Â 10 À7 (268 regions x 267 coefficients per regression model ¼ 71556 beta weights, p crit ¼ 0.05/ 71556). For the EW and PH contrasts, we examined only the coefficients with positive beta values (i.e., increased connectivity associated with increased activation), however due to the bidirectional nature of our EW > PH contrast, we examined both positive (i.e., increased connectivity associated with increased activation for EW > PH) and negative (i.e., increased connectivity associated with increased activation PH > EW) beta coefficients in separate models.
For visualization purposes, we then binarized the significant coefficients and represented the significant connections as bidirectional in a 268 Â 268 square matrix for each contrast. We did this both across the entire brain (i.e., 268 regions to 268 regions) as well as for the ROIs. Following this, we performed a k-core analysis on both the whole brain and ROI binarized coefficient matrices using the k_core function from NetworkX python toolbox (Hagberg, Schult and Swart, 2008). K-core decomposition is a graph theory technique that can be used to find the underlying 'backbone' of a network by recursively removing nodes until all remaining nodes have a degree (i.e., number of connections) of at least k (see Hagmann et al., 2008 for more information about this technique), whereby we identified the largest k-cores for each network.

Overall fMRI activation patterns across all participants
Results from the whole-brain analysis for each contrast are shown in Fig. 2 and cluster statistics are in Table 1. EWs and PHs activated an extensive network of brain regions that spanned attention and reading areas in both the left and right hemispheres. For the EW versus PH contrast, there were no regions that showed greater activation for EWs than PHs, however PHs showed significantly greater activation than EWs several key areas, most notably the bilateral IPS/SPL, left fusiform gyrus, left IFG.

Predicted neural responses from DTI modeling across all grey-matter voxels
Representative ROI tractography for two participants can be found on OSF (https://osf.io/3bq52/, 'Representative Participant Probabilistic Tractography' folder). Fig. 3 shows the fMRI activation of the most functionally specific voxels (i.e., top 5 percent) for each contrast of both the predicted (from the DTI and group average fMRI models) and actual results for a single participant, in which the predicted results from the DTI model show a strikingly similar pattern of results to the actual fMRI response (particularly in comparison to the fMRI model predictions). This suggests that individual subject activation patterns can be predicted from their DTI connectivity patterns. This is supported by our measures of prediction accuracy, reported below.

Lexical (EW) reading
EW reading typically elicits activation in a ventral occipito-temporal circuit of brain regions that includes the anatomical VWFA in the left fusiform gyrus, lateral extrastriate regions, and inferior temporal gyrus, as well as language and phonological representation areas including the left IFG and TPJ (Taylor et al., 2013;Cummine et al., 2012Cummine et al., , 2015Borowsky et al., 2006Borowsky et al., , 2007. Fig. 3 shows concordance between the predicted and actual responses, particularly in the left fusiform gyrus. Fig. 4a) shows significant correlations between predicted and actual activation for the overall model from Group 1 for all brain regions (see Appendix A of the supplementary material for predicted versus actual scatterplots for each contrast, ROI, and participant).
Comparison to group average activation. Group 1: Averaging across all grey matter voxels for all participants, the model created from DTI connectivity showed significantly lower MAEs than the model created from group fMRI activation, t(14) ¼ À5.01, p ¼ 1.90 Â 10 À4 (see Fig. 5 for means and standard deviations for the EW contrast).

Sublexical (PH) reading
Phonetic decoding has been shown to rely more strongly on dorsal stream regions including the left inferior frontal gyrus (IFG, i.e., Broca's area) as well as the posterior superior temporal gyrus, angular gyrus, and supramarginal gyrus (Taylor et al., 2013;Cummine et al., 2012Cummine et al., , 2015Borowsky et al., 2006Borowsky et al., , 2007. Fig. 3 shows highly similar activation in ventral stream areas including the left fusiform gyrus, and dorsal stream areas including the angular and supramarginal gyri, and posterior parietal regions for the actual vs. predicted results of an example participant. Fig. 4b) shows significant correlations between predicted and actual activation for the overall model from Group 1 across all brain regions.
Comparison to group average activation. Group 1: Similar to the EW task, averaging across all grey matter voxels, the model created from DTI connectivity showed significantly lower MAEs than the model created from group fMRI activation, t(14) ¼ À4.24, p ¼ 8.19 Â 10 À4 (see Fig. 5 for means and standard deviations).

Lexical (EW) > sublexical (PH)
Contrasts between the EW and PH conditions were also assessed in order to determine whether DTI connectivity patterns are able to capture individual differences in lexical versus sublexical processing. Typically, lexical reading has been shown to elicit greater activation than sublexical reading in ventral stream areas including the parahippocampal and fusiform gyri, and middle temporal gyrus (MTG), as well as the posterior cingulum and precuneus, angular gyrus, gyrus rectus, and medial orbitofrontal cortex (Taylor et al., 2013; see also Carreiras et al., 2014;Price, 2012; see Fig. 3). Fig. 4c) shows significant correlations between predicted and actual activation for the overall model from Group 1 across all brain regions.
Comparison to group average activation. Group 1: The model created from DTI connectivity showed significantly lower MAEs than the model created from group fMRI activation, t(14) ¼ À10.56, p ¼ 4.72 Â 10 À8 (see Fig. 5 for means and standard deviations).

Sublexical (PH) > lexical (EW)
Pseudowords have been shown to elicit greater activation than real words in the posterior fusiform gyrus, occipitotemporal cortex, precentral gyrus, left IFG (i.e., Broca's area, supplementary motor area, superior temporal pole, left insula, left parietal cortex, and right inferior parietal cortex (e.g., Taylor et al., 2013). Further, PHs elicit greater activation in the left inferior/superior frontal and middle temporal gyri, left insula, and left SPL than pseudowords (Braun et al., 2015), as well as the angular and supramarginal gyri, and IPL (Borowsky et al., 2006). Fig. 3 shows the voxels with t-statistics in the top 5 percent for the PH > EW contrast across the brain. Because this contrast is the reverse of the EW > PH contrast, the MAE comparison statistics are the same for both groups as the EW > PH contrast.

Connectivity-based predictions of neural responses within regions of interest
Results from our ROI analysis in the most functionally specific regions showed that, when examining MAEs using a paired-samples t-test across all participants, the connectivity-based model outperformed the group average model for the majority of ROIs for each contrast (Table 2). Interestingly, the DTI model was the most accurate at predicting functional activation for the Lexical vs. Sublexical (i.e., EW > PH/PH > EW) contrast, with all but four ROIs showing significantly better prediction accuracies (at a Bonferroni-corrected threshold of p < 3.57 Â 10 À3 ) for the DTI model than the group fMRI average model. Further, when examining voxel-wise comparisons for each participant for each ROI, we found that DTI connectivity models in general outperformed the group average models for the majority of participants for ROIs in both word reading (see Fig. 7) and attention areas (see Fig. 8), with the EW vs. PH contrast showing the greatest number of participants better predicted by the DTI connectivity model than group fMRI activation model.

Permutation tests
Histograms of the permuted and DTI model R 2 values for each ROI from the final model of all Group 1 participants can be found in Appendix Note. R ¼ right, L ¼ left, Z ¼ maximum Z score of the voxel in the cluster. Coordinates are in MNI space. Fig. 3. Representative actual versus predicted (from the DTI connectivity and fMRI models) activation for a single participant. Activation shows the most functionally selective voxels for each contrast (i.e., the top 5% of activation). The top two rows of brains for each contrast are the left hemisphere, the bottom two the right hemisphere.    6. Mean prediction errors and comparison to the group average benchmark model as a function of predictive model (i.e., connectivity versus group fMRI average) for Group 2. Prediction error represents the mean absolute errors across all cortical voxels, error bars represent the standard deviation. Predictions from the connectivity models were significantly more accurate than predictions from the group-fMRI models in all conditions. B of the supplementary materials and on OSF (https://osf.io/3bq52/, 'ROI Permutation Histograms' folder). To summarize, results from permutation tests showed that all ROIs were significant at p < 2.00 Â 10 À4 , with the DTI model R 2 value exceeding all of the permuted R 2 values for all contrasts.

Functionally relevant DTI networks for reading
The matrices of significant beta coefficients, binarized coefficients, and k-cores for each contrast can be found on OSF (https://osf. io/3bq52/, 'Connectivity Matrices' folder), and the binarized coefficient and k-core matrices can be visualized online using the Bioimage Suite Connectivity Viewer (https://bioimagesuiteweb.gith ub.io/webapp/connviewer.html#). This allows readers to explore these extensive networks to examine network architecture in a dynamic way. The k-core network for the whole-brain EW contrast identified 122 nodes with at least 22 connections, 116 nodes with at least 23 connections for the PH contrast, 87 nodes with at least 18 connections for the positive beta coefficient in the EW > PH contrast, and 96 nodes with at least 18 connections for the negative beta coefficient in the EW > PH contrast. The ROI connectivity profiles (i.e., patterns of connectivity using the ROIs as nodes) for each contrast are shown in Fig. 9. Of note, both networks show extensive connectivity between left hemisphere language regions (particularly the anatomical VWFA and the left IFG) and right hemisphere attentional regions in both the parietal and frontal lobes. Further, connectivity for the PH contrast appears to rely more strongly on connectivity between dorsal stream regions than the EW contrast, which may suggest increased attentional involvement during sublexical reading.
Results from our k-core analysis of ROI connectivity is shown in Fig. 10. This analysis identified a core for the significant positive beta coefficients for the EWs of 14 nodes with at least 4 connections. Of note, the right IFG (ROIs 16 and 19) showed connectivity to the right TPJ and putamen, and left IFG and anatomical VWFA; the right SPL/IPS (ROI 43) showed connectivity to the right TPJ and left IFG and anatomical VWFA; the right TPJ (ROIs 47 and 48) showed connectivity to the right and left IFG and right SPL/IPS and precentral gyrus; the left IFG (ROIs 151 and 156) showed connectivity to the right IFG, SPL/IPS, TPJ, and putamen; and the anatomical VWFA (ROI 199) showed connectivity to the right IFG, precentral gyrus, SPL/IPS, and putamen.
The k-core for the ROI PH contrast had a core of 26 nodes with at least 4 connections. Notably, the right IFG (ROIs 16 and 19) showed connectivity to the right SPL/IPS, TPJ, and putamen, and left IFG and anatomical VWFA; the right SPL/IPS (ROI 43) showed connectivity to the right IFG, TPJ, and precentral gyrus, and left IFG; the right TPJ (ROIs 47 and 48) showed connectivity to the right and left IFG and right SPL/IPS and precentral gyrus, and the left anatomical VWFA; the left IFG (ROIs 151 and 156) showed connectivity to the right IFG, precentral gyrus, SPL/IPS, Note. SD ¼ Standard deviation. *Significant at p < 0.05. **Significant at a Bonferroni-corrected threshold of p < 0.05/(number of ROIs) (14) ¼ 3.57 Â 10 À3 . TPJ, and putamen; and the anatomical VWFA (ROI 199) showed connectivity to the right IFG, precentral gyrus, and TPJ. The k-core for the ROI EW > PH contrast of positive beta weights had a core of 14 nodes with at least 4 connections. Notably, the right IFG (ROIs 16 and 19) showed connectivity to the right TPJ, left IFG, and left anatomical VWFA; the right TPJ (ROIs 47 and 48) showed connectivity to the bilateral IFG; the left IFG (ROIs 151 and 156) showed connectivity to the right IFG and TPJ; and the left anatomical VWFA showed connectivity to the right IFG. The k-core for the ROI EW > PH contrast of negative beta weights (i.e., PH > EW) had a core of 11 nodes with at least 4 connections. The right IFG (ROIs 16 and 19) showed connectivity to the right precentral gyrus and TPJ, and left anatomical VWFA; the right TPJ (ROIs 47 and 48) showed connectivity to the right precentral gyrus and IFG; the left IFG (ROIs 151 and 156) showed connectivity to the right precentral gyrus; and the left anatomical VWFA showed connectivity to the right precentral gyrus and IFG.

Predicted neural responses from DTI modeling across all grey-matter voxels
Together, our findings suggest that anatomical connectivity predicts fMRI activation during the cognitive process of reading. When examining MAEs across all cortical voxels we found that predictions from the DTI connectivity model were significantly more accurate than predictions from the group fMRI activation model across all contrasts. This suggests that voxel-wise fMRI activation of an individual across the cortex can be predicted using only their structural connectivity. This corroborates the findings of Osher et al. (2016) suggesting that structural connectivity fingerprints, in part, dictate functional activation, and extends them into the domain of skilled word reading. These results were found not only for LOOCV participants in Group 1, but also for an independent group of subjects (Group 2). Results from each of our contrasts show that models created from DTI connectivity outperformed those created from group-average fMRI activation across the entire cortex and thus were better able to predict voxel-wise fMRI activity during both lexical and sublexical reading. A particularly exciting finding was that our connectivity model was also sensitive to differences between lexical and sublexical processing, suggesting that this technique is sensitive to detecting the differential structural networks that underlie lexical and sublexical processing.

Connectivity-based predictions of neural responses within ROIs
We also examined the performance of our connectivity-based model at predicting neural responses in ROIs. First, we examined ROIs in brain areas known to be involved in reading and language, which included the left fusiform gyrus (i.e., VWFA), IFG, and TPJ and showed that, in general, our model was able to accurately predict fMRI activation in these regions to a similar degree, or better, than a group average fMRI model. As the integral role of visuospatial attention in word reading has recently been stressed in the research literature (e.g., Ekstrand et al., 2019aEkstrand et al., , 2019bFranceschini et al., 2012Franceschini et al., , 2017, we also chose to examine whether connectivity models could better predict reading task-based fMRI activation than group activation models in ROIs associated with spatial attentional processing in the right dorsal (SPL/IPS and FEF) and ventral (TPJ and IFG) spatial attention streams. Our results showed that reading related activation could be accurately predicted in the majority of these spatial attentional ROIs, thus highlighting the importance of spatial attention in reading, as activation in these regions can be accurately predicted from models of lexical and sublexical reading. Together, these results show that structural connectivity fingerprints to ROIs associated with reading and language play an important role in dictating subsequent fMRI activation during lexical and sublexical reading.
Results from permutation testing showed that estimates of fit (i.e., R 2 ) from the DTI model were robust, exceeding those found when the fMRI voxelwise data was randomized. This is an important finding, as one critique of ours and Osher et al.'s (2016) method as a whole is that the connectivity model has a much larger number of parameters compared to the fMRI model, and thus this discrepancy may be contributing to increased flexibility of the DTI model to fit the data in comparison to the fMRI model. Although shuffling the voxel-wise fMRI activation may potentially destroy hidden structure of the data, it still provides valuable insight into the importance of voxel-wise correspondence between activation and connectivity. Specifically, if participants have similar activation patterns in the ROIs that is driving model performance with this high number of parameters (and not individual differences), we would not expect the DTI connectivity model to outperform the permuted models (because fMRI activation would be similar across participants). However, this does not appear to be the case, suggesting that specific voxel-wise correspondence between activation and connectivity is driving the performance of our model, particularly in these small, localized ROIs.

Structural connectivity networks of skilled reading
In order to identify which connections are important for predicting fMRI activation from DTI connectivity, we examined the significant betaweights from the connectivity model. Results from the whole-brain analysis uncovered a complex network of DTI connectivity that underlies skilled reading that spans both hemispheres of the brain, including the anatomical VWFA, speech production regions, and voluntary and reflexive attentional orienting areas. When examining patterns of connectivity associated with our ROIs, several interesting patterns emerged. First, PH reading appears to recruit more dorsal stream regions than EW reading, whereby the connectivity networks for the PH contrast showed greater interconnectivity between the right frontal and parietal cortices than the EW contrast. Second, both types of reading appear to rely on connectivity between the anatomical VWFA and right hemisphere attentional regions, suggesting that direct connectivity to the attentional network is an important facet of both lexical and sublexical reading. Third, when examining the core of the ROI network for each contrast, we found that the right putamen is an important part of lexical and sublexical reading. Connectivity from this region to speech production and spatial attentional areas was found to be a core feature of the EW and PH contrasts, supporting findings that purport that the putamen plays an important role in phonological output (e.g., Gould et al., 2017;Gould et al., 2018;Oberhuber et al., 2013;Seghier and Price, 2010).

Implications for understanding networks of skilled reading
These results extend the work of Saygin et al. (2012) that examined face processing in the fusiform gyrus, and Osher et al. (2016) that examined face, body, scene, and object processing across the cortex, into the processing domain of skilled reading. They are also in concordance with the findings of Saygin et al. (2016), which showed that connectivity to the VWFA (even prior to reading development) can predict subsequent fMRI activation in that area. We extend these findings not only to the entire cortex (i.e., by modeling 268 different brain regions that span cortical grey matter), but importantly, to skilled, adult readers. In addition, the atlas used in our experiment (i.e., the Shen et al., 2013 268 node parcellation) provides a higher resolution parcellation than the Destrieux 148 node atlas (Destrieux et al., 2010) used by Osher et al. (2016), thus providing a more fine-grained examination of structure and function. Results from our ROI analysis of the left fusiform gyrus show that distinct connectivity patterns exist in adulthood that account for a significant amount of variance in reading fMRI activation in these regions. As this region is critically important for reading, this suggests that distinct structural connectivity patterns to this region underlies word processing. Further, we also show that other regions integral to language and reading (i.e., the left IFG and TPJ) have specific connectivity profiles that allow for accurate prediction from DTI connectivity.
A particularly exciting finding from this study is the ability of our model to predict relatively subtle differences in lexical and sublexical processing (i.e., the EW vs. PH contrast). Although there is some specialization for lexical and sublexical reading processes, there is also a large amount of overlap in the reading and language networks (for example, in phonological representation areas, as well as visual and semantic processing regions), particularly between real words and PHs. Comparatively, Osher et al. (2016) examined four different object categories that each have distinct specialization across the cortex and, in particular, the fusiform gyrus (i.e., faces, objects, scenes, and bodies). Further, processing these different object categories occurs naturally and does not require explicit instruction. Thus, the presence of unique structural connectivity patterns underlying these different types of object processing may be more intuitive from an evolutionary perspective.
In contrast, the distinction between lexical and sublexical reading is much more subtle, yet our connectivity model was able to identify different structural connectivity patterns that underlie this distinction.
Thus, there appears to be unique underlying architecture that subserves lexical versus sublexical processing. Based on our connectivity analyses, the core of this network appears to rely on connectivity between the left anatomical VWFA and right hemisphere attentional regions (particularly the rTPJ and rIFG), as well as attentional regions and phonological output areas (i.e., the left IFG). It is possible that this lexical versus sublexical system may be based on the structural development of other cognitive networks (consistent with Dehaene and Cohen's, 2007, cortical recycling hypothesis, whereby new cognitive processes overtake evolutionarily older brain circuits), including those for spatial attention.
Our findings support the idea that reading is reliant on adequate Fig. 10. Core of the ROI connectivity network as determined by k-core decomposition. Lines represent the binarized significant beta coefficients from the DTI model. development of underlying structural connectivity to regions that make up the language and reading networks. This conclusion is supported by the work of Vanderauwera et al. (2018) and Wang et al. (2016) who found that pre-reading tract integrity is an important predictor of subsequent reading outcomes, as well as Saygin et al. (2013) who found that white matter tract volume in key language pathways played an important role in reading development. Further, our model provides exciting insights into the nature of reading impairments by uncovering patterns in structural connectivity associated with skilled reading in adult readers that can serve as biomarkers for identifying reading deficits. Thus, our model may have the potential to help to identify those at risk for reading impairments based on their early structural connectivity fingerprints (similar to Saygin et al., 2012), which may have implications for targeting remediation strategies (e.g., through spatial attentional training, see Franceschini et al., 2015Franceschini et al., , 2017). Our connectivity model was also able to successfully predict activation in known attention areas in the dorsal and ventral stream, lending support to the idea that spatial attention is an integral component of reading. This is in concordance with the findings of Ekstrand et al. (2019aEkstrand et al. ( , 2019b showing that spatial attention is differentially associated with reading strategy (i.e., lexical versus sublexical). Our results provide evidence that the underlying structural network architecture to attention related regions predicts the involvement of attentional orienting regions (as indexed by fMRI activation in reading tasks) during lexical and sublexical reading. This supports research that has found that white matter connectivity in the right hemisphere plays an important role in reading (Catani and Mesulam, 2008;Horowitz-Kraus et al., 2014), as well as research using spatial attentional training as an effective reading intervention (Franceschini et al., 2015(Franceschini et al., , 2017. Further, based on our model coefficients, connectivity to right hemisphere attentional regions in the right frontal and parietal cortices was a core feature of the skilled reading networks, particularly their connectivity to the anatomical VWFA. Thus, future research should continue to examine the role that right hemisphere connectivity and the attentional system play in skilled reading.
It is important to note the purpose of this work is not to argue that connectivity models should replace group average fMRI analyses. Rather, we propose that modeling the relationship between connectivity and function provides valuable insight into the structural network architecture associated with specific cognitive activation patterns, and thus is highly complementary to group level functional analyses. We believe that analyses such as this will further advance exciting and novel lines of research focused on inferring task relevant structural network connectivity using network analysis approaches. Further, while absolute errors provide one window into the efficacy of DTI connectivity models for predicting brain function, future research should assess the efficacy of DTI connectivity models for predicting functional activation using other measures. In addition, this research only scratched the surface of the myriad of information that can be gleaned from complex task-related structural networks, and future research should incorporate additional forms of network and graph theoretical techniques to further characterize these networks.
Importantly, these connectivity-based models may better account for individual variability in fMRI activation than group models, which has implications for both basic research and clinical application. This technique of modeling DTI connectivity with task-based fMRI activation may help to uncover characteristic structural connectivity associated with specific cognitive functions that accommodates individual differences. For example, recent research has shown that skilled readers (who do not differ in reading ability) can be clustered into separable groups based on their brain's response to written stimuli (Fischer-Baum et al., 2018). Thus, even within neurotypical populations there is individual variability in processing, which is possibly accounted for by consistent differences in structural connectivity networks. Therefore, future research should examine whether differences in DTI connectivity can systematically account for individual differences in fMRI activation. This technique could also be used to develop universal models relating brain structure to function using large databases such as the Human Connectome Project (https://www.humanconnectome.org) to characterize consistent patterns of structural connectivity that underlie specific neural responses. Further, it may provide a valuable clinical tool for uncovering language and reading networks in patients for whom functional imaging cannot be performed (e.g., patients who require sedation in the MRI, who are unresponsive/comatose, or are unable to perform the tasks required for functional scanning) from their DTI connectivity. Future research should assess the efficacy of these models for predicting functional brain responses in patient groups, including those who may have irregular or compromised network connectivity.
In conclusion, we show that brain activation during both lexical and sublexical reading in skilled readers can be accurately predicted using DTI connectivity. This finding extends to known reading and language areas including the left IFG (i.e., Broca's area), left TPJ, and the anatomical VWFA, as well as important spatial attentional areas including the right TPJ, IFG, IPS/SPL, and FEF. Further, we identified the structural connectivity networks associated with different aspects of skilled reading using connectivity analyses, showing that interconnectivity between left hemisphere language areas and right hemisphere attentional areas underlies both lexical and sublexical reading. This research broadens our understanding of the structural connectivity fingerprint that underlies skilled reading and has important implications for understanding reading impairment. It may also have clinical implications for aiding localization of language and reading function in patients where functional neuroimaging is not possible. Thus, this research shows that there is a relationship between skilled reading and extrinsic brain connectivity, suggesting that functional organization of reading and language can be determined (at least in part) by structural connectivity patterns. We hope this work will serve as an impetus to examining the structural biomarkers of skilled reading to help broaden our understanding of this essential cognitive process.