Dynamic tractography: Integrating cortico-cortical evoked potentials and diffusion imaging

INTRODUCTION
Cortico-cortical evoked potentials (CCEPs) are utilized to identify effective networks in the human brain. Following single-pulse electrical stimulation of cortical electrodes, evoked responses are recorded from distant cortical areas. A negative deflection (N1) which occurs 10-50 ms post-stimulus is considered to be a marker for direct cortico-cortical connectivity. However, with CCEPs alone it is not possible to observe the white matter pathways that conduct the signal or accurately predict N1 amplitude and latency at downstream recoding sites. Here, we develop a new approach, termed "dynamic tractography," which integrates CCEP data with diffusion-weighted imaging (DWI) data collected from the same patients. This innovative method allows greater insights into cortico-cortical networks than provided by each method alone and may improve the understanding of large-scale networks that support cognitive functions. The dynamic tractography model produces several fundamental hypotheses which we investigate: 1) DWI-based pathlength predicts N1 latency; 2) DWI-based pathlength negatively predicts N1 voltage; and 3) fractional anisotropy (FA) along the white matter path predicts N1 propagation velocity.


METHODS
Twenty-three neurosurgical patients with drug-resistant epilepsy underwent both extraoperative CCEP recordings and preoperative DWI scans. Subdural grids of 3 mm diameter electrodes were used for stimulation and recording, with 98-128 eligible electrodes per patient. CCEPs were elicited by trains of 1 Hz stimuli with an intensity of 5 mA and recorded at a sample rate of 1 kHz. N1 peak and latency were defined as the maximum of a negative deflection within 10-50 ms post-stimulus with a z-score > 5 relative to baseline. Electrodes and DWI were coregistered to construct electrode connectomes for white matter quantification.


RESULTS
Clinical variables (age, sex, number of anti-epileptic drugs, handedness, and stimulated hemisphere) did not correlate with the key outcome measures (N1 peak amplitude, latency, velocity, or DWI pathlength). All subjects and electrodes were therefore pooled into a group-level analysis to determine overall patterns. As hypothesized, DWI path length positively predicted N1 latency (R2 = 0.81, β = 1.51, p = 4.76e-16) and negatively predicted N1 voltage (R2 = 0.79, β = -0.094, p = 9.30e-15), while FA predicted N1 propagation velocity (R2 = 0.35, β = 48.0, p = 0.001).


CONCLUSION
We have demonstrated that the strength and timing of the CCEP N1 is dependent on the properties of the underlying white matter network. Integrated CCEP and DWI visualization allows robust localization of intact axonal pathways which effectively interconnect eloquent cortex.


Introduction
A primary goal of current neurobiological research is to determine the structural and functional networks that support cognition in the human brain. Structural connectivity in the human brain is typically explored with non-invasive diffusion-weighted imaging (DWI)-based methods, including tractography. However, these data are limited in that they provide "static" measures such as the presence or absence of white matter pathways, spatial trajectory, and biomarkers for myelination. Further, current tractography-based network analysis is prone to false-positive tracking of fibers particularly at the regions that multiple fibers intersect (Maier-Hein et al., 2017). As a part of cortical mapping procedures during resective neurosurgery, clinicians in tertiary epilepsy centers use cortico-cortical evoked potentials (CCEPs), an invasive technique in which a single-pulse current is applied to the surface of the cortex to causally establish its effective projections (Matsumoto et al., 2017(Matsumoto et al., , 2004Nishida et al., 2017). By quantifying the electrocortical dynamics downstream from the stimulation site, CCEPs can be used to causally explore the time-varying dynamics of cortical effective network connectivity with high temporal resolution. Single-pulse stimulation of receptive language areas of the left superior temporal gyrus were reported to elicit CCEPs in speech-critical areas of the left inferior frontal and precentral gyri (Matsumoto et al., 2017(Matsumoto et al., , 2004Nishida et al., 2017).
CCEPs are limited in that they do not provide any direct information about the physical pathway taken by the connectivity. Furthermore, current methods to visualize CCEP data at the subcortical level remain suboptimal. In the present study, we sought to investigate a new complementary approach for advancing structural and effective network analysis which integrates data from invasive CCEPs with non-invasive DWI-based tractography to gain deeper insights into the intrinsic connectivity networks that support human cognition. This new "dynamic tractography" approach may provide improved visualization and utilization of the spatiotemporal profiles of cortico-cortical neuronal propagation elicited by electrical stimulation.
CCEP methodology has been employed as an intervention-based approach for mapping cortical networks in neurosurgical patients (Matsumoto et al., 2017(Matsumoto et al., , 2004Matsuzaki et al., 2013;Nishida et al., 2017;Tamura et al., 2016). CCEPs are typically recorded from intracranial cortical surface disk electrodes or stereo-EEG depth electrodes (Prime et al., 2018). During the procedure, a small current is applied to electrode pairs while simultaneously recording from all other electrodes (Matsumoto et al., 2017(Matsumoto et al., , 2004. Stimulation of the cortex can result in evoked potentials at distant electrode sites observable as a complex of negative and positive deflections beginning as early as a few milliseconds post-stimulus (Keller et al., 2014a;Matsumoto et al., 2017). Mapping nonepileptic cortical networks in this way provides a rare and unique window into the dynamics of human cortical networks with a spatiotemporal resolution typically exclusive to animal models. CCEPs consist of mono-and polysynaptic components which can be roughly distinguished by several features. Large, negative voltage deflections, increases in the signal root mean square (Shimada et al., 2017), high gamma (>70 Hz), power (Crowther et al., 2019;Usami et al., 2015), or broadband power (Logothetis et al., 2010) 10-50 ms post-stimulus are referred to as "N1" components and believed to reflect mostly monosynaptic postsynaptic potentials (Logothetis et al., 2010), whereas later negative potential deflections ("N2", 80-200 ms post-stimulus) in the 2-6 Hz range are thought to reflect a subsequent inhibition of action potentials (Alarc on et al., 2012;Keller et al., 2014b;Matsumoto et al., 2004).
Historically, detailed studies of structural cortical networks have been limited to animal tract-tracing studies. In recent decades, DWI-based tractography and associated measures, such as fractional anisotropy (FA), which may be related to neural bundle integrity (Scholz et al., 2013), have been extensively used to estimate properties and trajectories of white matter networks which support cortico-cortical communication. DWI-based tractography and FA have been successfully used to describe cortical networks in both health and disease (Jeong et al., 2015(Jeong et al., , 2011Martin et al., 2016;Park and Friston, 2013;Silverstein et al., 2018). Resection of networks identified via tractography have been shown to increase the likelihood of post-operative deficits in epilepsy neurosurgery (Lee et al., 2019).
However, while DWI can accurately map spatial trajectories of white matter pathways across the brain, it doesn't provide any functional or directional information about these pathways in vivo. Using tractography and structural MRI data alone, it is difficult to differentiate true and false positive pathways (Maier-Hein et al., 2017). CCEPs and DWI provide complementary lenses on structural networks in the brain. Integration of the two may provide a greater depth of insight into the networks which support healthy and pathological cortico-cortical communication. DWI-based tractography can reveal the biological plausible paths taken by the CCEP, while CCEPs, as an interventional technique with high specificity, can be used to filter likely true positives from the tractography data, thereby deepening our understanding of how white matter networks give rise to observed functional dynamics.
Recent interventional evidence has indicated a direct link between white matter pathways and CCEPs. Yamao et al. (2014) directly stimulated the arcuate fasciculus during an intraoperative recording and observed CCEPs at each endpoint which were similar to waveforms generated from cortical stimulation at arcuate endpoints. However, studies directly integrating CCEPs and DWI tractography have been rare. The initial study directly combining DWI tractography and CCEP data demonstrated a positive correlation between streamline counts and N1 peak amplitude and a negative correlation between count and latency along the arcuate (Conner et al., 2011). The relationship between CCEP connectivity and streamline count was later replicated in a study integrating stereo-EEG CCEP data with publicly available atlas-derived tractography (Donos et al., 2016). Similarly, there is a body of work demonstrating the effect of distance on CCEP connectivity. Matsumoto et al. established a correlation between inter-electrode cortical surface distance and N1 latency (Matsumoto et al., 2012). Euclidean distance between recording and stimulation electrodes has been shown to correlate with CCEP latency (Crowther et al., 2019), decreasing connection reciprocity (Keller et al., 2014a), as well as decreasing indegree, outdegree, and connection strength (Entz et al., 2014). Likewise, several studies have shown that DWI-based tractography streamlines are more likely to connect stimulation sites with electrode locations showing evoked potentials than those that do not (Elias et al., 2012;Javadi et al., 2016) as well as convergence between DWI-based structural networks and CCEP-based effective networks (Donos et al., 2016;Parker et al., 2018). A recent paper which related CCEP data from a patient population to a DWI tractography atlas derived from different subjects demonstrated decreased connection probabilities and increased latencies with increased distance as measured by atlas-based streamline lengths (Trebaul et al., 2018). This latest study provides intriguing insights into the relationship between tractography and CCEP connectivity but has yet to be replicated using CCEP and DWI tractography data collected from the same individual.
To improve the interpretability of CCEPs and better understand how CCEP propagation relates to the underlying white matter, we here build on the existing literature to establish a general model of CCEP propagation by integrating diffusion data and CCEP data recorded from the same individuals. We started with the basic hypothesis that the length of a tractography-derived streamline between the stimulation and recording sites should predict the CCEP N1 latency. Next, we investigated the relationship of white matter properties to N1 peak voltage. Insofar as myelination provides axons with insulation, it reduces the signalconduction resistance (Hartline and Colman, 2007). Given that long-range pathways are more likely to be susceptible to loss of phase synchrony across axons due to divergent pathways, this leads to the hypothesis that longer-range streamlines should predict lower N1 voltages. Lastly, we investigated whether the mean FA along a streamline predicts the N1 propagation velocity, as has been demonstrated in some peripheral (Heckel et al., 2015) and central nervous system pathways (Whitford et al., 2011), and sought to explore the spatial distribution of propagation velocities across the brain. To our knowledge, this is the first multimodal imaging case-series study to examine the relationship between streamline length, FA, and CCEP connectivity in the same patients.

Patient data
Patients were selected using the following inclusion criteria: (i) a history of focal epilepsy scheduled for extraoperative subdural intracranial EEG (iEEG) as part of the presurgical evaluation, (ii) patient underwent measurement of CCEPs, (iii) preoperative diffusion and T1weighted MRI sequences. Exclusion criteria consisted of: (i) presence of massive brain malformations affecting the lateral or central sulci (e.g., perisylvian polymicrogyria), (ii) history of previous epilepsy surgery. Twenty-three patients satisfied these criteria (age range: 3-20, mean ¼ 12.7 AE 3.6; 12 female). This study was approved by the Institutional Review Board at Wayne State University and performed in accordance with the approved guidelines. Informed consent/assent was obtained from the patients or guardians of patients. Patient demographics, pathology, and seizure-related data are summarized in Table 1.

Subdural electrode placement
Subdural platinum grid and strip electrodes (10 mm center-to-center distance; 3 mm diameter exposed) were placed as previously described . Electrode plates were stitched to adjacent plates or to the edge of the dura mater to prevent electrode plate movement after placement. iEEG recordings were obtained for 3-5 days, using a 192-channel Nihon Kohden Neurofax 1100A Digital System (Nihon Kohden America Inc., Foothill Ranch, CA, USA) at a sampling frequency of 1000 Hz and an amplifier bandpass of 0.016-300 Hz. The averaged voltage of iEEG signals derived from the fifth and sixth intracranial electrodes of the iEEG amplifier was used as the original reference (Wu et al., 2011). iEEG signals were then re-montaged to a common average reference (Sinai et al., 2005;Wu et al., 2011). Channels within 25 mm of the stimulation site and channels contaminated with artifacts or interictal epileptiform discharges were excluded from the common average reference.
2.3. Cortico-cortical evoked potential (CCEP) stimulation protocol As part of the clinical management, during extraoperative monitoring, single-pulse electrical stimulation (SPES) was delivered to adjacent pairs of subdural electrodes at a frequency of 1 Hz for 40 s by E.A., who was blinded to the DWI findings. Each stimulation consisted of a biphasic square wave pulse with 5 mA intensity and 0.6 ms pulse width including both phases (Grass S88 stimulator: Astro-Med, Inc, West Warwick, RI). Biphasic stimulation was utilized to minimize the chance of charge buildup in cortical tissue. These stimulation parameters have previously been demonstrated to be clinically effective and safe (Koubeissi et al., 2012;Matsumoto et al., 2005Matsumoto et al., , 2004. During a train of stimulations, patients lie comfortably in bed with their eyes closed in a dimly lit room. We did not perform SPES during postictal or non-resting periods. Patients were under Fosphenytoin during the CCEP data acquisition procedure (Arya and Crone, 2018); however all other oral anti-epileptic drugs were stopped or reduced prior to the CCEP procedure. No adverse effects (e.g., afterdischarges) were noted. During each stimulation, cortical activity was simultaneously recorded at all other electrode sites. Stimulation electrodes were removed from the amplifier during stimulation. Across the 23 patients a total of 1021 electrode pairs (24-82 pairs per patient) were stimulated according to the CCEP protocol outlined above. iEEG responses to each stimulation were recorded at 53-128 qualified nonepileptic electrodes, with a total of 3003 across patients (defined as those outside the SOZ, interictal spiking zone, and structural lesions; Fig. 1). Our clinical procedure aims to stimulate all non-artifactual, nonepileptic electrode sites (Motoi et al., 2018) during the rest period. We terminated SPES as soon as a given patient showed a symptom (e.g., twitching of a body part attributed to the excitatory responses following a motor cortex stimulation) and excluded such stimulus sites from further analysis. We do not extend the duration of extraoperative iEEG recording for research purposes.

Preprocessing
iEEG data were first segmented into 1000 ms trials time-locked to stimulus onset and bandpass filtered from 1 to 45 Hz using a zero-phase, hamming windowed, sinc FIR filter (Delorme and Makeig, 2004). The time-locked data were then baseline corrected by subtracting the mean of the À200 ms to À20 ms pre-stimulus period from each data point. Baseline corrected data were then averaged across trials to improve the signal-to-noise ratio.

Peak detection
Trial-averaged evoked potential voltage deflections greater than 5 standard deviations beyond the baseline period were taken as evidence for cortico-cortical connectivity (Trebaul et al., 2018). N1 peak voltage and corresponding latency were determined via the second derivative test within a 10-50 ms post-stimulus range (Fig. 2). We chose to start the N1 window at 10 ms to minimize potential stimulation artifact contamination (Conner et al., 2011;Matsumoto et al., 2017). We performed semi-automated peak detection using a modified form of the second derivative test. In addition to showing a local maximum, peaks were also required to be 5 standard deviations above the baseline period. Two expert board-certified clinical neurophysiologists/epileptologists (E.A. and M.S.) validated the peak detection algorithm results to eliminate the risk of misinterpretation of artifactual signals as true CCEPs.

Connectivity definition
CCEP data are inherently directionala pair of electrodes are stimulated, and evoked potentials recorded at distant individual electrodes. By contrast, tractography is non-directional. Each pair of electrodes therefore has at least two CCEP peak/latency pairs associated with it, whereas each tractography region of interest (ROI) pair has a single property (e.g., AB→C, AD→C vs just A→C). Therefore, to more straightforwardly related CCEP pair quantities with DWI pair quantities, the N1 peak and latency data were rearranged into a standard network connectivity adjacency matrix configuration. "Overlapping" electrode quantities were averaged such that if AB→ C exists and AD→C exists, then A→C ¼ (AB→C þ AD→C)/2. This reflects the logic for determining connectivity and involvement of a cortical region during clinical  Left: CCEP framework. iEEG data are recorded during single-pulse direct current stimulation. CCEPs are detected and quantified based on whether the maximum negative voltage deflection relative to the pre-stimulus baseline period surpasses a z-score threshold (grey dotted line). The N1 amplitude and latency were defined as the peak negative voltage within the 10-50 ms post-stimulus window (pink box). Only peaks that were greater than 5 standard deviations beyond the À200 to À20 ms pre-stimulus baseline were analyzed. Right: Representative CCEP cortical distribution from patient 15 ( Table 1) following stimulation of the right superior temporal gyrus (STG). Colored cortical areas indicate N1 magnitude. Dots indicate iEEG electrode locations. Electrodes with detected N1s are colored by latency (pink). Stimulation site on the STG is indicated by red dots. B.H. Silverstein et al. NeuroImage 215 (2020) 116763 evaluation (Arya and Crone, 2018;Mooij et al., 2018). Electrodes that were within 25 mm of the stimulation site were eliminated from the analysis to limit the effect of stimulation artifacts and potential volume conduction (Conner et al., 2011;Matsumoto et al., 2004;Shimada et al., 2017). To focus the analysis exclusively on nonepileptic tissue, electrodes that were classified as within the SOZ, showed interictal spiking, or structural lesions were eliminated from the analysis.
2.5. DWI analysis 2.5.1. Image acquisition MRI scans were performed on a 3T GE Signa scanner using an eightchannel head coil. DWI images were acquired with a multislice singleshot diffusion-weighted echo-planar imaging sequence at TR ¼ 12,500 ms, TE ¼ 88.7 ms, FOV ¼ 24 cm, 128 Â 128 acquisition matrix, contiguous 3 mm thickness axial slices in order to cover the entire brain using 55 isotropic gradient directions with b ¼ 1000 s/mm 2 , number of excitations ¼ 1, and single b ¼ 0 image. For anatomical coregistration, T1-weighted images were acquired using a fast spoiled gradient-recalled echo sequence at TR ¼ 9.12 ms, TE ¼ 3.66 ms, TI ¼ 400 ms, slice thickness ¼ 1.2 mm, and planar resolution ¼ 0.94 Â 0.94 mm 2 .

ROI definition
Subdural electrodes were coregistered to the T1-weighted images as previously described (Alkonyi et al., 2009;Muzik et al., 2007;Wu et al., 2011). The accuracy of the coregistration was validated by intraoperative photographs showing in situ electrode locations (Dalal et al., 2008). T1-weighted images were skull-stripped and converted to grey and white matter surfaces via FreeSurfer (Fischl et al., 1999). The center coordinate of each electrode was placed on the grey matter surface according to the above-cited methodology Electrodes were placed directly on the cortical surface, but tractography aims to generate streamlines that terminate at the grey matter-white matter interface. To facilitate this, the center coordinate of each electrode location was mapped onto the white matter surface using FreeSurfer (Desikan et al., 2006;Ghosh et al., 2010). Electrode coordinates were then transformed to the DWI space by applying an affine transformation derived from coregistering the DWI b0 images and T1 images using the ANTs library package (Avants et al., 2008). A 4 mm radius sphere was then grown around each electrode coordinate to generate ROIs. As mapping from the grey matter to the white matter surface causes electrode locations to shift, some electrode ROIs showed overlap. Overlapping ROI volumes were deleted and in the case of >50% overlap, one of the two ROIs was dropped (Fig. 3) as the remaining ROI segment was often too small to provide robust streamline tracking results. In case one of the two stimulating electrodes was dropped, only streamlines from the remaining ROI were considered.

Tractography, streamline length, streamline quality control, and fractional anisotropy
Anatomically-constrained whole brain probabilistic tractography (Patenaude et al., 2011;Smith et al., 2012Smith et al., , 2004Smith, 2002;Zhang et al., 2001) was performed using the iFOD2 reconstruction algorithm (Tournier et al., 2010) with a seed density of 3000 seeds/voxel, randomly distributed, minimum/maximum streamline length of 20/250 mm, and a maximum angle of 70 . Streamlines connecting electrode ROIs were then sorted for further analysis (Smith et al., 2015), resulting in electrode connectomes (Fig. 3). Streamline length for a given ROI pair was computed as the minimum length of all sorted streamlines connecting the two ROIs. To quantify the FA associated with each electrode pair, the mean FA across all voxels intersected by the shortest streamline was between each ROI was computed. All DWI tractography generation procedures were performed using the MRtrix3 package (htt p://www.mrtrix.org/). For group-level tractography plots, streamlines were normalized to template space using a nonlinear warp field Fig. 3. Upper: Representative electrode tractography from one patient overlaid on a render of a patient T1-weighted image. Whole-brain probabilistic tractography was estimated in the MRtrix framework. Electrodes locations were coregistered with the DWI image. Middle: Representative complete matrices for N1 amplitude, latency, and streamline length. Lower: Intersection networks computed from the intersection of the three complete matrices. This process leverages the interventional networks (N1 amplitude peak and latency) to reduce the false positive rate of the tractography network. F: frontal, P: parietal, T: temporal, O: occipital. generated by coregistering the patient b0 and T1 images with the Free-Surfer template image using the ANTs library package (Avants et al., 2008).
To minimize the presence of biologically implausible streamlines, we calculated the Euclidean distance between streamline endpoints and statistically inferred the length of the streamline. That is, a group-level linear regression model relating the two properties: distance vs. length was computed. Streamlines with model errors greater than 3 SDs of the mean were eliminated as outliers. This ad-hoc procedure eliminated biologically implausible streamlines, while preserving biologically plausible curvatures, e.g. U-fibers, the arcuate fasciculus, and the uncinate fasciculus, all of which are well-represented in our dataset (Fig. 6).

N1 propagation velocity
The N1 propagation velocity between the stimulation and recording sites was estimated according to the simple formula, velocity equals distance divided by time. Here, to estimate velocity, we divided the minimum streamline length between two electrode ROIs by the N1 latency between them.

Covariates
Due to the potentially heterogenous nature of epilepsy patient data, we first investigated whether factors related to epilepsy or development influenced our outcome measures of N1 peak amplitude and latency. Correlations between the primary outcome measures, N1 peak amplitude, N1 latency, streamline length, and velocity and the clinical covariates, patient age, sex, number of oral anti-epileptic drugs, proximity of the stimulation site to the seizure onset zone, and stimulation hemisphere were minimal to non-existent (Table 2). Therefore, these factors were excluded from the regression models.

True positive network
To focus the model on reliable connections and minimize the influence of tractography false positives, only electrode pairs with significant N1 peaks, interconnecting streamlines, and a mean FA > 0.2 were included in the analysis (Fig. 3). The resulting "intersection networks" were then aggregated across patients for group level analysis.

Regression model
To control potential type I and II errors resulting from tractography tracking algorithm and parameter selection, we performed an intervalwise regression analysis where the N1 peak amplitude and latency were averaged at every 10 mm interval of streamline length. The observed streamlines ranged between 25 and 175 mm, resulting in 15 length bins. To generate a comparable model for N1 propagation velocity, we computed the mean velocity across 15 equal-sized FA bins. The resulting values were used in a regression models predicting CCEP properties from DWI-based properties.

Streamline length predicts N1 peak latency and amplitude
Global analysis of N1 peak amplitude and streamline length revealed that N1 latency and streamline length revealed that N1 latency was weakly predicted by streamline length, where shorter streamlines predicted small latency values (R 2 ¼ 0.21, β ¼ 0.034, p ¼ 0.048). On the other hand, N1 peak amplitude was strongly predicted by streamline length (R 2 ¼ 0.84, β ¼ À0.38, p ¼ 8.22e-7; Fig. 4). These findings confirm our first two hypotheses that streamline length is a strong indicator of cortico-cortical connectivity between nonepileptic cortical areas.

Fractional anisotropy predicts N1 propagation velocity
Given that FA can be considered a proxy for white matter bundle integrity and myelination, we hypothesized that velocity would be related to the mean FA along a streamline trajectory. As a corollary, we further hypothesized that higher velocities would be mapped onto longrange, highly myelinated fiber bundles. Our analysis yielded a power law distribution of velocities across all patients and electrodes (Fig. 5). The mean FA along the streamline path strongly predicted the propagation velocity between the streamline end points (R 2 ¼ 0.85, β ¼ 7.73, p ¼ 7.14e-7; Fig. 5).

N1 propagation velocity distribution
Given the observed correlation between FA and N1 propagation velocity, we further investigated whether the distribution of propagation velocities can be further broken down into categories of pathways, in which association bundles, e.g. the arcuate fasciculus or inferior longitudinal fasciculus, would have higher velocities than short range connections and U-fibers. To explore the distribution of velocities across categories of white matter pathways, we first sorted streamlines into major association bundles and U-fibers. We identified substantial numbers of streamlines in the arcuate fasciculus, uncinate fasciculus, the combined superior longitudinal fasciculi I, II, and III (SLF), the combined inferior longitudinal fasciculus, middle longitudinal fasciculus, and inferior frontal-occipital fasciculus (ILF/MdLF/IFOF), the frontal aslant tract (FAT), and long-range corticocortical connections associated with the sensorimotor cortex. These pathways were visualized along with their N1 propagation velocity distributions (Fig. 6). Supporting our hypothesis, we found that while U-fibers had a low proportion of high velocity (e.g. !4 m/s) pathways (2.5%), association most association bundles had moderate to large proportions of high-velocity connections (arcuate fasciculus: 24.5%, uncinate fasciculus: 28.8%, SLF: 22.1%, ILF/ MdLF/IFOF: 16.4%, and sensorimotor pathways: 23.1%. Of the association bundles, only the FAT had a proportion of high velocity connections less than 10% (9.6%).

Dynamic tractography visualization
In order to visualize the cortical distribution of CCEP properties we mapped individual-level CCEP quantities onto tractography streamline bundles, (Fig. 7). As predicted by the observed relationships between streamline length, streamline FA, N1 latency, and N1 velocity, CCEP properties varied across the cortex, but were roughly clustered by bundle length. To visualize the directionality and spatial trajectory of specific cortical pathways, we created an animation rendering the propagation of N1 CCEPs along streamline pathways. The supplementary video S1 demonstrates all N1 trajectories for a stimulation site in the STG, while supplementary video S2 demonstrates all CCEP N1 trajectories inbound to the STG.

Table 2
Correlation coefficients between primary outcome measures, latency, amplitude, streamline length, and velocity and the five clinical covariates. All R < 0.12 (i.e., p-value < 0.05), therefore all five covariates were excluded from the subsequent regression models. Silverstein et al. NeuroImage 215 (2020) 116763 4. Discussion

Overview
This is the first multimodal imaging study to integrate CCEP amplitude and latency with DWI-based tractography estimates of white matter path length and fractional anisotropy in the same individuals. With this unique dataset, we were able to demonstrate that the underlying white matter structural network can predict the amplitude, latency, and propagation velocity of the CCEP N1 component. The streamline length interconnecting two electrode-based ROIs negatively predicted the N1 peak amplitude and positively predicted the latency. The N1 propagation velocity was positively predicted by the mean fractional anisotropy along the propagation streamline and N1 propagation velocity distributions tended to be higher in association bundles relative to U-fibers. These results are consistent with the hypothesis that CCEP N1s are primarily conducted by cortico-cortical white matter bundles and that higher myelination improves propagation velocity of axon bundles by decreasing the resistance. Our work replicates and extends the pioneering work in which CCEP N1 amplitude and latency were shown to correlate with streamline count (Conner et al., 2011). Similarly, our work is in keeping with a recent study in which N1 amplitude and latency were shown to correlate with streamline lengths drawn from an atlas (Trebaul et al., 2018). Here, we provide deeper evidence for such a relationship with a successful replication of the same analysis with CCEP and DWI data acquired from the same patients. These results will help ground further studies of integrated CCEP and DWI data and may help improve the utility of CCEPs in localization of functionally-important white matter pathways connecting eloquent cortices.

CCEP mechanisms and morphology
CCEP quantification methodology is a developing field and no clear gold standard yet exists. Similarly, the CCEP itself is a complex phenomenon; interpretation of CCEP data may depend on the particular quantitative approach elaborated in a given study. Early evoked responses, prior to 50 ms post-stimulus, such as the N1 (Conner et al., 2011;Keller et al., 2014bKeller et al., , 2014aMatsumoto et al., 2017;Nishida et al., 2017;Trebaul et al., 2018), P1 (Keller et al., 2014b;Matsumoto et al., 2004), high gamma (>70 Hz) power (Crowther et al., 2019;Matsumoto et al., 2017), and broadband power (Logothetis et al., 2010) are thought to reflect direct cortico-cortical connectivity. Late evoked responses, occurring 80-250 ms post-stimulus, such as the N2 (Matsumoto et al., 2004), may reflect subsequent cortical deactivation (Logothetis et al., 2010), whereas some do not rule out the possibility that N2 reflects cortical responses more associated with indirect connectivity, such as cortico-thalamo-cortical loops (Keller et al., 2014b) In this study we utilized an anatomically constrained tractography approach, therefore our streamlines were limited to cortico-cortical connections. We therefore limited our CCEP measures to maximize the involvement of direct measures. The N1 was chosen as a primary measure of cortico-cortical Fig. 4. Group level relationship between streamline length and CCEP measures. Streamline length weakly correlates with N1 peak latency (left) and has a strong negative correlation with N1 peak amplitude right. Pink line indicates best fit line, dotted lines indicate 95% confidence interval. Each point is the average streamline length within each equal-sized length bins. . The distribution takes a power law shape, often observed in biological systems. Mean FA along the streamline pathway predicts N1 velocity (right). Pink line indicates best fit, dotted lines indicate 95% confidence intervals. Each point on the regression is the average velocity value within equal-sized FA bins. Binning was designed to provide a model comparable to the latency and peak models above.
connectivity based on the robust evidence from both human (Yamao et al., 2014) and animal models utilizing both micro-and mesoscale electrodes (Logothetis et al., 2010) which suggest it likely reflects excitatory potentials mediated by white matter pathways. Although the influence of oligo-or polysynaptic connections cannot be completely ruled out, monosynaptic effects likely explain most of the variance in our data.

Factors which potentially influence N1 propagation velocity
Myelination is known to increase the conduction velocity of action potentials in both central and peripheral axons, and this is partially reflected as increases in FA (Heckel et al., 2015;Whitford et al., 2011). Evidence from animal models indicates that the N1 is a local field potential generated by the accumulation of post-synaptic excitatory potentials in neurons downstream from the stimulation site (Logothetis et al., 2010). Therefore, although FA is only a rough proxy for myelination and can be influenced by other non-myelin factors in the neuropil (Scholz et al., 2013), it is likely variations in myelination which best explain the relationship between FA and propagation velocity.
Interestingly, although mean FA predicts velocity, and longer white matter bundles tend to have higher values of FA, N1 propagation velocity was observed to have a somewhat heterogenous distribution across the brain (Fig. 6). This observation may suggest that there are 1) other factors which alter FA along the streamline path other than myelination (Scholz et al., 2013) and 2) additional macrostructural factors such as the number of axons in the bundle and the spacing of nodes of Ranvier could impact N1 propagation velocity, (Hartline and Colman, 2007). While we did not find correlations between velocity and the clinical variables noted above, this may be due to the relatively small sample size (n ¼ 23), limited data below the age of 10, and none above the age of 20. It is plausible to hypothesize that the size of developmental changes in effective or structural connectivity was too small to detect with the narrow age distribution available in this study. Further studies including larger number of young children are warranted to determine the developmental changes in dynamic tractography.
Furthermore, AEDs are reported to reduce the amplitude of motorand TMS-evoked potentials (Premoli et al., 2017). All patients in the present study were under Fosphenytoin during the CCEP data acquisition Lower plots indicate the velocity distribution for each association bundle as well as U-Fibers (far right). Histograms are plotted on a semilog axis. Association bundles have a greater proportion of high velocity pathways. In each plot, a dashed vertical line indicates the "high velocity" cut-off at 4 m/s. B.H. Silverstein et al. NeuroImage 215 (2020) 116763 procedure. Likewise, all patients were, up until intracranial electrode placement, were receiving at least one oral AED. Individuals not taking AEDs may exhibit different relationships between white matter and propagation velocity, but CCEPs cannot be measured from healthy individuals. Further work should explore the distribution of velocities across different white matter bundle types and lengths and incorporate more detailed quantification of macroscopic water diffusivity using NODDI (Winston, 2015;Zhang et al., 2012) and myelin water fraction (Wu et al., 2006).

Binned regression
We conducted a binned regression in order to improve the signal-tonoise ratio, particularly with respect to the tractography-based measures. A limitation of the binned approach is that it can reduce the variance in the data. Exploratory analysis of the unbinned data using a mixed model approach (including patient as a random effect) supported the binned peak/length and velocity/FA relationships. Both relationships, however, were better fit by a loglinear rather than a linear model (peak/length: linear adjusted R 2 ¼ 0.11, loglinear adjusted R 2 ¼ 0.16; velocity/FA: linear adjusted R 2 ¼ 0.12, loglinear adjusted R 2 ¼ 0.19). An unbinned mixed model failed to indicate a relationship between length and latency as observed in the binned approach (adjusted R 2 ¼ 0.03). However, this is not surprising given the weaker binned results and that N1 propagation velocity varies across connections. In the case of the relationship between length and latency, most of the variance comes, from the variance in velocity. The substantial variation in N1 propagation velocities across white matter bundles, explored in Figs. 5 and 7, likely explains the lack of a clear brain-wide relationship between length and N1 peak latency.

Relevance for behavior
CCEP networks have been shown to capture broader networks than resting state networks (Hebbink et al., 2018). This suggests that CCEP analysis reveals the full repertoire of a network, whereas resting state or task-based designs show how the brain may utilize a subset of the network to sustain a particular psychological process. Yet, although CCEPs do add an additional layer of information onto structural networks, neither CCEP networks nor DWI-based networks can provide direct information about the cognitive processes these networks support. This information can only come from study designs that specifically probe or interfere with a psychological process. Future work will integrate data from electrical stimulation mapping and task-induced high gamma data (Flinker et al., 2015;Hermes et al., 2015) with dynamic tractography to determine the precise relationship between cognition and observed structural networks.

Methodological considerations on iEEG sample rate and spectral responses
In this study, we focused our analysis on the N1 peak amplitude and latency as calculated from trial-averaged evoked potentials. In addition to evoked potentials, there is substantial evidence that high gamma power (Crone et al., 2001;Kambara et al., 2018bKambara et al., , 2018aNakai et al., 2017;Nishida et al., 2017) and broadband power (Bartoli et al., 2019;Hermes et al., 2015;Kucewicz et al., 2017) are complementary iEEG measures of cortical activation. Where evoked potentials (after averaging across trials) capture responses phase-locked to the stimulus, high gamma and broadband power reflect a summation of phase-locked and responses with a variable phase component (Crone et al., 2001;Flinker et al., 2015). Recently, high gamma (Crowther et al., 2019;Matsumoto et al., 2012) and broadband power (Logothetis et al., 2010) have also been demonstrated to reflect the strength of cortico-cortical connectivity as estimated by CCEPs. However, at a sample rate of 1 kHz, as was done in this study, the 0.3 ms electrical stimulation used for CCEPs generates a large degree of noise in the peri-stimulus window at and above approximately 60 Hz, (Sugiura et al., 2020), effectively obscuring biological activity in that frequency range. Sample rates upwards of 2 or 3 kHz may be necessary for effectively quantifying spectral power above 60 Hz during CCEP trials (Crowther et al., 2019). Future work will utilize these higher sample rates.

Methodological considerations for CCEP detection
No gold standard for CCEP detection currently exists. Generally, the procedure utilized in the field is to choose a measurement of CCEP "activity" such as peak amplitude, area under the curve, or spectral response magnitude, then compare it to a baseline (or other null) distribution, resulting in a Z-score. A Z-score threshold is selected, above which a CCEP is considered significant, or an FDR correction is employed. In this current study, we followed the results of Trebaul et al. (2018) which indicated a trade-off between sensitivity and specificity as the Z-score threshold is increased from 3 to 7, with 5 providing a reasonable mid-point. Testing alternative Z-score thresholds of 3 and 6 produced qualitatively similar relationships as those reported.

Generalizability
A common critique of research utilizing iEEG data is that it is not generalizable outside of the epilepsy population. While this is an important consideration, several factors allay concerns about a lack of generalizability. First, as with many iEEG studies, the sample in this study is heterogenous in terms of age, sex, etiology, seizure onset zone, and antiepileptic drugs. Therefore, it is likely that trends which emerge are Fig. 7. N1 peak amplitude (upper left), peak latency (upper right), propagation velocity (lower left) and mean FA (lower right) from all stimulation sites that resulted in significant CCEPs have been mapped onto streamlines for one patient. All figures are a lateral view of the right hemisphere. Longer streamlines, e.g. the arcuate fasciculus and longitudinal fasciculi tend to have higher velocities, but lower voltages and longer latencies than short range and U-fiber pathways. Data are from Patient 15 (Table 1). due to biological features common to all patients. In support of this theory, no clinical variables correlated with the outcome measures presented in this study. Furthermore, in this study all electrodes which showed interictal spiking activity or were colocalized with the seizure onset zone were excluded from analysis, thereby minimizing the influence of epileptogenicity on the results presented.

Individual differences
Our statistical approach combined all patients into a single analysis and therefore did not incorporate individual differences. In our view, this approach is justified due to the lack of correlations between individual covariates (age, sex, and number of antiepileptic drugs) and the outcome measures. Nevertheless, in neuroscience in general, and particularly in epilepsy, it is important to think one brain at a time. It is likely clinical variables will impact how individual brains deviate from the general trend observed here. Insofar as these results further clinical research, they should be adapted towards individualized analyses.

Surface versus depth electrodes
Our study made use of grids and strips of electrodes placed on the cortical surface. This technique provides excellent spatial resolution and broad coverage. However, the cortical surface electrodes are unable to sample deep cortical areas such as the cortex at the bottom of a sulcus. Likewise, we were unable to sample subcortical structures such as the thalamus, basal ganglia, or cerebellum. Therefore, these areas are not included in our analyses. However, our analysis combined pathways from diverse brain regions to highlight general relationships between structural and effective connectivity properties. A recent study utilizing stereotactic depth electrodes also found convergence between CCEP and DWI-based connectivity (Donos et al., 2016), thus, insofar as deep and subcortical brain structures are interconnected by similar white matter structures as those measured in our study, it is possible that the relationships observed here extend to those regions.

Conclusion
Evidence of the tight relationship between DWI-based tractography and the CCEP N1 may provide greater guidance in the use of DWI in surgical planning through visualization of the likely CCEP pathway. Furthermore, this study sheds additional light on the functional tractography connectome, which is of general interest to human neuroscience. Future research should further dissect the functional tractography networks to tease apart likely directional structures and local modulations of dynamic structures.

Data and code availability
All data and code used to generate the analyses reported in this manuscript are available upon request. All requests for data sharing will be reviewed by the corresponding author. Priority will be given to peerreviewed projects. If there is conflict or uncertainty regarding the merit of a request, a University Institutional Review Board representative will be consulted.

Ethics statement
The Institutional Review Board at Wayne State University has approved this study. We obtained informed consent from the guardians of patients and assent from children older than 13 years.

Funding
This work was supported by grants from the National Institutes of Health (R01 NS089659 to J.J. and R01 NS064033 to E.A.).

Declaration of competing interest
None of the authors have disclosures or conflicts of interest regarding this study.