Towards linking diffusion MRI based macro- and microstructure measures with cortico-cortical transmission in brain tumor patients

We aimed to link macro- and microstructure measures of brain white matter obtained from diffusion MRI with effective connectivity measures based on a propagation of cortico-cortical evoked potentials induced with intrasurgical direct electrical stimulation. For this, we compared streamline lengths and log-transformed ratios of streamlines computed from presurgical diffusion-weighted images, and the delays and amplitudes of N1 peaks recorded intrasurgically with electrocorticography electrodes in a pilot study of 9 brain tumor patients. Our results showed positive correlation between these two modalities in the vicinity of the stimulation sites (Pearson coefficient 0.54±0.13 for N1 delays, and 0.47±0.23 for N1 amplitudes), which could correspond to the neural propagation via U-fibers. In addition, we reached high sensitivities (0.78±0.07) and very high specificities (0.93±0.03) in a binary variant of our comparison. Finally, we used the structural connectivity measures to predict the effective connectivity using a multiple linear regression model, and showed a significant role of brain microstructure-related indices in this relation.


Introduction
Diffusion Magnetic Resonance Imaging (dMRI) provides measures of brain macro-and microstructure noninvasively, although their accuracy is still undetermined ( Jeurissen et al., 2019;Jones, 2010;Maier-Hein et al., 2017;Thomas et al., 2014 ). In particular, dMRI-based tractography is considered insufficient to guide brain tumor resection due to inaccuracies in localizing cortical terminations of fiber bundles ( Duffau, 2005;Duffau et al., 2008;Jbabdi and Johansen-Berg, 2011;Jeurissen et al., 2019;. From the clinical perspective, invasive electrophysiological mapping with Direct Electrical Stimulation (DES) ( Mandonnet et al., 2010;Ojemann, 1983;Penfield and Boldrey, 1937 ) remains the "gold standard " for probing connectivity ( Friston, 2011;Parker et al., 2018 ) when resecting a tumor located in eloquent brain areas ( Duffau, 2005;Keles and Berger, 2004 ), as it allows to optimize the extent of resection ( Hamer et al., 2012 ) and accommodate inter-patient variability ( Vigneau et al., 2006 ). However, a clinical use of DES is mostly based on empirical knowledge and hardly reveals any structural information about neural connections affected by the stimulation ( Duffau, 2005;Duffau et al., 2008 ). In this paper, we study the link between complemen- Fig. 1. A schematic shape of a cortico-cortical evoked potential, which typically consists of four consecutive voltage peaks named P1, N1, P2, and N2, where N are negative peaks while P are positive ones. An alternative naming convention is also used in the literature, according to which the positive peaks are enumerated from zero, hence P0 and P1 (in parentheses). The plot conforms to the commonly used convention of drawing ECoG recordings with the reversed y -axis.
we propose to engage both macro-and microstructure information into this joint study of brain connectivity.
From the technical perspective, CCEPs typically consist of four consecutive voltage peaks named P1, N1, P2, and N2 ( Matsumoto et al., 2017;Terada et al., 2012 ), where N are negative peaks while P are positive ones 1 ( Fig. 1 ). Nonetheless, hitherto studies of CCEPs mostly concentrate on monitoring N1, which is attributed to excitation of pyramidal cells ( Alarcón et al., 2012;Keller et al., 2014 ). Furthermore, N1 is usually more pronounced in the recorded signal than the other peaks ( Vincent et al., 2020;, which designates it as the most distinctive feature to study. In our work, we measured both delays and amplitudes of N1 in order to probe the propagation of CCEPs along WM tracts. Our goal was twofold. First, we aimed to validate our dMRI-based structural connectivity measures using the CCEP-based effective connectivity as reference. For this, we correlated the measures of fiber pathways, including aggregated streamline lengths and log-transformed counts of streamlines that connected the stimulation sites and the recording electrodes, with the delays and amplitudes of N1s. Earlier studies ( Conner et al., 2011;David et al., 2013;Keller et al., 2014;Silverstein et al., 2020 ) showed a linear relation between these quantities, so we fitted a linear regression model to predict the effective connectivity with our macrostructural measures.
As the second goal, we addressed the variability of neural conduction at the microstructure level, which depends on axon diamater, myelin sheath thickness, and axonal membrane properties ( Sen and Basser, 2005 ). For this, we incorporated into our regression model a set of microstructure indices that quantify the tissue composition. Whenever possible, we preferred the clinically feasible metrics derived from Diffusion Tensor Imaging ( Basser and Pierpaoli, 2011 ) (DTI) or Diffusional Kurtosis Imaging ( Jensen and Helpern, 2010 ) (DKI), e.g. axonal water fraction volume or tortuosity of the extra-axonal geometry ( Fieremans et al., 2011 ). For those requiring long multishell acquisition, we used Mean Apparent Propagator with Laplacian regularization ( Fick et al., 2016a ) (MAPL) signal representation.
Our work contributes to the study of the structural underpinnings of the CCEP-based effective connectivity. Despite being used in the clinical practice, little is known about the propagation of DES-induced evoked potentials along WM tracts ( Vincent et al., 2020 ). Hitherto research in this domain ( Conner et al., 2011;David et al., 2013;Keller et al., 2014;Matsumoto et al., 2017;Parker et al., 2018;Trebaul et al., 2018 ), performed on epileptic patients, led to discoveries in neuronal conduction ( David et al., 2013;Matsumoto et al., 2004;Trebaul et al., 2018 ) and pathological zone characterization ( Parker et al., 2018 ). However, a generalization of findings from patients with a single type of pathology to normal brain networks raised reasonable criticism ( Keller et al., 2014 ). Yamao et al. (2014) addressed this issue by repeating the same experiment on brain tumor patients, reaching similar results. Here, we adapted their methodology to common clinical conditions of a brain tumor surgery, which implied a use of low-current DES (2-5 mA instead of 10-15 mA) and small-sized stimulating electrodes (5 mm between poles rather than 10 mm). Thanks to these modifications, an integrated study of effective and structural brain connectivity, suggested by Silverstein et al. (2020) , would become available in a broader set of patients. More importantly, though, our endeavour to predict CCEPs with dMRI-based data aimed to explain the role of brain macro-and microstructure in the cortico-cortical transmission. Having this, we could apply our structural model to anticipate the clinically valuable effective connectivity measures prior to craniotomy or use this information when treating patients that need no surgical intervention.
We report here our pilot study of 9 patients. For each of them, we correlated presurgical dMRI-based structural connectivity measures, including streamline counts and lengths, with delays and amplitudes of N1. In addition, we considered binary variants of the above quantities to assess the structural connectivity thresholds for CCEPs. Finally, we used the macro-and microstructure measures of WM to predict the CCEP propagation using a multiple linear regression model. Our results show a significant role of microstructure indices in this relation.

Methods
This study, approved by the French national ethics committee and registered in the Clinical Trials database (NCT 03503110) ( NIH U.S. National Library of Medicitne, 2009 ), did not modify the usual surgical nor brain mapping protocols.
Our data set consisted of newly acquired anonymized presurgical dMRI and intrasurgical ECoG signal recordings, both of which we describe in detail later in this section. They are not publicly available due to restrictions imposed by the administering institution. The authors will share them by request from any qualified investigator after completion of a data sharing agreement.
Each patient signed an informed written consent to participate in our study. All of them received the same antiepileptic medication (Levetiracetam). Only Patient 6 had a history of previous anticancer treatment ( Table 1 ).

Surgical settings
All patients in this study were operated in the Asleep-Awake-Asleep protocol ( Szelényi et al., 2010 ). This created an opportunity for us to observe clinical responses to DES applied in the awake condition and to record brain electrophysiological activity resulting from DES under general anesthesia.

Data acquisition
Diffusion MRI Table 1 Summary of patients data including type, location, and volume of the pathology (II, III, IV -WHO grades; L, R -left/right brain hemisphere), epilepsy, antiepileptic drug administered, and the history of previous anticancer treatment.  Table 2 Summary of parameters used for the intrasurgical electrical stimulation and the simultaneous acquisition of the ECoG signal in 9 patients. The square waveform of biphasic current was composed of pulses having constant widths of 1 ms each. The stimulation frequency was fixed at 2 Hz in all cases except for Patient 1 having 5 Hz. The current intensities ranged from 2 to 5 mA. The ECoG montages comprised of one or two short strips with 4 electrodes each and/or one longer strip with 6 electrodes, totaling 8, 10, or 14 recording electrodes. The stimulation sites were located in the proximity of accessible electrodes, i.e. all electrodes except for those slipped under the dura matter. Their number is provided as sites num . The fixed sampling rate of 2 kHz was used for all acquisitions. We acquired presurgical diffusion MR images with 2 × 2 × 2 mm voxel size resolution. Our protocol comprised 4 shells having ∈ {400 , 800 , 1550 , 3100} [s/mm 2 ] with {6 , 13 , 29 , 51} directions, respectively, interleaved with 6 images without diffusion weighting ( = 0 ). The pulse width was = 0 . 0322 s and the diffusion time Δ = 0 . 0450 s. We used the above acquisition parameters since they maximized image reconstruction accuracy of the MAPL signal representation ( Fick et al., 2016a ) under the constraint of 25 min of imaging time on the clinical GE 3T-scanner. The MAPL helped us obtain numerous brain microstructure indices ( Özarslan et al., 2013 ) that we present in Section 2.6 .
Direct electrical stimulation A spatial distance between anode and cathode of a stimulation probe, and a current intensity of DES are two crucial parameters that determine a volume of brain tissue that is excited in our experiment ( Vincent et al., 2017 ). However, neither of them can be freely altered. The former parameter is fixed by the probe manufacturer. In our study, we used common-type probes with 5-mm separation between poles. The latter parameter, current intensity, must be adapted to each patient individually by the neurosurgeon ( Duffau et al., 2008 ), leaving very little room for modifications. In our experiment, the current intensities lied in the range of 2-5 mA ( Table 2 ).
Other DES parameters offer more flexibility. In particular, the clinically used stimulation frequencies are 50-60 Hz ( Vincent et al., 2017 ). Nonetheless, for the purpose of recording CCEPs, much smaller values must be used since a typical electrical response to DES lasts about 250 ms or more ( Conner et al., 2011;Keller et al., 2014;Matsumoto et al., 2004 ). We started with 5 Hz for Patient 1. Then, we decreased the frequency to 2 Hz for all the remaining patients, aiming to extend the observed time window from 200 to 500 ms ( Table 2 ). In each case, we applied biphasic square wave pulse of the width 2 ×1 ms.

Electrocorticography
We used sterile strips of ECoG subdural electrodes manufactured by DIXI Medical®. The electrodes had 4 mm contact diameter with 10 mm center-to-center distance between electrodes within each strip. Depending on the size of a bone flap in a given patient's skull and the accessible area of the cortex, we placed 2 or 3 strips. Their configuration comprised one or two short strips with 4 electrodes each and/or one longer strip with 6 electrodes, totaling 8, 10, or 14 recording electrodes ( Table 2 ). For every patient, we positioned the strips in alignment with the cortical terminations of Arcuate Fasciculus (AF) and Superior Longitudinal Fasciculus III (SLF3) obtained from tractography ( Fig. 2 a).
We recorded the ECoG signal right after tumor resection, while the patient was under general anesthesia. During the acquisition, the neurosurgeon stimulated consecutive cortical sites located in the proximity of all accessible electrodes. Note that parts of strips could be inserted beyond the bone flap, slipped between the dura matter and the cortex, allowing to record the signal in areas inaccessible to stimulation ( Fig. 2 b). We used a 32-channel signal amplifier manufactured by TM- Fig. 2. ECoG electrodes for recording CCEPs in Patient 7: (a) presurgical planning of the locations of electrodes (red circles) and stimulation sites (green crosses) using the tractography-based dissections of Arcuate Fasciculus and Superior Longitudinal Fasciculus III (both pictured jointly as blue streamlines), (b) intrasurgical positions of the ECoG electrodes and the stimulation sites (marked with labels S1-S10).
Si® with our recording electrodes set to the common reference mode, the ground electrode placed on the patient's back, and all the remaining inputs being short-circuited. Finally, the recordings were made using ASALAB software by ANT Neuro® with a 2 kHz sampling rate.

Structural connectivity
For each patient, we ran probabilistic tractography algorithm iFOD2 ( Tournier et al., 2010 ) twice, using the same dMRI data acquired prior to the surgery. The first run was destined for surgical planning and neuronavigation. To this end, we computed a whole-brain tractogram containing one million streamlines. From that, we dissected AF and SLF3, as mentioned earlier, using MI-Brain software by Imeka®. Next, we overlaid the streamlines of interest on a T1-weighted image that we transferred into the neuronavigation system produced by Brainlab®. In the operating room, the neurosurgeon placed the ECoG electrodes on the cortex in the precise locations of the WM bundle terminations indicated with tractography and localized spatially with the neuronavigation system ( Fig. 2 ).
Our second tractography was performed retrospectively, after a surgery. Here, we considered two types of Regions Of Interest (ROIs): seeding ROIs centered in the stimulation sites and target ROIs centered in the ECoG electrode locations. We obtained their spatial coordinates from the neuronavigation system that had been used during the surgery. Note that these coordinates referred to the points on the cortex, whereas tractography algorithms typically produce streamlines originating at the Fig. 3. Sample spherical ROIs surrounding two endpoints -a stimulation site (seeding ROI) and an ECoG electrode location (target ROI). The shortest streamline (showed as thick black line) found between them by a tractography algorithm begins and terminates near the boundaries of the spheres, which leads to underestimation of the streamline's length. As a compensation, the streamline's length is increased with the average distances to the centers of ROIs, as illustrated with dotted lines.
white/gray matter boundary ( Jeurissen et al., 2019 ). We accounted for this by defining our ROIs as spheres with 10-mm diameter, centered at the cortex ( Fig. 3 ), which is a common practice in such scenarios ( Reveley et al., 2015 ). Next, for each seeding ROI, we ran the DIPY implementation ( Garyfallidis et al., 2014 ) of the probabilistic tractography algorithm iFOD2 using 5000 random seeds per voxel. The step size equaled 0.5 mm with the maximum angle of 30 degrees. We set the maximum streamline length to 500 mm and defined the threshold on Generalized Fractional Anisotropy (GFA < 0.25) as a stopping criterion. From all the streamlines produced, we only retained those that touched any part of the target ROIs. As illustrated in Fig. 3 , the shortest streamlines thus obtained often originated and terminated near the boundaries of spherical ROIs, which led to underestimation of the streamline lengths. As a compensation, we increased the streamline lengths with the average distances to the centers of ROIs on both ends ( 2 × 3 . 7 mm in this study 2 ).
Having these, we computed the log probabilities of succeeding the Monte Carlo connectivity experiment. Technically, we divided the streamline counts by the total number of seeded tracts and applied log 10 transform to account for the streamlines distribution which is known to be biased toward shorter links ( Calabrese et al., 2015;van den Heuvel et al., 2015;Parker et al., 2018 ). Then, we took reciprocal values of our log-transformed streamline counts ( -log count ) to ensure the same monotonicity type as in the length-based macrostructural measures and thus obtained coherent presentation of results.

Effective connectivity
Our ECoG recordings were multi-channel time series with manually triggered time tags marking stimulation events. In post-processing, we epoched the data and zeroed out stimulation artifacts. Next, we subtracted the common average (i.e. an average of all recording electrodes at a given time) from the signal to reduce the noise induced by shortcircuited inputs and remove ECG artifacts. Finally, we averaged the epochs to increase signal to noise ratio.
A typical CCEP consists of four consecutive voltage peaks named P1, N1, P2, and N2 ( Matsumoto et al., 2017;Terada et al., 2012 ), where N are negative peaks while P are positive ones ( Fig. 1 ). In this work, we studied the propagation of N1s, which we were able to identify with the highest confidence among the three voltage peaks recorded in our ECoG signals. To this end, we implemented a Python script that recognized characteristic U-shaped patterns of N1s based on the derivative of the signal. In order to reduce noise, we applied a 30 Hz low-pass filter prior to identifying N1s. As a result, we were confident about the peaks that we found, although we might have potentially disregarded some of the less reliable ones due to their low amplitudes.

Correlation between structural and effective connectivity
After post-processing of our dMRI and ECoG data, we sought to evaluate a relation between the tractography-based macrostructural information and the effective connectivity measures. We approached this in two manners. First, by correlating streamline lengths and streamline counts ( -log count ) with N1 delays and N1 amplitudes. Second, by comparing binary variables which determined either a presence or an absence of structural and/or effective connection between endpoints. Such a comparison at the binary level helped us assess the structural thresholds above which CCEPs were observed and approximate the rate of false positives among selected tractography results, as we will explain later (in Section 4.2 ).
In the first approach, we defined the following three distance measures (given in millimeters) to quantify the streamline-based information: • minimum streamline length ( min str ) -a length of the shortest streamline connecting a given pair of endpoints, as in Silverstein et al. (2020) , plus the correction described in Section 2.4 ; • median streamline length ( med str ) -as above, using median instead of minimum; • distance along WM surface ( wm dist ) -a length of the shortest path connecting a given pair of endpoints traced along the surface of WM using Freesurfer ( Fischl, 2012 ).
In the second approach, we defined binary effective connectivity by labeling the pairs of endpoints between which we observed a propagation of CCEPs as effectively connected and all the remaining ones as effectively disconnected .
The case of binary structural connectivity was not so straightforward, since almost all pairs of ROIs in our data sets were joined with dMRIbased streamlines. In accordance with the previous studies, including Maier-Hein et al. (2017) , Thomas et al. (2014) , and de Reus and van den Heuvel (2013) , we assumed that a considerable number of links associated with lowest connectivities were false positives. Hence, we introduced a cut-off threshold on the macrostructural connectivity measures that would filter out potentially spurious links. All the pairs of endpoints for which the respective measure surpassed the threshold would be considered as structurally connected , whereas all those below the threshold would be considered as structurally disconnected . However, there is no common agreement on which cut-off values are most appropriate in such scenarios. Parker et al. (2018) used an arbitrary threshold for the whole study. Here, we adjusted the cut-off values per subject to account for the inter-patient variability and the differences in the DES parameters. To this end, we defined connectivity matrices = [ ] representing stimulation sites as rows and recording electrodes as columns. In the length-based connectivity matrices, we computed the coefficients by dividing min str distances between > 0 source and > 0 target ROIs by a sum of all such distances as follows (1) In the count-based one, we used the log-transformed streamline counts instead of min str . This way, we obtained two probability distributions, a length-based and a count-based one.
Next, we sought cut-off values ensuring best agreement between the binary structural connectivity maps and the binary effective connectivity maps. Following Parker et al. (2018) , we used Jaccard Index (JI) for this purpose. Let us recall that JI measures a ratio between an intersection and a union of two sets. When applied to binary maps, it computes their mutual similarity score. Hence, by maximizing JIs, we obtained cut-off values on the count-and length-based probabilities that maximized agreement between the macrostructural and effective connectivity measures.
Finally, we estimated the size of area affected by DES using macrostructural measures. For this, we quantified the streamlines that linked any two (binary) effectively connected regions.

Prediction of effective connectivity empowered with microstructure indices
Taking into account that propagation of evoked potentials is related to tissue microstructure (e.g. properties of nodes of Ranvier, Tasaki and Takeuchi, 1941 ; or axon diameter, Goldman and Albus, 1968;Rushton, 1951 ), we extended our set of macrostructural measures obtained from dMRI by quantities related to the brain WM microstructure, namely: • DTI metrics ( Basser and Pierpaoli, 2011  We computed the above metrics in all voxels containing one or more streamlines connecting our pairs of endpoints. Next, we ran the MRtrix implementation of the SIFT2 algorithm ( Smith et al., 2015 ) in order to assess a contribution of each such streamline and, with those in hand, we computed weighted averages of the microstructure indices.
Our approach in this regard was strictly data-driven. We applied stepwise regression method ( Efroymson, 1960 ) on the full set of indices for a feature selection. Also, we arranged the streamlines in ascending order with respect to their lengths and tested various subsets of streamlines restricted with low-pass and high-pass filters with cut-off values defined by percentiles of lengths , ∈ {0 , 10 , 20 , … , 100} . Note that, particularly, the case of = = 0 used no microstructure information at all. We will refer to it as macrostructure only and present for comparison.
Having relatively few data, some of which were noisy and incomplete (i.e. Patients 2 and 6), we trained our model with a leave-onepatient-out cross-validation using the data from Patients 1, 3, 4, 5, 7, 8, and 9.

Validity of the effective connectivity measurement
Our ECoG recordings were free from significant distortions other than stimulation artifacts. Despite the possible impact of volumeconducted potential ( Shimada et al., 2017 ), we were able to identify N1, P1, and N2 peaks in many recordings, whereas P0 was most often covered by a stimulation artifact. Our cortical responses preserved the variability in timing and strength of response at regions equidistant from the stimulation site, which is not likely to be due to volume conduction, as explained by Keller et al. (2014) .
All observed N1s occurred between 10 and 60 ms after the onset of DES, with most of them falling into the time interval 20-40 ms ( Fig. 4 a). Fig. 4. In most of the cases, CCEPs emerged 20-40 ms after the stimulation onset and propagated 10-60 mm away from the stimulation sites. The first two columns present the histograms of (a) delays and (b) amplitudes of N1 peaks. The remaining three columns show the histograms of our macrostructural measures: (c) minimum and (d) median lengths of dMRI-based streamlines connecting source and target ROIs, and (e) distances between the respective pairs of ROIs measured along the white matter surface. The histograms in blue present all structural connections, while the ones in orange refer to those among the structural connections for which we observed propagation of CCEPs (effective connectivity).
Note that an incomplete montage of electrodes had a visible impact on the characteristics of recorded N1s. In Patient 2, with 6 + 4 = 10 electrodes ( Table 2 ), the fraction of early N1s, occurring approximately 20-30 ms after DES, was over-represented with respect to all the other cases ( Fig. 4 a). Also, the range of amplitudes was much narrower for this patient, as well as for Patient 6 with 4 + 4 = 8 electrodes ( Fig. 4 b).
From the length-based perspective, the distances between stimulation sites and recording electrodes, measured either as shortest streamline lengths or WM surface distances, spread between 10 and 130 mm, whereas median streamline lengths lied in the range 30-130 mm ( Fig. 4 c-e). These ranges visibly changed after filtering out all the pairs of ROIs for which we did not observe effective connectivity. As a result, mostly the short-distance connections and a few long-distance ones were preserved, suggesting that most of our observed CCEPs traveled along U-fibers. Such modified min str and wm dist measures peaked at around 20 mm, while most of the med str values lied between 60 and 80 mm ( Fig. 4 c-e). We will consider only these filtered macrostructural measures in our correlation analysis below.
Finally, let us note that DES performed during the awake part of a surgery provoked clinical symptoms in all patients. In the representative example of Patient 1, we observed speech arrest, anomia and phonological paraphasias while applying DES in the designated stimulation Recorded propagation of CCEPs while applying DES to Patient 1. All plots illustrate the locations of the 6-electrode strip and the two parallel 4-electrode strips placed on the cortex. The schemes given at the top present the clinical symptoms, observed when applying DES, and the propagation of CCEPs from each of the stimulation sites (green crosses) shown as N1 delays (quantities are proportional to the diameter of the red circles) or tiny black boxes marking the electrodes where we observed no N1. The field at the bottom shows a sample denoised ECoG signal acquired during electrical stimulation as seen before (gray plots) and after averaging of event epochs (blue plots). The two zoomed plots show the N1 delays and N1 amplitudes recorded in the proximity of the stimulation (lower plot) and distally (upper plot). sites ( Fig. 5 ). The presence of such observable symptoms confirmed our choice of stimulation sites and DES parameters for inducing CCEPs.

Correlation of the connectivity measures
The individual patient features such as age, sex, pathology type, tumor volume, and presence or absence of epileptic seizures, turned out uncorrelated or very weakly correlated with the propagation of CCEPs ( Fig. 6 ). Hence, we disregarded these variables from our study.
We observed a positive correlation between the macrostructural and effective connectivity measures. The Pearson correlation coefficients for the count-based connectivity versus the N1 delays equaled on average 0.51 ± 0.17 ( Table 3 ). In the distance-based measures, the highest coefficients spanned between 0.34 and 0.71 ( Table 3 ). Note that in all cases except for Patient 2, the N1 delays were strongly positively correlated with the shortest streamline lengths. Even though the wm dist measure gave relatively high Pearson correlation coefficients in Pa-tients 1 and 3, the remaining scores of wm dist were lower, particularly in Patients 2, 4, and 6 (around or below 0.0). The medians of streamline lengths also showed positive correlation with the N1 delays, yet often the lowest among the three distance-based measures, with one notable exception of Patient 2.
We emphasized the role of uncontrollable or partially controllable factors in the experimental setup by providing (in Tables 3 and 4 ) the four parameters that inevitably influenced our outcomes: current intensity of DES, number of recording electrodes used, number of samples with identified N1 peaks, and average N1 delays/amplitudes with the respective standard deviations. Note that the reduced ECoG montages, composed of 10 or even 8 electrodes (Patients 2 and 6, respectively) delivered fewer data samples than the 14-electrode montages, i.e. 7-9 samples as opposed to 26-44 ( Table 3 ). They also resulted in the shortest observed average N1 delays (about 23-24 ms).

Table 3
The N1 delays were strongly positively correlated with the shortest streamline lengths ( min str ) in most of the cases. The three rightmost columns show the Pearson correlation coefficients computed for the N1 delays paired with the three studied distance measures. The respective p -values are given in parentheses. The highest scores per patient, spanning between 0.34 and 0.71, are printed in bold. Note that a positive correlation was also held for the reciprocal of log-transformed streamline counts ( -log count ). Additionally, in order to emphasize the role of uncontrollable factors in experimental setup, the following are specified: current intensity of DES, number of recording electrodes used, number of samples with N1 peaks identified, and average N1 delays with the respective standard deviations.  Table 4 The square roots of N1 amplitudes were strongly positively correlated with the inverse of the shortest streamline lengths ( min str −1 ) in most of the cases. The three rightmost columns show the Pearson correlation coefficients computed for the square roots of N1 peak amplitudes paired with the inverse of three studied distance measures. The respective p -values are given in parentheses. The highest scores per patient, spanning between 0.32 and 0.79, are printed in bold. A positive correlation was also held for the inverse of minus log-transformed streamline counts ( -log count −1 ) in all cases except for Patient 6. Additionally, in order to emphasize the role of uncontrollable factors in experimental setup, the following are specified: current intensity of DES, number of recording electrodes used, number of samples with N1 peaks identified, and average N1 amplitudes with the respective standard deviations. We made analogous observations from the study of correlations between the N1 amplitudes and the macrostructural connectivity measures. In fact, we correlated the square roots of N1 amplitudes with the inverse of minus log-transformed streamlines count -log count −1 and inverse distances, i.e. min str −1 , med str −1 , and wm dist −1 . We dis-covered empirically that the N1 amplitudes were rather inversely than directly proportional to the distances, and that the square root function better fitted the above proportion. The Pearson correlation coefficients in the case of N1 amplitudes were lower than the ones in N1 delays. They reached on average 0.37 ± 0.26 for the count-based measures and spanned between 0.32 and 0.79 in the best cases among distance-based measures ( Table 4 ). The inverse min str was strongly positively correlated with N1 amplitudes in all the cases except for Patients 3 and 6, which received the lowest current intensities (2.0 and 2.2 mA, respectively). The remaining two distance measures, i.e. med str and wm dist , most often gave lower Pearson correlation coefficients than min str . The role of stimulation parameters was also visible. The N1 amplitudes were larger and more dispersed when the current was relatively high (3.5 mA or above).

Binary connectivity measures
In this two-class approach, we considered our endpoints as either connected or disconnected. Aiming to find the upper bounds of accuracy between structural and effective connectivity, we sought cut-off thresholds on our structural measures that would maximize accordance with the propagation of CCEPs.
Among the spectrum of potential thresholds, we chose for further study the ones that gave maximum Jaccard index on their respective structural measures. This way, we obtained two-class classifiers of the structural information, reaching the highest agreement with the binary effective connectivity maps, as illustrated for Patient 9 in Fig. 7 .  Fig. 7. A comparison of the structural and effective brain connectivity measures at the binary level allowed to assess the structural thresholds above which CCEPs were observed and approximate the rate of false positives among selected tractography results. Column (a) shows Jaccard indices computed for all possible tresholds on the count-based (top row) and the length-based (bottom row) binary structural connectivity in Patient 9. The maximum (marked with a filled circle) is reached at 0.007, which means that a high-pass filter with the cut-off value 0.007 ensures highest agreement between the structural and effective connectivities. The associated ROC curves are presented in column (b), while column (c) shows the confusion maps. In the latter, dark green squares represent true positives, light green -true negatives, dark red -false positives, and light red -false negatives.

Table 5
Summary of quantitative results of the maximum Jaccard Index ( JI ) classifiers based on the log-transformed streamline counts. The thresholds ( Thresh ) on the normalized count-based connectivity measure range between 0.006 (Patient 3) and 0.017 (Patient 6). The column labeled max len presents maxima of the shortest lengths among endpoints after applying a threshold. The next three columns show the maximum Jaccard Index together with the sensitivity and specificity of the count-based connectivity classifiers. The last four columns hold the coefficients of the respective confusion matrices (i.e. TP -true positives, FP -false positives, TN -true negatives, FN -false negatives). Note that the inter-patient averages gave 6 ± 3% false positives and 4 ± 2% false negatives. The maximum Jaccard index ranged between 0.43 and 0.68 for the log-transformed streamline counts ( Table 5 ), and between 0.39 and 0.74 for the minimum streamline lengths ( Table 6 ). The corresponding sensitivities averaged over all patients equaled to 0.79 ± 0.08 in the case of -log count and 0.78 ± 0.07 for min str , while specificities were considerably higher reaching 0.92 ± 0.03 and 0.94 ± 0.03, respectively.
Considering the binary structural connectivity matrices referenced to the binary effective ones, the average rate of false positives in the countbased approach was equal to 6 ± 3%, while 4 ± 2% were false negatives ( Table 5 ). In the minimum streamline lengths connectivity ( Table 6 ), the inter-patient averages gave a bit lower percent of false positives (5 ± 2%) and the same percent of false negatives (4 ± 2%) as in the count-based case.
It is worth mentioning that the above false positive and false negative rates corresponded to relatively small brain regions located in the proximity of the stimulation sites. Note that the thresholds that we put on the structural connectivity measures restricted our sets of streamlines to those below 17-55 mm in the count-based approach ( Table 5 ), and 14-24 mm in the length-based approach ( Table 6 ). In other words, our structural measures reached such a good agreement with the CCEP propagation on the streamlines shorter than or equal to the above limits.
Prediction accuracy The first class of our linear regression models, which used macrostructure information only, produced similar mean residuals for each of the four input measures ( Table 7 ). The root mean squared errors (RMSEs) of N1 delays spanned between 8 ± 9 ms ( min str ) and Table 6 Summary of quantitative results of the maximum Jaccard Index ( JI ) classifiers based on the shortest streamline lengths. Here, we simply disregarded all the structural links which were longer than a certain limit length ( limit len ) determined by our threshold ( Thresh ) in a length-based connectivity matrix. The thresholds ranged between 0.008 and 0.024, which corresponded to the limit lengths between 14 and 24 mm. The confusion matrix coefficients (i.e. TP -true positives, FP -false positives, TN -true negatives, FN -false negatives) were comparable with the ones in the count-based binary connectivity. Although, the inter-patient averages gave a bit lower percent of false positives (5 ± 2%) and minimally higher percent of true negatives (76 ± 6%) than in the other measure.  Table 7 Microstructure information improved the accuracy of prediction of the effective connectivity measures, especially N1 delays. For each dependent variable (showed in the first column), i.e. N1 delay and N1 amplitude, a set of four prediction models was built using the respective macrostructural measures as base explanatory variables (second column). The mean R 2 scores and root mean squared errors (RMSE) obtained with the leave-one-patient-out cross-validation are given (with standard deviations) in the columns entitled macrostructure only . The same prediction models were extended by adding microstructure indices selected through stepwise linear regression. This modification resulted in a systematic increase of R 2 scores and a decrease of RMSEs, as presented in the columns entitled macro-and microstructure . In each case, the best results are printed in bold.  9 ± 11 ms ( med str ). Nonetheless, a dispersion of residuals was relatively high. Particularly, when predicting N1 amplitudes, the RMSEs varied from 240 ± 384 μV ( -log count ) to 267 ± 418 μV ( med str ). Variances of the effective connectivity data were best explained by min str in the case of N1 delays ( R 2 = 0 . 11 ± 0 . 27 ) and -log count for N1 amplitudes ( R 2 = 0 . 06 ± 0 . 14 ). The stepwise regression method produced consistent results regarding microstructure features selection ( Table 7 ). Tortuosity of the extraaxonal space was chosen in all the four models aimed at predicting N1 delays. Also, mean and radial kurtosis metrics, non-Gaussianity, and NODDI indices (ND and ODI) appeared repeatedly, while AWF was chosen only once to improve the prediction of N1 delays. Among these metrics: Tort, ND, ODI, and NG increased with the N1 delays, while MK and RK were negatively related.
Analogously, for predicting N1 amplitudes, the stepwise regression method selected FA for -log count and wm dist , while FA and NG ⟂ for min str and med str . Out of these two metrics: FA was positively, whereas NG ⟂ negatively related to the N1 amplitudes.
As we observed, the class of multiple linear regression models, based on a combination of macro-and microstructure data, reached higher prediction accuracy than the one using macrostructure only ( Table 7 ). However, their performance varied depending on the length of streamlines along which the microstructure indices were computed ( Fig. 8 a and  d). The average R 2 score of the model based on min str extended with the microstructure indices visibly increased and became less dispersed after filtering out the shorter half of streamlines ( Fig. 8 ). Interestingly, the R 2 rose much more for N1 delays (reaching 0.33 ± 0.09) than for N1 amplitudes (0.12 ± 0.19), which suggests that our microstructure metrics better captured the variability in conduction velocity than the signal decay.
Finally, let us point out that the macrostructure only models visibly overestimated the low values of N1 delays and N1 amplitudes, yet understimated the high ones ( Fig. 8 b and e). An inclusion of the microstructure indices helped reduce this bias to some extent ( Fig. 8 c and f).

Discussion
We found a large accordance between the structural and effective connectivity measures in brain tumor patients and showed a contribution of microstructure indices in explaining the cortico-cortical transmission. Particularly, the CCEPs that we recorded instrasurgically were linearly related to the indices quantifying axon dispersion and WM tissue composition, which we discuss further in this section.
Our study adds to the growing evidence that the propagation of CCEPs corresponds with the structure of WM bundles ( Conner et al., 2011;Keller et al., 2014;Matsumoto et al., 2017;Parker et al., 2018;Yamao et al., 2014 ). = 50 , = 100 (c and f). In each of the four plots, an overestimation of low values and an underestimation of high values is seen, however the cases with the microstructure information included partially compensated this bias.

Large accordance between structural and effective connectivity measures
We observed a positive correlation between the structural connectivity measures and the propagation of CCEPs quantified as delays and amplitudes of the N1 peaks. The Pearson correlation coefficients equaled on average 0.54 ± 0.13 (for N1 delays) and 0.47 ± 0.23 (for N1 amplitudes).
Our results were solid compared to the other approaches described in the literature. In a similar study conducted on epileptic patients, Parker et al. (2018) obtained values of between 0.05 and 0.21 when correlating the absolute amplitudes of CCEPs with the normalized streamline counts between endpoints in the frontal and/or temporal lobes. However, their areas of interests were considerably larger than ours due to the use of grids rather than strips of electrodes. In many of our patients, the size of brain cortex exposed during a surgery would not allow to place grids.
In an animal experiment, van den Heuvel et al. (2015) correlated the log-transformed streamline counts from in vivo dMRI scans of macaque brains with two tracer-based connectivity atlases, reaching an accordance of 0 . 25 ≤ ≤ 0 . 31 . In the mouse brains, Calabrese et al. (2015) compared the tractography-based structural connectivity with the tracer data on three levels of granularity. The authors obtained the Spearman correlation coefficients of 0.42 on finegrained ROIs, 0.71 on midle-sized ROIs, and 0.99 on coarse-grained ROIs ( ≪ 0 . 05 in all cases). However, the voxel-wise colocalization revealed rather weak correlation (0.23 ± 0.11) between neuronal tracer and probabilistic tractography. Also note that Spearman and Pearson correlation coefficients are defined differently and as such they cannot be compared directly.

Double role of the binary connectivity experiment
In our binary approach, we sought the highest agreement between the structural and effective connectivity measures in distinguishing connected from disconnected regions. This experiment allowed us to address two issues that we will discuss below:

Can binary effective connectivity help validate tractography?
Major criticism concerning tractography focuses on its tendency to produce false connections ( Jeurissen et al., 2019;Jones, 2010;Maier-Hein et al., 2017;Thomas et al., 2014 ). Even 36 ± 17% of streamlines generated by contemporary algorithms may link regions that are not actually connected ( Maier-Hein et al., 2017 ). Our binary connectivity experiment indicated 6 ± 3% of false positive and 4 ± 2% of false negative connections when comparing the log-transformed streamline counts with the propagation of CCEPs.
It is tempting to consider the CCEP-based effective connectivity as a method to minimize a fraction of false positive structural links. Note that effectively connected and effectively disconnected regions seem natural candidates as a reference measure of connectivity for tractography ( Trebaul et al., 2018 ). However, one must be aware of the flaws associated to tracing CCEPs. Our study showed that lower current intensities of DES (about 2.0-2.5 mA) came with lower N1 amplitudes and smaller numbers of effective connections found. This implied that the CCEP-based connectivity was more reliable when identifying false negatives than false positives. Indeed, a missing structural link between two effectively connected regions was very likely to exist (a false negative scenario), whereas a lack of effective connectivity between two regions joined with a streamline might not necessarily mean a false positive, yet also be caused by wrong DES parameters, low signal to noise ratio in ECoG recordings, or an interference with other electrophysiological activities. Jeurissen et al. (2019) claimed recently that tractography required filtering methods to decrease the rates of false positives. We believe that CCEP could deliver crucial priors in this regard, although the methodology related to recording and tracing the effective connectivity induced by DES still needs improvements.

Which structural connectivity is needed to ensure effective connectivity?
It is known that DES of the brain cortex primarily affects pyramidal cells at the stimulation site ( Keller et al., 2014 ). However, the propagation of CCEPs towards the target pyramidal cells through WM and its relation to the structural connectivity are still not sufficiently studied ( Keller et al., 2014;Matsumoto et al., 2017 ). For instance, it is not clear whether some measures of structural connectivity can discriminate between effectively connected and effectively disconnected regions. And if so, then where is the structural borderline between these two classes? Matsumoto et al. (2004) and Yamao et al. (2014) observed longranged effective connectivity between anterior and posterior perisylvian language areas. On the contrary, the propagation of CCEPs in our experiments was mostly local. It is worth noticing that in their experiment the current intensities ranged from 10 to 15 mA and the stimulating electrode had a 10-mm distance between poles, as opposed to 2-5 mA and 5-mm distance in our case. According to Vincent et al. (2017) , these two variables mostly affect the scope of DES. We thus believe that it is more appropriate to determine structural boundaries of the effective connectivity given certain DES parameters.
In our binary connectivity experiments, highest agreement between the structural and effective measures was reached when the streamline count-based probability surpassed 0.009 ± 0.004 or when the shortest streamline lengths remained below 20 ± 3 mm, which roughly approximated the limits of effective connectivity in terms of the structural measures. Additionally, this suggests that our CCEPs might have traveled along U-fibers rather than long fascicles like AF, since the typical lengths of the former lie between 3 and 30 mm ( Schüz and Miller, 2002 ).

Prediction of CCEP propagation with structural connectivity measures
The type and strength of correlation between the studied structural and effective connectivity measures suggested that we might predict the CCEP propagation by applying linear regression on the structure-based variables. Indeed, the residuals of the proposed prediction models were relatively low on average, although somewhat dispersed.
Our experiments showed that the linear regression of the logtransformed streamline counts or the minimum streamline lengths reached higher accuracy than the one based on med str or wm dist . It is then very likely that our CCEPs traveled along the shortest structural links between the endpoints rather than long-range fibers. Interestingly, the fact that min str outperformed wm dist in this comparison suggests that the trajectory from a stimulation site to a recording electrode as measured on the WM surface was not accurate enough to capture the CCEP propagation pattern. This supports the argument of communication via U-fibers, advocated by Conner et al. (2011) , rather than the impact of volume-conducted potential signalized recently in the literature ( Shimada et al., 2017 ).
In this work, we hypothesized that the propagation of CCEPs, measured as a spatial accumulation of N1 delays and a decay of N1 amplitudes, may implicitly relate to the properties of neuronal tissue connecting the source and target pyramidal cells. The observed increase of the R 2 score and decrease of its dispersion after including the microstructure-related indices in our linear regression models imply tangibly that these indices contained information related to CCEPs.
As a matter of fact, neural conduction largely depends on axon diameter, myelin sheath thickness and axonal membrane properties ( Sen and Basser, 2005 ), however an approximation of these microstructural features is challenging in vivo ( Assaf et al., 2008;Basser and Pierpaoli, 2011;Fick et al., 2016a;Özarslan and Basser, 2008;Veraart et al., 2020 ). One solution is to capture them through non-specific markers that quantify the tissue composition ( Pierpaoli et al., 1996 ). For instance, FA was already shown to correlate with conduction velocity ( Silverstein et al., 2020 ) and functional connectivity ( Boorman et al., 2007 ). In this work, we considered FA and other DTI-based indices, i.e. MD, RD, and AD, which often serve as surrogate measures of tissue microstructure ( Basser and Pierpaoli, 2011;Zhang et al., 2012 ). In addition, we included the metrics of DKI-based quantification of non-Gaussian diffusion, namely MK, RK, and AK, whose clinical significance has been shown in relation to the age-related demyelination ( Falangola et al., 2008 ) and brain tumor staging ( Raab et al., 2010 ). Simultaneously, Sen and Basser (2005) argued that tortuosity of the extra-axonal space, which also can be quantified with DKI ( Fieremans et al., 2011 ), limits the conduction in a restrictive geometry of fibers.
Numerous microstructure features with known impact on conduction velocity provide more specific information, although require longer acqusition than DTI or DKI. Jespersen et al. (2010) showed that neurite density estimates correlated with the myelination level, which led us to include the NODDI-based indices of ND and ODI derived from our multi-shell acquisition. The use of MAPL signal representation allowed us to incorporate a group of less common indices of tissue composition like QIV, considered more sensible than FA ( Fick et al., 2016b ), or NG. Surprisingly, RTOP, RTAP, and RTPP, belonging to the group of indices typically atributed to volume, cross-sectional area, and length of pores (respectively) ( Fick et al., 2016a;Özarslan et al., 2013 ) turned out less significant than other indices in our study.
On the other hand, the presence of FA and tortuosity among the features selected with our stepwise regression is by far the most intuitive. The impact of both these metrics is consistent with the previous reports ( Sen and Basser, 2005;Silverstein et al., 2020 ), similarly to ND and ODI, as we mentioned earlier. Interestingly, the N1 delays, reflecting de facto the conduction velocity, appeared visibly more related to the studied microstructure indices than the N1 amplitudes. The latter apparently decayed proportionally to the streamline lengths and counts, although maintained the variability that we were unable to explain with the WM microstructure. Alternative approaches to weighting the q-space indices or filtering tracts could probably improve those results.

Conclusions
In this work, we studied the relation between the structural connectivity measures obtained from dMRI and the effective connectivity measures based on the propagation of CCEPs in the brain tumor patients. Our experiments showed positive correlation between streamline lengths and counts with the delays and amplitudes of N1 peaks. When comparing the binary variants of the structural and effective brain connectivity, we discussed the potential of CCEPs to filter out the results of tractography and estimated the structural criteria of the CCEP propagation. Finally, we showed that brain tissue microstructure features help explain the propagation of CCEPs. The N1 delays and N1 amplitudes measured instrasurgically were linearly related to the indices quantifying axon dispersion and WM tissue composition. Particularly, the tortuosity of the extra-axonal space, mean and radial kurtosis, and neurite densitity indices considerably improved the prediction of effective connectivity measures. We believe that our findings extend the clinical significance of microstructure indices and contribute to the goal of understanding the propagation of CCEPs.