Whole-Brain Network Connectivity Underlying the Human Speech Articulation as Emerged Integrating Direct Electric Stimulation, Resting State fMRI and Tractography

Production of fluent speech in humans is based on a precise and coordinated articulation of sounds. A speech articulation network (SAN) has been observed in multiple brain studies typically using either neuroimaging or direct electrical stimulation (DES), thus giving limited knowledge about the whole brain structural and functional organization of this network. In this study, seven right-handed patients underwent awake surgery resection of low-grade gliomas (4) and cavernous angiomas. We combined pre-surgical resting state fMRI (rs-fMRI) and diffusion MRI together with speech arrest sites obtained intra-operatively with DES to address the following goals: (i) determine the cortical areas contributing to the intrinsic functional SAN using the speech arrest sites as functional seeds for rs-fMRI; (ii) evaluate the relative contribution of gray matter terminations from the two major language dorsal stream bundles, the superior longitudinal fasciculus (SLF III) and the arcuate fasciculus (AF); and (iii) evaluate the possible pre-surgical prediction of SAN with rs-fMRI. In all these right-handed patients the intrinsic functional SAN included frontal, inferior parietal, temporal, and insular regions symmetrically and bilaterally distributed across the two hemispheres regardless of the side (four right) of speech arrest evocation. The SLF III provided a much higher density of terminations in the cortical regions of SAN in respect to AF. Pre-surgical rs-fMRI data demonstrated moderate ability to predict the SAN. The set of functional and structural data provided in this multimodal study characterized, at a whole-brain level, a distributed and bi-hemispherical network subserving speech articulation.


INTRODUCTION
Speech articulation is the process of producing individual sounds to compose a word (Levelt, 1993). According to the articulatory loop model by Baddeley (1990), speech articulation is composed by a planning phase in which sounds sequences are retrieved followed by an execution phase which involves the corresponding motor plans to produce desired words. The speech articulatory loop thus constitutes the final and essential stage for encoding articulatory programs into language production. Understanding the network of brain areas involved in the speech articulatory loop is crucial for basic modeling of speech production (Hickok and Poeppel, 2004;Saur et al., 2008;Fridriksson et al., 2016) and its disorders (Duffy, 2016;Boschi et al., 2017;Brabenec et al., 2017;Cordella et al., 2017), as well as for reliable identification of key speech articulation areas to be avoided during surgical resections .
Speech arrest is defined as discontinuation in number counting but with voluntary movements of tongue being possible (anarthria) (Sanai et al., 2008;. Because a counting task minimally requires to process semantic, syntactic, or phonological information, the detection of this deficit provides a very specific tool to disrupt speech articulation (Price, 2010). The counting test is, in fact, commonly adopted worldwide to map and preserve the speech articulation functional network during surgery for resection of intra-cerebral lesions in critical regions independently from the language lateralization (Pallud et al., 2018). Considering the high-reproducibility of speech arrest localization by direct electrical stimulation (DES) in the ventral pre-motor cortex (VPMC), it is routinely used to assess the amplitude threshold for intra-operative mapping during awake surgery (Freyschlag and Duffau, 2014;van Geemen et al., 2014;. It has been demonstrated that combining brain stimulation with neuroimaging techniques in a "perturb and measure" approach allows mapping distributed brain activity related to the perturbation of a cortical region (Paus, 2005). Therefore measuring patterns of activation related to speech arrest may provide a tool to map the network of functionally and structurally connected cortical regions (hubs) supporting speech articulation (Sporns, 2013).
With respect to functional connectivity of speech articulation, a recent invasive study characterized the temporal evolution of this network using intracranial recording and stimulation on 10 right handed pediatric epilepsy patients responding a tone listening and repetition task (Nishida et al., 2017). The authors recorded cortical-cortical evoked potential in the inferior frontal gyrus (IFG) and superior temporal gyrus (STG) following stimulation of the precentral gyrus (PCG) that caused speech arrest in the investigated hemisphere (five left, five right). In a different non-invasive study, areas subserving speech arrest were investigated with repetitive transcranial magnetic stimulation (rTMS) in combination with language task fMRI (Kononen et al., 2015). The authors found a bilateral role of the PCG for speech arrest with rTMS regardless of subjects' handedness (Kononen et al., 2015). However, the covert (internal speech) language fMRI task yielded brain activation overlapping with the areas of speech arrest only in the dominant hemisphere (Kononen et al., 2015). The limited brain coverage of intracranial recording (a few regions within one hemisphere) together with the choice of a strong language lateralizing fMRI task (internal speech) combined with the TMS measures may have not allowed to detect all the cortical areas recruited for speech articulation and, in particular, to detect whether mirror areas across the two hemispheres are involved. This last hypothesis (i.e., bi-hemispherical organization of speech articulation functional connectivity) is indeed supported by the findings of intra-operative mapping studies that have consistently detected speech arrest on both right and left VPMC in right handed patients (Mayka et al., 2006;Chang et al., 2011;Tate et al., 2014).
With regard to structural connectivity of speech articulation a series of studies have reported that intra-operative subcortical stimulation of both the anterior portion of the superior longitudinal fasciculus (SLF III) and the arcuate fasciculus (AF)cause speech arrest (Duffau et al., 2003;van Geemen et al., 2014;Sarubbo et al., 2016). However, it is not clear which among these two tracts plays the most critical functional role for speech articulation. In a study aiming to provide a subcortical atlas of human brain functions (Sarubbo et al., 2015) the authors compared the distance between the centers of mass of several functional response errors collected at DES and the white matter (WM) bundles included in a DTI Atlas (Thiebaut de Schotten et al., 2011). They found that the center of mass of subcortical speech arrest points collected in 22 patients was closer to the SLF III than the AF. This result suggests that SLF III may provide the main tract subserving speech articulation. However, to prove this hypothesis it remains to be demonstrated that this bundle is the main pathway of structural connectivity between the cortical areas recruited in speech articulation at a single patient level.
Based on the findings of above mentioned studies we hypothesized that the structural and functional speech articulation network (SAN) could be characterized at the level of the whole brain in single-subject data by combining pre-surgical resting state fMRI (rs-fMRI) and diffusion-weighted MRI tractography data together with information about speech arrest sites obtained with DES during routine awake surgery. Therefore, we collected and performed a combined analysis of these data in a group of patients with brain tumor and cavernous angiomas with the aim of: mapping the SAN at a whole brain level calculating the rs-fMRI functional connectivity of the speech arrest sites detected with DES with all the other regions of the brain; determining which among the SLF III and AF provides the major associative tract for speech articulation by comparing their density of cortical terminations in the regions of the SAN; evaluating the extent to which pre-operative rs-fMRI networks obtained with independent component analysis (ICA) can predict the SAN as determined using the intra-operative speech arrest seed. This information, in particular, may be helpful if rs-fMRI would be used for pre-surgical mapping of speech articulation areas.

Participants
MRI and intra-operative DES data were acquired and retrospectively analyzed on seven right-handed patients (four males; age range 36-48 years, mean ± SD = 41.4 ± 5.4 years) submitted to routinely surgical removal of low grade gliomas (4) and cavernous angiomas with DES cortico-subcortical mapping. Table 1 reports age, gender, lesion location and histology (in case of gliomas) for each patient included in the study. All the surgical (i.e., awake surgery with DES) and neuroradiological (i.e., diffusion and rs-fMRI imaging) procedures performed and described in this study are part of the standard protocol in our Center for preoperative planning and resections of lesions harboring eloquent areas. All the participants provided an informed consent to the collection of the data and the use for scientific purposes.

MRI Data
All the patients underwent to a routinely protocol at the "S. Chiara" Hospital of Trento APSS for pre-operative planning for resection of lesions in critical areas, including diffusionweighted imaging (DWI) for tractography and rs-fMRI. Images were acquired using a clinical Optima MR450w GE 1.5 T scanner (GE Healthcare, Milwaukee, WI, United States) equipped with an 8-channel receive head RF coil.
For tractography, a 60-direction DWI scheme was acquired (one acquisition) using a single-shot multislice spin echo-echo planar sequence with the following attributes: 50 slices; square FOV: 240 mm; voxel size = 2.4 × 2.4 × 2.4 mm 3 ; TR/TE = 13000/95.8 ms; flip angle: 90 • ; b values of 0 and 1,000 s/mm 2 . For additional tumor characterization, at the end of the MR protocol a second 3D T1-weighted image was acquired (identical to the previously described) following injection of gadolinium. During the procedure of DES, the spatial coordinates of the speech arrest sites are saved using this post gadolinium structural image as reference.

Resting State fMRI Preprocessing
Resting state fMRI data preprocessing was performed using SPM12, FSL (version 5.09), and custom based MATLAB code (version 8.1.0) following Minati et al. (2014). The preprocessing steps are outlined in Supplementary Figure S1 and consisted of: removal of the first four EPI volumes to allow the signal to reach steady state magnetization, slice timing and head motion correction, median (r = 2) filtering, de-trending with a fourth order polynomial and low pass filtering with second-order Butterworth filter (f < 0.1 Hz). Head motion shift and rotation parameters and the average time series of WM and cerebrospinal fluid (CSF) tissue masks, obtained from the segmentation and co-registration of the pre Gad T1-weighted images to the rs-fMRI data, were temporally filtered as above and regressed out. Volume outliers were inspected (after head motion correction) and removed using the ArtRepair software version Version 5b 1 . A volume was classified as an outlier either if its frame by frame head motion was greater than 0.5 mm or its whole brain average time series value was above 2.5 standard deviation the average global signal intensity and it was replaced by interpolating the values from the preceding and subsequent volumes. Finally before iFC analysis rs-fMRI data was spatially smoothed using an 8 mm Full Width Half Maximum Gaussian filter.

Diffusion Processing and Tractography
The processing of DWI was performed concatenating a step for pre-processing, a step for voxel model reconstruction, and a step for deterministic tractography. The analyses pipeline was implemented using FSL and Dipy, an open source library for the analysis of diffusion MRI data (Jenkinson et al., 2002;Garyfallidis et al., 2014).
Following eddy current and head motion correction, images were co-registered to the structural post-Gad 3D T1 and the diffusion tensor model was applied to reconstruct the main orientation of diffusivity for each voxel. The subsequent step of tracking was based on the Euler Delta Crossing (EuDX) method, with 106 seeds. The resulting whole-brain tractograms consist of approximately 150 thousand streamlines. The bundles included in this study (i.e., SLF III and AF) were tracked with TrackVis 2 and a multiple (inclusion/exclusion) regions of interest (ROIs) approach by an expert brain anatomist (author, SS), as part of the pre-surgical planning and for intra-operative support in neuronavigation. For AF an inclusion ROI was placed at the stem (in the anterior and lateral para-trigonal space), posteriorly to the posterior sulcus of the insula and below the posterior thirds of STG and middle temporal gyrus (MTG) and an exclusion ROI was placed at the level of the stem of SLF III, in the lateral WM at the most ventral border between frontal and parietal lobe . Vice versa, the AF stem ROI was used with "no part" logical operator for tracking the SLF III. U-fibers (defined as fibers connecting adjacent gyri) and the posterior indirect component of the SLF were removed.
Surgical Procedure: Intra-operative Functional Mapping and Monitoring All patients underwent asleep-awake-asleep surgery with total intra-venous anesthesia using Remifentanil and Propofol infusion stopped before the awake mapping and neuropsychological monitoring. Cortico-subcortical DES was performed to obtain the most reliable mapping of eloquent structures (Coello et al., 2013). A bipolar electrode with 7 mm spaced tips delivering a biphasic current (pulse frequency of 60 Hz; single-pulse phase duration of 1 ms; amplitude range: 2-4 mA) was used. During the mapping, the intensity threshold was set by evoking a speech arrest (with alternated and regular tongue movements being possible) while the patient was performing a counting task (from 0 to 10), accordingly to previous reports and protocols (Tate et al., 2014;Sarubbo et al., 2015;Riva et al., 2016;. Following speech arrest in all the cases motor mapping was performed for searching for a motor response (facial or tongue) to identify the primary motor cortex and then a full cortical mapping was performed with denomination object and Pyramids and Palm Tree Test (Bello et al., 2006;Coello et al., 2013). In the four patients operated on the right hemisphere after localizing speech arrest sites any other aspect of language responses (i.e., anomia, semantic and phonological paraphasia) was tested. After the cortical mapping the brain tissue resection started, with patients performing continuously the tasks, and it was stopped when the expected functional limits were reached. The different functional response errors collected during intra-operative mapping were associated with a specific tag, stored with high-resolution pictures and subsequently included on the post-gadolinium preoperative T1-weighted images by comparing the anatomical landmarks of the single pictures with the axial volumetric T1-weighted and the reformatted sagittal and coronal images oriented according to the intra-operative lateral position of the head (Tate et al., 2014;Sarubbo et al., 2015). Each functional response was saved as a separate ROI.

Speech Arrest Seed-Based Functional Connectivity Maps
The goal of this analysis was to generate a group intrinsic (i.e., at rest) functional connectivity (iFC) map of the SAN using for each subject the speech arrest sites as ROI seeds for the functional connectivity. For each patient a seed based iFC map of the SAN was generated according to the following pipeline: a 7 mm radius sphere was drawn from the geometrical center of the speech arrest ROI aligned with the T1-weighted post-gadolinium images. The six degree of freedom rigid body matrix to register the T1-weighted post-gadolinium to the T1-weighted (pre-gadolinium) images was determined using FSL-FLIRT. This matrix was subsequently composed with the one calculated to co-register the pre-gadolinium T1-weighted images to the rs-fMRI data and applied it to the 7 mm radius speech arrest sphere. The Pearson's correlation coefficient was calculated between each voxel of the pre-processed rs-fMRI dataset and the speech arrest ROI average time series and converted to z-score using the Fisher transform. For each patient the obtained cross correlation map was thresholded at z = 0.21, p = 0.01 single voxel and family wise error corrected at a significant level of p < 0.05 using cluster-based thresholding based on cluster-size parameters estimated by AFNI's 3dClustsim (Forman et al., 1995).
In order to assess the consistency and variability of the speech articulation iFC across patients, a voxel-wise map of iFC agreement frequency across subjects was created on the Montreal Neurological Institute (MNI) space. To this purpose each subject's thresholded speech articulation iFC map (z > 0.21) was warped into MNI space (FSL T1 template, 2 × 2 × 2 mm 3 ). The co-registration matrix was determined by multiplying the inverse of the matrix obtained to co-register the pre-gadolinium T1-weighted images to the motion corrected rs-fMRI data with a 12 degree of freedom matrix calculated using FSL-FLIRT to register the pre-gadolinium T1-weighted images to the FSL T1 template, 2 × 2 × 2 mm 3 . The quality of the co-registration was visually assessed and deemed adequate for each patient by the first author. Finally, as a measure of SAN consistency across subjects, the proportion of patients for which each voxel of the FSL template brain was included in the speech articulation iFC map was calculated (Seghier and Price, 2016). The average group network consistency within the gray matter (GM) ROIs of the anatomical automatic labeling atlas (AAL) atlas was also derived (Tzourio-Mazoyer et al., 2002).

Relationships Between Functional Regions of SAN and Tractography Reconstructions
The goal of this analysis was to determine, for each subject, the topographic relation between the cortical terminations of the SLF III and AF WM bundles and the subject's specific functional SAN. The underlying question was to test whether one of these two key language tracts was more strongly associated with the articulatory loop than the other. For each patient a 12 parameters affine transformation was determined using FSL-FLIRT to register the T1 post-gadolinium weighted images to the MNI space and applied to the reconstructed SLF III and AF tracts. Following Hau et al. (2016) for each region of the AAL atlas included in the group SAN (iFC cluster sites, see Results) a termination percentage score (TPS) was calculated, dividing the total number of tract streamlines ending in it by the total number of streamlines ending in the GM regions of the AAL atlas (multiplied by 100) for both SLF III and AF. Subsequently, because each patient had only the lesion hemisphere SLF III and AF pre-operatively reconstructed, the left and right portion of each iFC cluster sites were combined and the median TPS across the patients in each of the iFC cluster sites where the frequency of tract terminations was >0 (i.e., at least one subject) was calculated. A Wilcoxon test was run to assess whether the distribution of cortical terminations for the two bundles in the iFC cluster sites was statistically significantly different.

Pre-operative Determination of the Functional SAN
It would be useful if pre-surgical rs-fMRI data could be used to help the neurosurgeon anticipate prior to surgery, the subject's specific functional network related to speech arrest sites. The goal of this analysis was therefore to compare the spatial correspondence between the functional SAN derived from the speech arrest seeds with the ICA networks that were predicted prior to the intervention.
For each patient the preprocessed rs-fMRI dataset was decomposed into 20, 30, and a number of components (range 29-39) determined by probabilistic ICA implemented in MELODIC-FSL (Beckmann et al., 2005). For each component a voxel-wise z-score map was obtained reflecting the agreement between the time series of each voxel with the time series of the specific component. Each z-score map was warped in MNI space using the previously determined registration matrix to warp the seed based iFC maps. Subsequently, following Van Dijk et al. (2009), an automated template matching procedure was applied to select the component that best fit a binary template of the SAN. This template was defined selecting the voxels of FSL template brain that were included in the SAN group frequency iFC map (see Results). Specifically, for each component the difference between the average z-score of all the voxels inside the template and the average z-score of all voxels outside the template was calculated and the one in which this difference (the goodness of fit-GOF) was the greatest was selected. The Spearman's correlation coefficient between the GOF of the component best matching the speech articulation template and the average correlation strength of the speech arrest seed based iFC map inside the template was calculated. Finally, a Kruskal-Wallis's test, was run to assess whether the degree of similarity (highest GOF) between the ICA and seed-based speech articulation iFC maps was affected by the number of ICs used (20, 30, pICA determined).

Intrinsic Functional Connectivity Mapping of the SAN
The speech arrest seed-based iFC maps for each patient are shown in Figure 1. The speech arrest sites were mainly located in right and left frontal inferior and rolandic opercula, in VPMC, according to previous literature findings on the intra-operative detection of anarthria (i.e., interruption of speech articulation not due to tongue movement impairment) (Chang et al., 2011;Tate et al., 2014). We report the coordinates ( Table 2) of each speech arrest ROI center in MNI space (Figure 2). The SAN was consistently and bilaterally observed across all patients in the fronto-ventral, temporal, sensorimotor and inferior parietal FIGURE 1 | Single-subject pre-surgical resting state functional MRI connectivity networks (z > 0.21, hot color scale) using as seed each patient's speech arrest location (green ROIs, Table 2 and Figure 3) obtained during awake surgery. Data are shown in each subject's space. L/R, left/right brain hemispheres.
lobule regions, regardless of whether speech arrest sites were localized in the left (three cases) or right (four cases) brain hemispheres.
In Supplementary Figure S2 the group functional connectivity SAN are shown as voxel-wise frequency maps in MNI space at four different frequency threshold (10, 20, 30, and 40%). We determined that the group maps thresholded at 40% (i.e., at least three or more patients from the group of See also Figure 2 for a visualization of the seeds in MNI space. L/R, left/right brain hemispheres. Frontiers in Human Neuroscience | www.frontiersin.org FIGURE 2 | MNI space localization of each patient's speech arrest ROI obtained during awake surgery using direct electric stimulation (see Table 2 for MNI coordinates). These ROIs were used as seeds in pre-surgical resting state fMRI data to calculate the speech articulation functional connectivity map (Figure 1). The speech arrest ROIs are denoted with a different color for each different patient. R, right hemisphere (four patients) and L, left hemisphere (three patients).
seven, Figure 3) could reliably represent a group functional connectivity SAN because it includes on both hemispheres cortical regions known to be involved in speech articulation (Basilakos et al., 2017) while the additional areas appearing at lower thresholds may be attributed to the inter-subject variability from the various noise sources in rs-fMRI signal (Birn, 2012;Murphy et al., 2013;Power et al., 2015). The GM AAL atlas regions included in this map are listed in Table 3. Group consistency of the functional SAN was observed bilaterally in speech-specific auditory processing (middle and posterior STG), working memory/language articulation (inferior frontal cortex, supramarginal gyrus [SMG]) and motor articulation (sensorimotor cortex) areas. Figure 4 shows a 3D surface representation of the overlap between the SLF III and AF cortical terminations maps (median across subjects) and the group iFC of speech articulation thresholded at 0.4, as in Figure 3. The analysis of the distribution of cortical terminations across the iFC sites overall showed (Figure 5) a much higher density of the SLF III (range 0.22-14.60%) in comparison with the AF (0.30-4.35%), according also to the results of the statistical analysis (Wilcoxon test z = 2.02, p = 0.04, and n = 7). Specifically the SLF III terminated predominantly in the most ventral and anterior parietal cortices (14.6% for SMG and 13.6% for PostCG) followed by the most ventroposterior cortices of the frontal lobe (13.29% for PCG and 8.75% for rolandic opercula). Minor branches projected toward the posterior third of the STG (2.83%), frontal inferior operculum and insula (0.47 and 0.22%). The terminations of the AF in the same regions were instead observed predominantly in PCG (4.35%), Rolandic Opercula (3.59%), and the Insula (2.2%) followed by the PostCG (3.03%) and SMG (0.3%). We excluded from the analysis of the Heschl's gyrus because median PTS was 0 for both SLF (three subject only with cortical termination range 0.07-4%) and AF (one subject PTS = 0.1% with cortical terminations).

Pre-operative Resting State fMRI Determination of the SAN
Prediction of the SAN by pre-surgical rs-fMRI ICA gave promising results as evidenced by the moderate (although not significant due to the small size of our sample N = 7) degree of correlation between the connectivity of the SAN seed regions and the GOF of the speech articulation component obtained with ICA (Figure 6): r = 0.57, p = 0.2; r = 0.57, p = 0.2; and r = 0.61, p = 0.16, respectively, for 20, 30, and pICA determined number of components. Supplementary Figure S3 shows for each patient the z-score map of the component best matching the SAN template obtained including all the voxels in the speech articulation frequency iFC map with a value higher than 40% (at least three or more patients from the group of seven). The similarity between the automatically selected (highest GOF) ICA and seed based speech articulation iFC maps was not dependent on the number of components the rs-fMRI signal was decomposed into (Kruskal-Wallis's test H = 0.32, p = 0.85, dof = 2).

DISCUSSION
This study combines, for the first time to our knowledge, presurgical rs-fMRI and tractography data with intra-operatively defined speech arrest locations to characterize the functional and structural connectivity topography of the SAN. The main findings of the study are: (i) the intrinsic functional connectivity of the SAN includes a set of speech-specific auditory processing (middle and posterior STG), working memory/language articulation (inferior frontal cortex, SMG) and motor articulation (sensorimotor cortex) areas mirrored across the left and right hemisphere; (ii) these cortical regions are structurally connected with higher density by WM SLF III fibers; and (iii) pre-surgical rs-fMRI data provides moderate prediction of the articulatory loop network without knowing a priori the localization of cortical areas that, if stimulated, would impair speech articulation.

Intrinsic Functional Connectivity of the Speech Articulation Network
One of the goals of this study was to characterize at the whole brain level the spatial pattern of the functional SAN, in particular its inter-hemispheric connectivity. To accomplish this, we used a seed-based functional connectivity analysis on pre-surgical rs-fMRI data using as seeds the speech arrest sites determined intra-operatively during standard awake surgery. The cortical areas showing group intrinsic functional connectivity included the PCG, PostCG, the rolandic opercula, pars opercularis of the IFG, SMG, insula, STG (particularly the posterior third), and Heschl's gyrus both left and right hemisphere. Among the set of cortical areas with the highest frequency within the group (>40%), two regions, the SMG in the parietal lobe and the Heschl's gyrus in the temporal lobe had not been previously shown to elicit neuronal response for articulation tasks in fMRI studies (Basilakos et al., 2015(Basilakos et al., , 2017. The recruitment of these two FIGURE 3 | Group speech articulation functional connectivity network shown as voxel-wise group consistency map of spatial overlap in MNI space (color bar represents % of subjects with common SAN voxel). The frequency map of spatial overlap (thresholded at >40%, or at least three subjects from the group of seven) is shown on the MNI cortical surface (upper row) and axial slices (lower row). Regions with highest speech planning network consistency across patients are distributed bilaterally in fronto-temporal and parietal areas. R, right hemisphere and L, left hemisphere. areas, however, had been hypothesized by previous subcortical mapping studies that assigned a role to SLF III for speech articulation because they demonstrated how DES of this bundle induces speech arrest (Duffau et al., 2003;van Geemen et al., 2014;Sarubbo et al., 2015). The Heschl's gyrus and SMG are indeed among the cortical regions connected by the SLF III (Duffau, 2015). Our findings therefore confirm in a multimodal perspective the actual involvement of these two cortical areas in speech articulation.
With regards to the inter-hemispheric intrinsic functional connectivity of the articulatory loop we found a consistent symmetrical bilateral segregation of the network, across right-handers explored, regardless the brain hemisphere of the speech arrest seeds. This evidence suggests, firstly, that the articulatory network is integrated only partially with the classical left-hemisphere dominated language network, as classically emerging from rs-fMRI studies (Liu et al., 2009). The absence or minimum requirement of semantic, syntactic and phonological processing for pure articulation may explain why its activation pattern is not necessarily biased on one hemisphere (Price, 2010). Secondly, our findings imply a role of the right hemisphere in speech articulation. Previous task-based fMRI and intracranial cortico-cortical studies reported activation of regions of the right hemisphere for speech articulation in right handed subjects (Basilakos et al., 2015;Nishida et al., 2017), particularly of the STG and premotor and motor cortices. However, the multimodal approach used in this work (i.e., DES sites based rs-fMRI) allows extending the results from previous studies demonstrating, consistently across subjects, iFC of the articulatory loop between homologous regions across brain hemispheres. These results are also important because they suggest, in agreement with electrophysiology studies of speech production (Brown et al., 2012(Brown et al., , 2014Cogan et al., 2014), that the right hemisphere should be in fact considered part of the articulatory network, even in right-handers. This finding is in contrast with the dual-stream model of auditory language processing predicting that the auditory information is transferred to the left IFG and left PCG through the dorsal pathway (Hickok and Poeppel, 2007). However, the reason why the right hemisphere had not been considered before to be involved in this model could be attributed to the covert language tasks used in the fMRI studies FIGURE 4 | 3D reconstruction of correspondence between the functional connectivity speech articulatory network, SAN (shown in Figure 3, here hot scale) and cortical terminations of two major WM language tracts in MNI space. The upper panel shows the cortical terminations of the SLF III (green). Terminations not overlapping the SAN appear in dark green, while those overlapping are shown in light green. The lower panel shows the corresponding results for the terminations of the AF, in blue. For visualization purposes the degree of overlap was calculated voxel/vertex wise first binarizing the density cortical terminations maps and then adding the speech arrest network frequency value in the same voxel. The results of the quantitative analysis of the cortical termination density maps in the regions of the SAN are presented in Figure 5. the dorsal stream model is based on. Indeed, in right-handed individuals these tasks have been demonstrated to provide strongly lateralized patterns of activation (Zaca et al., 2012).

Structural Connectivity in the Speech Articulatory Network
We found that the cortical regions more frequently mapped in the speech articulatory network are structurally connected predominantly by the SLF III relative to the AF in this patient's series.
On one hand, these findings thus provide an in vivo and multimodal evidence that the SLF III fibers support a functional connectivity pattern for integration of sensorial input and speech output between the temporo-parietal junction (SMG) and the ventroposterior frontal cortices (PCG and Operculi) in both left and right brain hemispheres, as previously reported with DES studies. A role in speech articulation for the SLF III has been hypothesized, in fact, starting from the seminal work by Duffau et al. (2003) demonstrating in a case report that DES of this bundle induced speech arrest. The functional role of the SLF III has been investigated more in detail by recent studies Papagno et al., 2017). Papagno et al. (2017) investigated the neural correlates of verbal short term memory by mapping item and order error in a digit span task and found that stimulation of the SMG and SLF III in the left hemisphere induced more order than item errors whereas stimulation of Broca's areas had the opposite effect. They concluded that the SLF III transfers in the left hemisphere phonological store information from the SMG to speech output areas whereas the AF (not stimulated in this study) would transfer item information from Wernicke's to Broca's area. In a different study post-mortem dissection of SLF III was performed and the authors found that its terminations were located predominantly in the PCG, PostCG, Operculi of the IFG, and SMG . Therefore, in relationship with these previous findings our multimodal study further supports a critical role of the SLF III for speech articulation in both brain hemispheres. In addition, these results can be considered reliable because of the correspondence between the cortical terminations of the SLF III determined from the in vivo diffusion MRI data and from the ex vivo tissue analysis .
On the other hand, our analysis of the overall matching between these bundles' terminations and cortical territories involved in the network is biased by the fact that AF terminate in the temporal and not in the parietal lobe. In this light our results cannot exclude a role of the AF, or at least of the more ventral fibers (Fernández-Miranda et al., 2015;Sarubbo et al., 2016) connecting the ventral and posterior cortices of the frontal lobe with the dorsal and posterior STG, as possible integration of the auditory and phonological information in the SAN.

Resting State fMRI Prediction of the Speech Articulatory Loop
It would be potentially useful, for example in pre-surgical planning studies, to predict the SAN from the intrinsic rs-fMRI networks without prior knowledge of speech arrest sites. To this aim we assessed the degree of spatial correspondence between the pre-surgical resting state networks obtained using ICA (using different choices for the number of IC) and the intrinsic functional connectivity maps seeded in speech arrest sites found FIGURE 5 | Quantitative correspondence between the cortical terminations density (termination percentage score, TPS) of SLF III and AF into the gray matter regions of the SAN cluster sites with highest consistency across patients (>40%, Figure 3 and Table 3). Only cluster sites with density > 0 are represented. See section "Spatial correspondence between the functional connectivity SAN and GM terminations of SLF III and AF white matter" for more details.
FIGURE 6 | Pre-surgical resting state fMRI: similarities between ICA and speech arrest seed-based measures of the group SAN template (Figure 3). For each subject the ICA with the highest goodness of fit to the seed-based network template is defined as well as the average z-score within the template. The correlation between these is shown for various ICA analyses: probabilistic ICA (pICA) and deterministic ICA with 20 or 30 independent components. The highest correlation was found when the number of IC was determined by probabilistic ICA (r = 0.61).
intra-operatively. We found moderate correlation (r = 0.57-0.61) between the two methods for calculating functional connectivity networks. Compared to the findings reported by Van Dijk et al. (2009) performing identical analysis in a study using a 3T scanner, we obtained smaller values for both the GOF of the speech articulation component obtained with ICA and the connectivity of the SAN seed regions. The explanation for this difference may be due to the lower signal to noise ratio and spatial specificity of the BOLD signal of our data acquired using a 1.5 T static magnetic field (Donahue et al., 2011).
The limited agreement between functional connectivity networks derived from ICA and seed-based results has been seen in several other studies (Joel et al., 2011;Rosazza et al., 2012). The way seed-based correlation and ICA determine functionally connected regions across the brain are conceptually different. With seed-based intrinsic connectivity analysis, networks are constructed detecting all the voxels in the brain whose signal is highly temporally correlated with the seed region average time series. In ICA, data-reduction algorithms are used to generate a set of time-courses and spatially independent maps, subsequently classified as noise or functional networks with various manual or automatic techniques (Marchitelli et al., 2016). Furthermore the differences of the results between the two techniques could be also due to a different sensitivity to the physiologically determined systematic fluctuations (Beckmann et al., 2005). Nevertheless, the comparison of these two approaches to obtain intrinsic functional connectivity maps can provide independent confirmation of findings.
To summarize, our results suggest that the SAN could be potentially obtained from pre-surgical rs-fMRI using ICA, although confirmation with a larger sample size is necessary. Furthermore, considering that speech arrest is routinely adopted as the first negative control task to define the current amplitude threshold for the whole cortico-subcortical mapping (Freyschlag and Duffau, 2014), providing SAN maps pre-operatively has the potential to guide and reduce time for intra-operative mapping and consequently to decrease the risks of co-related morbidity (Petrella et al., 2006).

Limitations
The main limitation of this study is the small sample size (N = 7). However, the population of patients we investigated in this work can be considered representative of a larger sample for the following reasons: (a) the balanced distribution of the intra-operatively recorded speech errors across both hemispheres (four right and three left) in right handed patients suggests that the represented bilateral patterns of speech articulation iFC should not be due to a bias in the location of speech arrest points and (b) the inclusion of a mixed sample (i.e., LGGs and CAs), including lesions not harboring cortical and subcortical structures supposed to support this network in order to guarantee that the detected areas functionally connected are not the result of recruitment of new cortical hubs due to plasticity phenomena. In addition, none of the mapped speech arrest points was located in spatial proximity with a lesion. Therefore we do not expect that any of the iFCmap may be affected by false negative due to neurovascular uncoupling or lesion related susceptibility artifacts (Zaca et al., 2014;Agarwal et al., 2016).
In this study, the structural and functional MRI protocols, both acquisition and analysis, were kept constant for all subjects, despite the fact that we study both patients with low grade gliomas and patients with cavernous angiomas. These pathologies have different bleeding patterns and the resulting magnetic susceptibility effects may give different sensitivity to local signal T2 * loss and geometric distortions.
With reference to the analysis of the WM tracts supporting speech articulation we limited our analysis to the cortical terminations of SLF III and AF. Considering the strong inter-hemispherical connectivity of the SAN found in this study, future work may investigate whether inter-hemispherical fibers, such as the anterior third of corpus callosum may have a role in speech articulation . Finally, we used for the tractography analysis a deterministic approach. However, DTI is still now the most diffused algorithm for reconstruction in the clinical and surgical fields. Moreover, our diffusion MRI protocol includes a 60 encoding directions. We preferred to utilize the original tractograms of our patients, rather than merging a healthy tractography atlas to patient data, in order to maintain the patient-specific WM topography for each tumor case.

CONCLUSION
Characterization of the cortical network supporting speech articulation is important to improve our understanding of this fundamental component of the language function involving planning and executing sequences of sounds to form words. Moreover, it may as well help to determine the neural correlates of speech articulation impairment following brain injuries and to minimize the chance of its resection during surgical intervention. In this study we described a multi-modal imaging approach that uses pre-surgical and intra-operative data to map at the whole-brain level, for the first time, the iFC of the SAN. This functional network includes mostly ventroposterior frontal and parieto-temporal cortical territories in both hemispheres, some of them had not been reported in previous task based fMRI studies of speech articulation, such as the SMG and Heschl's gyri. Furthermore, considering the two major WM tracks involved in the dorsal stream language model, AF and SLF III, we determined that SLF III gives the largest contribution of fibers providing the axonal subcortical circuits structurally connecting the cortical hubs of the articulatory loop network. Finally, we found that pre-surgical rs-fMRI data offers a possible prediction of the speech articulatory network. This, however, needs further validation on a larger dataset.
Given the potential utility of these results for pre-surgical planning and for future studies investigating disturbances of the articulatory loop with different neuroimaging techniques, we provide speech articulation spatial frequency maps and speech arrest seeding points (MNI coordinates) for download 3 .

ETHICS STATEMENT
This study was carried out in accordance with the standard preoperative and operative protocol of the Division of Neurosurgery of "S. Chiara" Hospital for resection of brain lesions harboring critical areas. All subjects gave written informed consent in accordance with the Declaration of Helsinki.

AUTHOR CONTRIBUTIONS
DZ processed the data, interpreted the results, and wrote the manuscript. FCo, MD, and GF collected the data, processed the data, and interpreted the results. UR, LA, LZ, and FCh collected the data, processed the data, interpreted the results, and revised the manuscript. PA processed the data, interpreted the results, and revised the manuscript. JJ and SS processed the data, interpreted the results, and wrote and revised the manuscript.

FUNDING
This work was supported by the CARITRO Foundation within the "Bando 2015 per progetti di ricerca scientifica svolti da giovani ricercatori" (University of Trento, Trento, Italy).