Causal Cortical and Thalamic Connections in the Human Brain

Causal Cortical and Thalamic Connections in the Human Brain | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Article Causal Cortical and Thalamic Connections in the Human Brain Josef Parvizi, Dian Lyu, James Stieger, Zoe Lusk, Vivek Buch This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-4366486/v1 This work is licensed under a CC BY 4.0 License Status: Published Journal Publication published 15 Jul, 2025 Read the published version in Nature Neuroscience → Version 1 posted You are reading this latest preprint version Abstract The brain's functional architecture is intricately shaped by causal connections between its cortical and subcortical structures. Here, we studied 27 participants with 4864 electrodes implanted across the anterior, mediodorsal, and pulvinar thalamic regions, and the cortex. Using data from electrical stimulation procedures and a data-driven approach informed by neurophysiological standards, we dissociated three unique spectral patterns generated by the perturbation of a given brain area. Among these, a novel waveform emerged, marked by delayed-onset slow oscillations in both ipsilateral and contralateral cortices following thalamic stimulations, suggesting a mechanism by which a thalamic site can influence bilateral cortical activity. Moreover, cortical stimulations evoked earlier signals in the thalamus than in other connected cortical areas suggesting that the thalamus receives a copy of signals before they are exchanged across the cortex. Our causal connectivity data can be used to inform biologically-inspired computational models of the functional architecture of the brain. Biological sciences/Neuroscience/Computational neuroscience/Dynamical systems Biological sciences/Neuroscience/Diseases of the nervous system/Epilepsy Biological sciences/Neuroscience/Cognitive neuroscience/Consciousness Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Introduction The brain’s dynamics are shaped by electrophysiological interactions throughout its subregions, encompassing not only cortical but also subcortical regions. Advancements in neuroimaging have provided significant insights into the global architecture of the brain’s functional connectivity 1 . However, we know considerably less about the dynamic, fast-paced and causal electrophysiological relationships in the human brain. This knowledge gap is even more pronounced when it comes to understanding the role of subcortical areas such as the thalamus, which is known to play a key role in modulating the global dynamics of the brain 2 . A classic method of studying causal electrophysiological relationships in the brain is by sending repeated single electrical pulses to intracranially implanted electrodes while recording the presence or absence of electrophysiological changes in all other areas of the brain where recording electrodes are present 3 . This method has been traditionally referred to as the study of cortico-cortical evoked potentials 4 (CCEPs) since it has primarily been used to study causal electrophysiological connectivity between pairs of cortical regions. Similar studies of thalamocortical and corticothalamic connections have been rarely conducted since direct recordings from and stimulations in the human thalamus are extremely rare in human neuroscience research 5–7 . As a result, it remains unknown if stimulation of the thalamus will have a different, or the same, effect on other regions of the brain. Moreover, with a few exceptions 7,8 , previous studies have been largely reliant on simple univariate measures detecting large peaks or the time-to-peak in the evoked signals, and limit responses to fixed windows of interest. These traditional methods are not able to capture the complex dynamics of physiological responses generated by electrical stimulations (Fig. 1 ); and they can vary significantly depending on the chosen cutoff z-scores for determining an effect, current intensity and the distance between stimulating and recording electrode contacts 9 . These problems become more significant when we lack prior data to inform our hypotheses about evoked potentials caused by human subcortical structures such as the thalamus. In the current study, we aim to use a novel approach to provide a map of causal electrophysiological connections in the human brain including the thalamus. Results Our results are based on intracranial recordings and stimulations in twenty-seven participants (40.7% female; mean age ± SD: 34.9 ± 10.0 years) diagnosed with focal epilepsy. Each participant provided data from an average of 180 (± 46) electrode sites. The entire dataset encompasses 4864 sites across all subjects, with each adjacent pair stimulated approximately 45 times (with sufficient inter-trial intervals to avoid overlapping effects), while recordings were taken from all other available sites. We leveraged a recently developed clinical method that involves The brain's functional architecture is intricately shaped by causal connections between its cortical and subcortical structures. Here, we studied 27 participants with 4864 electrodes implanted across the anterior, mediodorsal, and pulvinar thalamic regions, and the cortex. Using data from electrical stimulation procedures and a data-driven approach informed by neurophysiological standards, we dissociated three unique spectral patterns generated by the perturbation of a given brain area. Among these, a novel waveform emerged, marked by delayed-onset slow oscillations in both ipsilateral and contralateral cortices following thalamic stimulations, suggesting a mechanism by which a thalamic site can influence bilateral cortical activity. Moreover, cortical stimulations evoked earlier signals in the thalamus than in other connected cortical areas suggesting that the thalamus receives a copy of signals before they are exchanged across the cortex. Our causal connectivity data can be used to inform biologically-inspired computational models of the functional architecture of the brain. thalamic recordings in neurosurgical patients undergoing invasive stereotactic electroencephalography (sEEG) 5 . In this procedure, intracranial electrodes were implanted in three different sites of the thalamus — corresponding to anterior, mediodorsal, and pulvinar nuclei — by extending electrodes initially placed in the cortical areas. This procedure was designed for extended clinical indication without increasing the number of implanted electrodes. Identifying Common Stable Patterns in Stimulation Evoked Potentials To determine true physiological effects upon electrical stimulation, we used a non-linear manifold learning algorithm 10 integrated with human judgement, based on the spectral information of evoked potentials (EPs) across time and frequencies (Fig. 2 a-c). In this approach, we first labeled the data partially and tentatively, i.e., subject to subsequent data-driven corrections, based on two measures of connectivity: (i) trial-averaged power of EPs as a measure of strength, and (ii) inter-trial phase coherence (ITPC) across repeated stimulations as a measure of consistency for the evoked responses 11 . We applied a semi-supervised machine learning approach to decipher the inherent data structure from this preliminary labelling, and re-classified the data, within each subject, into two separate clusters ( Fig. S2 ): one cluster of data where stimulation of the seed site caused significant EPs in the target site, and another cluster of data where stimulation of the seed site did not cause any significant EPs in the target site. For simplicity, we called these two clusters “ activated ” and “ non-activated” , respectively. As a result, EPs in the activated cluster were differentiated from the spontaneous activity based on prominent changes in both power and ITPC relative to the baseline (pre-stimulation intrinsic neural activities), and also differentiated from the unnatural signals that we have manually labelled as noise (e.g., stimulation artifact or bad channels) by providing those typical cases. Once all subjects’ EPs were successfully classified as “activated” vs. “non-activated”, we ensured the cluster’s meaningful representation of the data by manually inspecting the original and reconstructed data for every subject. With the machine-learning refined labelling (activated vs. non-activated) as a benchmark, the preliminary activation criteria were re-evaluated for further sanity checks ( Table S1 ). We next aimed to characterize the prominent electrophysiological properties of the EPs in the “activated” cluster. Dimensionality reduction with uniform manifold approximation and projection (UMAP) revealed that the inherent structure of the data was different for stimulation-recording pairs within the same hemisphere (i.e., ipsilateral) compared to those across the two hemispheres (i.e., contralateral connections). Surprisingly, we also discovered that the inherent structure of the data was different for EPs generated by stimulation of the thalamic sites vs. cortical sites (Fig. 2 d; more label mapping can be found in Fig. S3 ). The thalamic vs. cortical recording sites, however, do not exhibit such obvious differences ( Fig. S3c vs. Fig. S3d ). We refer to this anatomical features (1. Stimulation from the cortex or thalamus, 2. Ipsilateral or contralateral connections) as UMAP localizers. We then searched for the significant features that set the data apart in the embedding space. To dissociate the data that contain the features of interest, we applied supervised learning with the UMAP localizers as labels. As expected, it resulted in four clusters in the embedding space with distinct spectral patterns, i.e., clusters of thalamic, cortical, ipsilateral, and contralateral EPs ( Fig. 2 e ) . Dissociating Distinct Neural Features in Evoked Responses Once the four clusters were identified, we used the power and ITPC spectrograms from each cluster for further statistical testing. Cluster-based permutation testing identified the significant time/frequency boundaries of three distinct neural features in the data that we labeled as Feature 1 (F1), Feature 2 (F2), and Feature 3 (F3) as tentatively shown in Fig. 2 e by visual identification. Details of statistical testing procedures are presented in the Methods . We summarize the three features as follows (Fig. 3 a, b): F1 is characterized by significant increases in power of high frequency activity in the gamma range and the strength of ITPC within the first 10-60ms. Of note, F1 was clearly dissociable from the artifact which was showed as a high-frequency (> 100 Hz) peak power immediately after the stimulation (< 10ms) and lacking signal continuation in both frequency and time dimensions ( Fig. S5) . F2 is characterizable by an increased power in the theta to alpha range, without phase locking, happening ~ 120 ms upon stimulation. Later than F2 (~ 200 ms post-stimulation), F3 is featured by a prolonged power increase in the theta band, with a phase-locking feature which indicates oscillation (Fig. 4 a). Within-subject comparisons showed that the ipsilateral cortical stimulations compared to contralateral cortical stimulations (i.e., COR-ipsi > COR-contr) caused significantly stronger F1 and F2 signals (Fig. 3 b). Similarly, ipsilateral thalamic stimulations compared to contralateral thalamic stimulations (i.e., THAL-ipsi > THAL-contr) caused stronger F1 and F2 signals (Fig. 3 b) i.e., stimulations of the thalamic and cortical sites lead to a stronger F1 and F2 signals in ipsilateral (compared to contralateral) brain sites. Of note, F1 was lacking in the contralateral recordings to the thalamic stimulation (i.e., THAL-contr). Importantly, the EPs generated from thalamic stimulations were distinct from the cortical stimulations by the presence of a strong F3 signal (i.e., higher theta frequency change starting emerging around 165 ms and persisting for about 250 ms (Fig. 3 b). By checking the signals in the time domain, combining the evidence of the evoked power and ITPC, we confirmed (i) the presence of oscillatory activity embedded in the F3 signal, and (ii) that the F3 signal is not a delayed F2 or continuation of it. In the time domain, F3 indicates a ~ 5 Hz oscillation, but F2 indicates a wider curve with jittered peaks across sites. In the time-frequency domain, both F2 and F3 share similar power spectrum, but their ITPC pattern differs significantly in that F2 does not have a phase-locking feature (Fig. 4 a). Moreover, unlike F2, F3 is lacking in cortical EPs, instead, it is specific to thalamic EPs with both THAL-ipsi and THAL-contr clusters showing this delayed and long-lasting F3 activity. We then investigated the spectrograms separately for different recording sites. Profiles of THALàCOR, CORàTHAL and CORàCOR in both ipsilateral and contralateral pairs reflected the main results of the UMAP localizers, which showed that the stimulation sites, rather than the recording sites, dominate the spectral features. Notably, multi-site thalamic coverage enabled us to examine intra- and inter-thalamic causal connectivity at the individual level (Fig. 3 c ) . The results provided unprecedent data from the thalamus where we could compare the connectivity of different thalamic subregions with the available brain recording sites and examine the effect of electrical stimulation of a given thalamus upon the activity of sites in other thalamic subregions within the same or across the other hemisphere ( Fig. S7 ). Timing and Global Patterns of Connectivity Indicated by Feature Presence Next, we aimed to understand how each of the three neural features unfolds in time within individual instances of connectivity for the four categories of interest (i.e., COR-ipsi, COR-contr, THAL-ipsi and THAL-contr). The group-averaged spectral information within the significant time windows from the COR-ipsi, COR-contr and THAL-contr categories was used to create the feature template for F1, F2 and F3, respectively (Fig. 4 a). Using a sliding-window correlation approach (see details in Methods and illustration in Fig. 4 b), we estimated the strength of feature representations (r*) and their presence in time. Corresponding to the previous group-level cluster-based permutation testing on the spectrograms, the decoding results also showed that both F1 and F2 peaks were seen more strongly for ipsilateral than contralateral EPs (Dr* = 0.11, t = 48.05, p FDR-corr < < 0.01, n-connection [conn.] = 57881, n-subject [sbj.] = 26 for F1, and Dr* = 0.09, t = 44.90, p FDR-corr < < 0.01, n-conn. = 57881, n-sbj. = 26 for F2, using the stimulation category “COR/THAL” as a covariate, corrected for the nested individual, regional and site effects in a hierarchical linear model [HLM]). Additionally, the decoding method was able to provide more precise time information than the previous test, which showed that both F1 and F2 were detected earlier in the EPs responding to cortical versus thalamic stimulation (Dlatency = 5.60 ms, t = 27.20, p FDR-corr < 0.4, latency in the 10–100 ms range], and Dlatency = 8.18 ms, t = 22.49, p FDR-corr < 0.4, latency in the 70–200 ms range]; using the cross-hemispheric category “ipsi/contr” as a covariate in an HLM). The decoding results also echoed the previous finding that F3 was seen mainly in evoked potentials pertaining to the thalamic stimulations, with F3 peaking higher in ipsilateral than in contralateral EPs (Dr* = 0.11, t = 48.05, p FDR-corr = 0.02, n-conn. = 6497, n-sbj. = 25), but without significant difference in the latency between ipsilateral and contralateral EPs (Dr* = 5.03 ms, t = 1.86, p = 0.06, n-conn. = 3106, n-sbj. = 25 for robust F3 with r* > 0.5, latency > 200 ms). Additional post-hoc comparisons among sub-groups are presented in Fig. 5 and SI Table S9 . The correspondence of the decoding results with the previous group-level spectrogram analysis suggests that our devised feature indices can reliably indicate the presence of the multivariate features with single values. We further examined the whole-brain spatiotemporal patterns indicated by the three neural features. To do so, we constructed connectivity matrix with each cell filled with the connectivity index of R for indicating the signal strength between two brain regions (R = \(\stackrel{-}{{r}^{*}}\) , averaged across the sites where the feature is present in the given regions). By organizing brain areas in the order of anatomical proximity, interesting patterns emerged: (i) the F1 matrix demonstrated a profile of modular connectivity, i.e., higher F1 representation was seen for connections between pairs located within the same hemisphere or lobe and within the ipsilateral thalamus (Fig. 6 a, b); (ii) the F2 matrix revealed wide-spread representations across regions and hemispheres and demonstrated that the global architecture of connections based on F2 did not differ much from the F1 matrix. This seems to suggest that F2 is a combination of jittered rebounds following the same paths for the F1 signals ( Fig. S10 ); (iii) unlike the F1 and F2 matrices which have a modular architecture, the F3 was sparsely present among the cortical connections, but largely present in the cortical responses (i.e., widespread and bilateral) associated with thalamic stimulation (Fig. 6 f, g). In a closer examination, we found that F1 corticothalamic and F3 thalamocortical connectivity features had the following unique characteristics in common: (i) F1 responses recorded in the thalamus had widespread and bilateral origins across the brain (i.e., stimulation of many cortical regions bilaterally could evoke responses in the thalamus) (Fig. 6 a), and vice versa for F3, the stimulation of a thalamic site would cause widespread cortical areas bilaterally to show delayed slow oscillations (Fig. 6 b); and (ii) both corticothalamic F1 and thalamocortical F3 responses were stronger and earlier than corticocortical F1 and corticocortical F3 responses, respectively. The statistical details are presented as below. We compared the cortical versus thalamic targets responding to the same cortical seed of stimulation in the same participant, and confirmed the F1 responses in the thalamus were earlier and stronger (Dlatency = 2.79 ms, t = 12.74, p FDR-corr < < 0.01, n-conn. = 48446, n-sbj. = 26; Dr* = 0.06, t = 18.20, p FDR-corr < < 0.01, n-conn. = 30965, n-sbj. = 26, with HLM corrected for the covariance effect of cross-hemispheric connections and the nested random effect of 1378 different stimulation sites of different brain areas in all subjects) (Fig. 6 a for R matrix, SI Fig. S10 for latency matrix). In other words, stimulation of a given cortical site changed the activity of the thalamus before it changed the activity of its cortical targets (Fig. 5 a and SI Table S9 for post-hoc latency comparisons). Moreover, the F3 matrix revealed that the stimulation of thalamic sites caused delayed low-frequency oscillatory effects (i.e., F3) in a large mantle of the cortex in not only ipsilateral but also contralateral hemispheres (Fig. 6 b). By contrast, the F3 feature was only sparsely present with the stimulation of non-thalamic seeds. Importantly, the latency of F3 responses caused by thalamic stimulations was significantly earlier than the non-thalamic seeds (36.27 ms earlier, t = 21.60, p FDR-corr < < 0.01, n-conn. = 20823, corrected for the covariance of cross-hemispheric effect, and the nested random effect of 3041 recording sites of different areas in all subjects with an HLM) (Fig. 5 c). Differential connectivity among thalamic subregions Up to this point, we have remained agnostic to the specific anatomical location of the thalamic sites where evoked responses were recorded or stimulations were seeded. However, each subject had more than one electrode in the thalamus targeting the anterior, mediodorsal and pulvinar nuceli 5 . Given that the precise anatomical boundaries of each thalamic nucleus within each subject may be variable and difficult to ascertain with neuroimaging data, and that the field of evoked responses or electrical stimulations may cross the nuclear boundaries, we refrained from labeling our thalamic sites according to specific nuclei, and instead, referred to them as anterior (antTH), middle (midTH), and posterior thalamus (pstTH). Leveraging the anatomical information of recordings and simulation sites - at the individual level - we compared the connectivity profiles of the three different thalamic sites. This analysis revealed unique site-specific profiles of connectivity between the thalamus and the cerebral cortex (Fig. 6 c, d, e ) . As describing the details of differential connectivity of thalamic sites are outside the scope of this report, we only highlight some of the most salient findings. The full results of pair-wise comparisons among the thalamic subregions for their bi-directional connections with each of the brain regions are provided in the supplementary Table S2-7 . Next, we compared the delayed thalamocortical F3 signal representation map with the ones for early thalamocortical F1 and early corticothalamic F1 maps, with pair-wise comparisons within the same pair of thalamic vs. non-thalamic sites. Using the measure of normalized feature representation \(\stackrel{-}{zscore\left({r}^{*}\right)}\) , as shown in Fig. 7 and detailed in Tables S6 and S7 , between a subregion of the thalamus and a given brain region of interest, we found significant asymmetry between thalamocortical F3 representations and thalamocortical F1 representations. Many regions had stronger F3 than F1 representations suggesting that the delayed thalamic (F3) outflow is represented across a more widespread cortical areas compared to the fast thalamic (F1) outflow signal (see brain regions with more brown colors in Fig. 7 a and positive estimate numbers in Table S6 ). The same applies to the comparison between thalamocortical F3 and corticothalamic F1 representations. In this comparison, the hippocampus (HPC) was an exception among all other examined brain regions: representation of F1 signals evoked by HPC and recorded in anterior and posterior thalamic sites was stronger than the representation of F3 signals recorded in HPC and evoked by the thalamus, with pair-wise comparisons for the same sites in the regions of interest (Dz = 0.30, t = 2.51, p FDR-corr =0.04, n-conn. = 190, n-sbj. = 21 for antTH; Dz = 0.62, t = 4.72, p FDR-corr < < 0.01, n-conn. = 194, n-sbj. = 18 for pstTH) ( Table S7 ). Discussion In our study, we addressed the limitations of past studies by leveraging a novel approach to characterize the complex features of electrophysiological responses that are evoked by not only cortical but also thalamic stimulations in the human brain. This is important because, as we documented in Fig. 1 , reliance on the magnitude of recorded neural responses or expecting typical N1 and N2 profiles in the evoked responses may greatly bias our views of causal connections in the brain. The magnitude of evoked responses can be affected by several confounding factors such as distance from the source of stimulation, volume conduction, and electrode impedance 12 and the time-to-peak may be elusive when there is no clear N1 or N2 peaks. While typical N1 and N2 waveforms reported in the CCEP literature are abundant, it is important to note that the past reports were based on group level analysis of data from cortical stimulations 13,14 . As noted in Fig. 1 , individual evoked signals may be far more complex than simple N1 and N2 waveforms and the magnitude of evoked responses may not always pass an arbitrary threshold we set for significance. This issue is even more pressing when waveforms from subcortical structures such as the thalamus are included. Therefore, in line with recent attempts by other investigators 7,8 , we felt it was imperative to employ a multivariate and unbiased approach for sorting the distinct profiles of connectivity across cortical and thalamic sites. In our encoding method, we relied on the machine learning algorithms of UMAP to identify the inherent variabilities among the datasets. We emphasize that our data-driven approach was not blind to biologically plausible features. Instead, it was directly guided by the 2/3 of the data that had already been labelled based on established neurophysiological criteria for significant evoked responses (i.e., significant amplitude above the noise level, significant inter-trial phase coherence concerning reliability of responses across trials, and more importantly, labeling what constitutes a stimulation artifact). Our post-hoc analysis also confirmed that the machine-learning refined labeling largely corresponded to the activation detection results based on the preset biological criteria, among which, the criterion based on significant ITPC (i.e., consistency in the evoked responses over trials) had the highest correspondence (0.80). In contrast, the traditional activation detection criterion of large peak detection (> 5 z-score) had a high precision (0.99) but low sensitivity (0.53), which is in line with our observation: many thalamic evoked responses (especially the contralaterally recorded ones) would be missed if based solely on the amplitude criterion as those responses usually have lower amplitudes on the voltage traces, though not necessarily possessing less robust evoked features. Our new approach yielded important results suggesting that the complex electrophysiological signals evoked by repeated single pulse electrical stimulations have three distinct features (i.e., F1, F2, and F3). In the context of past literature, unitary F1 and F2 features may correspond to clear N1 and N2 peaks, but we emphasize that changes in the neural activity of a given brain site evoked by the stimulation of another site cannot always be depicted as a uniform waveform consisting of N1 and N2, but rather, a complex signal with different components of the evoked waveform which may in turn reflect distinct physiological substrates with distinct sources. While we postulate that F1 and F2 may be related to the well-known N1 and N2 components of the EP signal described in the literature 12,14,15 , our data suggests that F3 is a unique feature of its own and not a prolonged F2 – even though F2 and F3 features have seemingly similar power spectrogram peaking at a low frequency range. Our rationale is as follows: (i) F3 has a peak oscillation at about 5 Hz, but F2 has wider frequency range with jittered peaks across trials. (ii) F3’s ITPC spectrogram differs significantly from the F2’s, as F3 has a phase-locking feature at the theta frequency but F2 does not (Fig. 4 a). The oscillatory feature of F3 was also observed at the same time-frequency cluster as long-lasting theta oscillations ( > = 2 cycles). The oscillatory feature of F3, unlike F2, was observed in the time domain: it was strongly time-locked and exhibited at least two cycles of oscillation - even when averaging the EPs recorded across all recording areas and across all thalamic stimulation sites. (iii) The temporal window of significant time-frequency clusters in F2 and F3 do not overlap (i.e., F3 always happens after F2, peaking > 200 ms post-stimulation). (iv) Lastly, F3 feature is present specifically when the thalamus is stimulated (including both THAL-ipsi and THAL-contr clusters), while F2 does not have this anatomical specificity. We documented that the F3 feature was detected with the stimulation of all three thalamic sites . This suggests that the F3 signal is not produced by the stimulation of only one specific thalamic nucleus. However, it's important to acknowledge a limitation of our study, which stemmed from the absence of data from a number of other thalamic nuclei. For instance, we have no information about the profile of electrophysiological responses evoked by the stimulation of thalamic nuclei such as the geniculate bodies. This limitation was primarily due to the clinical constraints and ethical considerations inherent in studies involving human subjects. We also note that the F3 oscillatory evoked responses were also seen with the stimulation of a few non-thalamic brain regions - especially the hippocampus and amygdala – both of which are non-cortical structures [Figure 6 b]). interestingly, the latency of F3 signal evoked by the stimulation of these non-thalamic sites was significantly delayed compared to the latency of F3 signal evoked by the stimulation of the thalamic sites, and they were mostly present in ipsilateral recording sites while the thalamic stimulations evoked F3 in far more widespread areas involving bilateral hemispheres. As the thalamus is considered to play a key role in cortical rhythms 16 , and as the hippocampus is known for its role in generating theta rhythms 17 , we hope our findings will motivate future systematic studies to determine how theta oscillations observed throughout a broad surface of higher association areas in the cerebral cortex during cognitive activities tasks or rest 18,19 are related to the delayed modulatory effects exerted by the thalamic (or medial temporal lobe) structures. It also remains to be determined if the bilateral presence of the F3 signal evoked by thalamic stimulations is functionally related to, or is crucial for, the integration of information across the two hemispheres. Our approach with signal decoding of the stimulation and recording data revealed that stimulation of a given cortical region evoked earlier and stronger responses in the thalamus compared to the cortical targets. A hypothetical interpretation of this finding is that the thalamus acts as a universal receiver of information being exchanged across cortical regions (i.e., a “listener of corticocortical dialogue”) – which is a necessary for it to function as a key structure for information integration as detailed recently 2 . Additionally, the multisite thalamic stimulation and recording within the same individual offered us an unprecedent opportunity to examine the intrathalamic and interthalamic causal connections. We demonstrate that the stimulation of a given thalamic site led to fast high frequency modulation of activity in other thalamic sites ipsilaterally as well as delayed low frequency modulation of thalamic sites contralaterally. These findings provide a plausible electrophysiological mechanism for intra-thalamic and inter-thalamic connectivity, which may have important implications and clinical relevance: Stimulating a specific thalamic nucleus with DBS (e.g., for treatment of epilepsy 20 ) is likely to impact the activity of other thalamic nuclei (beyond the one directly targeted) and the cortical networks they are associated with. Thus, the therapeutic benefits of DBS might be related to a broader network of brain regions that are modulated rather than a focal network targeted. However, this does not imply that the effect of stimulation of each targeted thalamic site is diffuse and all-encompassing, as shown in Fig. 5 , stimulations in three subregions of the thalamus do not elicit anatomically identical responses, but instead show relative anatomical preference. We are mindful that our data was acquired from patients with a neurological disease (epilepsy) and one may ask to what extent they are generalizable to normal human brains. While we acknowledge that molecular and cellular changes accompany local circuit rewiring in patients with chronic epilepsy 21 , there is no firm evidence suggesting changes in the fundamental topology of connections at the system level in these patients 22 . In the extant literature, unless studying the circuit level and local changes in a region-of-interest type of analysis 23–25 or graph theoretic measures for deriving high-level estimation of network capacity 26,27 , the global network topology of the brain has been shown not to be significantly different between epileptic cohort and the normal control 28 . Several studies focusing on epileptic brain reorganization have shown that the same sets of canonical brain networks can be discovered by the data-driven methods, just as in normal controls 29–31 . We are also mindful that electrode coverage in our subjects was limited and sparse. This is unfortunately a limitation of the intracranial approach in human participants that we, as researchers, have little control over since the implantation of electrodes in all participants was clinically driven. We made the best effort to reduce the confounding effects by maximizing the number of subjects, and diversifying the cohort with different epilepsies so as to avoid a system bias of encountering pathological tissues concentrated in specific areas. Despite the limitations of our study, we hope that our findings can serve as the steppingstone towards a wider understanding of the architectural topology of interactions across the cerebral cortex and between those and subcortical structures such as the thalamus. In closing, we acknowledge that the causal connectivity maps presented in our current study provide only one aspect of the functional architecture of the brain, and as such, they serve as rough approximations for the dynamic causal interactions occurring in real-life scenarios and during cognitive processing and are complimentary to those connectivity measures generated by structural diffusion MRIs or resting state or task specific imaging studies 3,4,32–34 . To highlight the functional significance of RSEPS-based connectivity maps, we recently confirmed that the strength of RSEPS connectivity between sites correlates with the strength of their co-activations during a cognitive task, i.e., if a specific pairs of neuronal populations across the cortex and thalamus had high RSEP-based connectivity strength they had higher correlation of high-gamma activations during a memory task 11 . Additionally, it has been shown that seizures propagate from cortical seizure onset zones to thalamic sites 5 or other cortical areas 35 that show faster and stronger RSEPS connectivity measures with the seizure onset zone. Online Methods Resource Availability Materials availability Raw data, including electrophysiology collected in this study will be shared publicly after publication. Data and code availability: The electrophysiological data have been deposited at Mendeley and are publicly available as of the date of publication. The DOI will be listed in the key resources table. All original codes have been deposited at Zenodo and is publicly available as of the date of publication. DOIs will be listed in the key resources table. Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request. Participants In this study, twenty-seven participants (40.7% female; mean age ± SD: 34.9 ± 10.0 years) diagnosed with focal epilepsy were recruited (see Table 1 ). Each participant had 180 ± 46 (mean ± standard deviation) SEEG contacts implanted. The total number of electrode contacts being studied is 4864. All participants underwent invasive electrophysiological monitoring at our medical center as part of their treatment for refractory epilepsy. The placement of electrodes was determined exclusively based on clinical requirements. Prior to their involvement, all participants provided written informed consent, and the study protocol was approved by the Institutional Review Board for human experimentation. Table 1 Participant demographics Participants Age Gender Race Hemisphere coverage Number of electrodes Seizure onset zone Epilepsy (year) S01_166 38 F White B 176 right hippocampus 2 S02_167 23 F Mixed R 168 right posterior insula 4 S03_169 23 M White B 186 inferior occipital gyrus 19 S04_170 40 M Asian L 116 left posterior inferior temporal lobe 24 S05_171 52 M Asian B 252 left lateral temporal and right orbitofrontal region 35 S06_172 47 M Asian B 174 mesial temporal and insula 14 S07_176 35 F White L 156 left hippocampus 2 S08_177 19 M White B 168 hippocampus 4 S09_178 28 M White B 252 right anterior temporal lobe 9 S10_181 23 F White B 248 right occipital lobe 5 S11_182 37 F White B 260 hippocampus 27 S12_183 33 M White B 174 right anterior temporal lobe 5 S13_185 36 F White B 142 left mesial temporal lobe 26 S14_188 31 M White B 122 hippocampus 10 S15_189 27 M Hispanic B 184 left posterior insula and left hippocampus 19 S16_190 50 F White R 116 right primary somatosensory cortex 4 S17_192 20 M Hispanic L 88 left hippocampus 5 S18_193 40 M Hispanic B 182 right periventricular nodular heterotopia, right occipital region 30 S19_194 37 M Hispanic B 256 left medial hippocampus, right posterior insula and pulvinar 33 S20_195 47 F Black B 138 mesial temporal lobe and insula 17 S21_196 40 F White B 158 left orbital frontal and left mesial temporal lobe 20 S22_197 20 M Black B 214 right supramarginal area 15 S23_198 38 M Black B 204 mesial temporal lobe 4 S24_199 36 M White B 178 right superior temporal gyrus 20 S25_201 27 M Hispanic B 174 mesial temporal lobe 7 S26_202 48 F Asian B 194 left mesial occipital area 16 S27_205 48 F White B 184 right anterior hippocampus 8 M, male; F, female; R, right; L, left; B, both. Patient Safety & Ethics As noted, we did not implant extra electrodes for thalamic recordings. We only extended the electrodes that were clinically planned for cortical monitoring to reach the thalamus. Our procedure was based on patients’ informed consent and all our research related activities were IRB approved. To minimize injury, we used a reduced diameter obturating stylet and reduced diameter electrodes with 0.86mm diameter (AdTech, Inc). With this surgical approach, we have observed no complications, no thalamic hemorrhage or edema, and no neurological symptoms post operatively in any of the patients examined. We note that the thalamic recording in epilepsy patients has been widely practiced in some centers in Europe for two decades yielding clinically important information 36–40 . The method has recently gained ground in the U.S. 41–44 , and a recent neurosurgery editorial recommended and encouraged more centers in the U.S. to do the same 45 . Additionally, iEEG implantation strategy is decided by a group of clinicians who meet weekly to decide the clinical needs of the patient, and not by this research project. Also, as noted, all patients provide informed consent prior to being enrolled in any research studies. Electrode Implantation MRI sequences for imaging the thalamus and its subnuclei included BRAVO, fGATIR and MP2RAGE. We co-registered the post-implant CT with a pre-implant MRI. High resolution T1, fast gray matter acquisition T1 inversion recovery (FGATIR), and T1 post-contrast imaging were used for planning. Trajectories were planned to traverse in an orthogonal plane to capture the cortical frontal or temporal operculum, insula, and to be extended into specific sites of interest in the thalamus. The priority for all trajectories was safety and avoidance of blood vessels. To achieve implementation of our multisite sampling approach, the cortical electrode trajectories were extended to reach three different thalamic sites: anterior, mid, and posterior thalamus as detailed in our recent publication 5 . Approximate locations and number of electrodes along with their trajectories were all planned in a multidisciplinary surgical epilepsy conference with detailed review of presurgical data leading to a clinical hypothesis of the most likely seizure onset zones. Electrode placements were at the discretion of the clinical care team, who were not specifically involved with this research project. Co-localization of electrodes A thin cut head-CT, obtained after electrode implantation, was co-registered to pre-operative 3Tesla MRI data for verification of the trajectory. Precise electrode positioning in the FreeSurfer surface space, voxel space and MNI space were automatically extracted by the iElVis toolbox 46 . The T1-weighted MRI scan was used to generate 3D cortical volume and subcortical segmentation using recon-all command of Freesurfer v7.3.2 47 . The post-implant CT scan were co-registered to the pre-implant MRI using the flirt from the Oxford Centre for Functional MRI of the Brain Software Library 48,49 or using bbregister from Freesurfer 50 to get the best results. We manually labelled each electrode on the T1-registered CT image using BioImageSuite 51 . The electrode coordinates in the native anatomical space were inspected and manually labeled by the experienced neurologist J.P., based on the individual brain’s morphology and landmarks. Thalamic parcellation To check the concordance between manual and automatic labeling of thalamic electrode sites, individual contact centers of electrode masses were defined in native T1 space for each subject. The center of mass of each contact was then converted into an X,Y,Z coordinate in the MNI space. A recently developed and widely used MNI thalamic atlas 52 was utilized to parcellate the nuclear location of each contact within the thalamus. For this, a 1mm cubic voxel region was created around each contact center of mass (contact neighborhood). For each contact neighborhood, the fraction of voxels that overlapped with each thalamic nucleus in the THOMAS atlas was calculated. Contact neighborhoods that fell solely within a single THOMAS nuclear mask would have a value of 1 for that specific nucleus, while contact neighborhoods that had no overlapping voxels with a given nuclear mask would have a value assigned to 0 for that specific nucleus. Thus, for every contact neighborhood anatomically inside the thalamus, the fractional overlap with each THOMAS nuclear mask was calculable. For each site, this within-subject fractional overlap was summed across subjects and across all thalamic contacts to generate overall contact neighborhood fractional overlap values. Identifying Regions of Interest All electrodes contact locations were manually labelled with brain regions of interest (Fig. 1 a) by the experienced neurologist J.P. on the team. This process involves carefully inspecting the CT-MRI fused images in axial, sagittal, and coronal planes for the implanted electrodes, and their surrounding anatomical landmarks in the brain morphology (e.g., basal ganglia including putamen, and pallidum etc.). Contacts entirely or partially located in the white matter (e.g., internal capsule, or cingulum etc.) were labelled as white matter. We emphasize that these regions of interest were not identified in the standard atlas space, but in each participant’s own high-resolution T1 scan to ensure anatomical precision of electrode localization. Repeated single electrical pulse stimulation (RSEPS) Stimulations were performed in a bi-polar manner stimulating between a pair of electrodes. Each RSEPS trial entailed an instantaneous (pulse duration = 0.2 ms) 6 mA biphasic square-wave pulse. For electrode contacts near the seizure onset zones (identified by the clinical team), we stimulated at 4 mA. Repeated pulses (n = 42 ± 2 trials) were delivered in a row with an interval of 2 seconds. We excluded data from (i) bipolar channels with both contacts located in the white matter (see above for details of anatomical localization), (ii) bipolar contacts not within the same anatomical region (e.g., one contact in insula and the other in orbitofrontal cortex), (iii) recording bipolar sites too close to the stimulation site (Euclidean distance between their middle points = 7 SD of the baseline) only in the artifact window (< 10 ms), and (v) bad channels based on more specific criterion as detailed below. Intracranial EEG data acquisition and preprocessing The data was collected using the Nihon Kohden recording system with a sampling rate of 1000 Hz. We used an in-house programmed preprocessing pipeline to denoise the intracranial EEG data before any statistical testing. The pipeline included notch filtering at 60, 120 and 180 Hz, data exclusion, and re-referencing. Excluded data includes pathological, noisy channels, and bad trials, the identification criteria of which are described as follows. We used a recently developed algorithm to identify pathological channels by the presence of pathological high-frequency oscillations (HFOs) 53 . Noisy channels were identified by the presence of extremely high raw amplitude (> 5 standard deviations [SD] across all channels) and/or the prevalence of spikes (> 3 median of the distribution across all channels), i.e. jumps between consecutive data points larger than 80 µV. Trials with extreme variations in the signal, i.e. having extreme mean (> 4 SD across all trials) of the absolute magnitude of the signal values, and extreme standard deviation (> 3.5 SD) of the variance of the signal values along the sampled timepoints, were excluded for being trial-averaged to create the time-domain evoked potential. Following data screening, we applied bi-polar re-referencing instead of common averaging, given that the stimulations were performed in a bi-polar manner, i.e., injected between a pair of electrodes. Common average reference of all non-noisy channels for the intracranial EEG recording analysis has also been experimented and compared with the bi-polar referencing scheme to make sure that the bi-polar subtraction did not remove signals. After manual inspection, we chose to use the bi-polar referencing for the further analyses. To extract temporal-frequency information, the continuous bipolar voltage traces were decomposed using Morlet wavelet filtering (at log-spaced frequencies between 1 and 256 Hz [59 total frequencies]; each wavelet having a width of 5 cycles). The output instantaneous power timeseries has a sampling rate of 200 at each frequency. The timeseries was then epoched, baseline corrected and averaged across good trials. Values in the power spectrogram were then log-transformed to approximate normal distribution and z-scored across time and frequency. To measure how EPs are consistent across trials, we calculated the inter-trial phase coherence (ITPC) by taking the complex phase component of the wavelet decomposition, averaging across trials, and then taking the absolute. Values in the ITPC spectrogram were then square-root transformed to approximate normal distribution and then z-scored across time and frequencies. The remaining text of the paper will refer to the evoked potentials (EP), the evoked power and ITPC spectrograms as their preprocessed version by default. Activation detection using semi-supervised Uniform Manifold Approximation and Projection (UMAP) Preliminary activation labelling The key to determining a casual connectivity lies on correctly detecting a true physiological activation among the volume conductions of an electrical stimulation. We employed a semi-supervised learning approach utilizing the Uniform Manifold Approximation and Projection (UMAP) algorithm with the python toolbox umap-learn ( https://umap-learn.readthedocs.io/en/latest/ ) 54 . Unlike traditional dimensionality reduction techniques, semi-supervised UMAP allows for integration of labeled and unlabeled data, using our partial knowledge to guide the manifold structure of a true activation. Before training, we produced a preliminary activation label based on the “consensus” of the pre-defined 3 criteria: (1) “Uniform-Cutoff”: a hard threshold of EP z-score > 7 at any time point in the [10, 600] ms time window. (2) “Channel-Adaptive”: stimulation-channel specific thresholds (SD > 2) for peak heights and peak dominances of the EPs generated from the same bipolar stimulation. EP peaks were detected using the MATLAB “findpeaks” function for all peaks (positive or negative), and the peak and prominence distributions were constructed from all the detected peaks, upon which we applied the thresholds. (3) “Phase-locking”: individual specific thresholds with a confidence level of 95% in the distribution of connectivity index scores across all stimulating-recording pairs of the same subject. The connectivity index is based on ITPC, recently created to measure causal connectivity as discussed in a recent study 11 . Apart from these criteria reported from the literature, we also imposed some criteria that should be met to be physiologically alike: the voltage trace of the EP should have smooth and fast rises in the early phase (10 ~ 60 ms), later peaks (after 250 ms) can have lower peak heights (z-score > 5), and the signal should approximate to baseline after 800 ms. When all the criteria agree, we labelled a stimulating-recording instance as either “1” (activated) or “0” (non-activated); otherwise (~ 1/3 of the cases), we labelled it as “-1” (unlabeled), which is left for the algorithm to apply learned rules and decide for us. Training data preparation We applied a semi-supervised UMAP with power and ITPC spectrograms as input data, as the spectrograms explicitly quantified the rich waveform information (components and consistency) of the EPs. We prepared the training data with the following steps. We down sampled and smoothed the spectrograms by averaging the nearest neighbors to reduce the data redundancy, save working memory for the intensive computation, and to denoise and highlight the features of interest. We used characteristic time-frequency kernels to smooth (averaging the nearest neighbors) and down-sampled the spectrograms, based on Leland McInnes k-means clustering from all stimulation-evoked potential peaks, and our prior knowledge that the first 10 ms upon a stimulation onset is mainly artifact ( Fig. S5 ). Therefore, we down sampled the temporal domain into 5, 25, 20, and 10 points to cover the windows of [-30, 10] ms (artifact window), [11, 117] ms, [118, 283] ms, and [284, 800] ms. The artifact window was kept in the training data for the algorithm to learn to recognize. We used more data points for representing the [11, 283] ms because most evoked potentials happen in this time window. On the frequency domain, they were down sampled into 4 datapoints to cover the [0.5, 5] Hz (delta), [5, 8] Hz (theta), [8, 15] Hz (alpha), [15, 30] Hz (beta), [30, 70] Hz (gamma) and > 70 Hz (high). The down sampled power and ITPC spectrograms were then vectorized and concatenated as a single vector of multiple features (number of features = timepoints*frequencies of power + timepoints*frequencies of ITPC) for each stimulating-recording instance. Each stimulating-recording vector was then concatenated to be the input training data. The input data was curated using the sklearn toolbox in python ( https://scikit-learn.org/stable/ ), including missing data imputation with mean statistics (negligible cases, about 0.089% of the datapoints), and quantile transformation to distribute the data normally. When both transformed and un-transformed data were experimented, the results did not change but showed a clearer clustering effect using the transformed data; hence, the results using transformed data are reported in this paper. Subject-level semi-supervised UMAP learning The number of latent dimensions for the UMAP mapping was set to 2, the size of the local neighborhood was set to 5, the minimum distance of the points in the latent space was set to 0.01. The rest of the parameters were kept at default. As a result, each participant’s EP data was classified into activation and non-activation clusters that were grouped at the far ends in latent space, while the few floating data points that were not clustered were considered outliers. Each person’s clustering was manually checked by (1) visualizing the spectrograms of the centroid and boundary points of each cluster, and (2) visualizing the inverse projection of the original space from the latent space. 26 out of 27 participants’ data were grouped into two meaningful clusters, while one participant’s data (S21_196) with extensive amount of noisy data failed to cluster and was excluded from the following analyses. Neural feature extraction for activated signals Group-level UMAP learning To gain an overview of the data structure of the activated spectrograms, we applied an unsupervised UMAP with input data of denoised power and ITPC spectrograms of the activated channels (number of latent dimensions = 4). To minimize the influence of electrical artifacts, we created a non-activated spectral information (power and ITPC) template for each stimulating channel by averaging all the non-activated pairs stimulating from this channel. This non-activated spectral template contains mostly the artifact of volume conduction ( Fig. S5 ), as the physiological events would be averaged out when not synchronized by an activation. We used the non-activated spectral template as a control and subtracted it from all the activated spectrograms from the same stimulation channel within the same participant. When the non-activated cases for a stimulation channel are too few (< 10) and thus unable to generate a stable mean as a contrast, this happened only once, we excluded the case from the subsequent analyses. The remaining traces of the activated spectral information were considered denoised and fed into the group-level analysis. We applied a supervised UMAP with input data of denoised power and ITPC spectrograms of the activated channels, to map to the anatomical labeling of interest. The provided anatomical labels were “THAL-ipsi”, “THAL-contr”, “COR-ipsi” and “COR-contr”, corresponding to stimulating from thalamus/cortex and recording from the ipsilateral/contralateral hemisphere. The same training data preparation as mentioned before was applied. The number of latent dimensions for the UMAP mapping was set to 2, the size of the local neighborhood was set to 15, the minimum distance of the points in the latent space was set to 0.1. These parameters varied from the within-subject semi-supervised UMAP learning because the data variability within-subject is assumed to be smaller than the group data, while we aim to find commonalities that override individual differences. Therefore, the local neighborhood and the latent space are relaxed to bigger numbers (As a general principle, we experimented with different sets of parameters [within a reasonable range] to achieve the best clustering results from visual inspection). The rest of the parameters were kept at default. We further employed a hierarchical-density clustering algorithm to define the boundaries of the embedding clusters against the background noise. Finally, we investigated the centroid and averaged the presentation of the data in each cluster to verify that the clustering preserved meaningful and distinct features within the data. To ensure the data-driven clusters are generalizable and meaningful, we fed the resulting embedding to a series of feature evaluation procedures, including using a support vector classifier to test their generalizability to the ¼ unseen testing data. We showed that the embedding for separating the cortical and thalamus evoked responses has 75.9% predictability to the unseen data, much higher than chance rate (25%). Especially, the clusters of thalamus and cortex stimulated data are mostly well preserved, suggesting these two types of connectivity show most distinctive features as encapsulated by their power and ITPC spectrograms. We further employed a hierarchical-density clustering algorithm 55 to define the boundaries of the embedding clusters against the background noise, using the HDBscan toolbox ( https://hdbscan.readthedocs.io/en/latest/ ). We then investigated the averaged spectrograms (i.e. the cluster centroid) for each cluster, for a sanity check of the clustering results. Cluster-based permutation test Category-specific spectral features were determined by cluster-based permutation significance testing: (1) one-sample T-test for determining the significant clusters for each category, and (2) paired T-test for determining the significant differences between categories (with a 2 [THAL/COR]-by-2 [ipsi/contr] design). Specifically, two-sample t-tests were conducted within subjects and the individual statistics (t-scores) were fed into a group-level one-sample T-test for final inference. The testing was performed using the toolbox MNE-python (v1.5.0) ( https://mne.tools/stable/index.html ). The number of permutations was set to 5000, the initial cluster forming threshold was bigger than 6 t-scores, and the cluster-wise permutation confidence interval was set to 99% (i.e., P cluster < 0.01). Neural feature decoding Based on the cluster-based permutation test, we identified three distinct neural features (i.e., significant clusters) that are specific to each anatomical category of interest (i.e., UMAP localizers). They are represented by the power and ITPC patterns in the window of ([10, 50] ms of the COR-ipsi category (averaged across all instances in the corresponding category), in the window of [70, 160] ms of the COR-contr category, and the [165, 280] in the THAL-contr category. The temporal-frequency relationships in these time ranges and anatomical categories were considered as neural features. They were decoded in each instance of stimulated/record pairs. To do this, we utilized the spectral information from the significant cluster in the group-level spectrograms as a template. We employed a sliding-window cross-correlation to match this template with the specific instances of the spectrograms, point by point in time. Different similarity (or distance) measures were experimented and turned out to produce similar results; therefore, we used the Pearson correlation r for its easy interpretability. The resulting timeseries of r has a samping rate of 200 per second, which is the same as the spectral dataset. A sliding-window cross-correlation with the 3 identified neural features generated three temporal similarity curves for each stimulating-recoding instance, each corresponding to the neural feature’s temporal evolvement in this specific pair of stimulated/recorded sites. With the temporal similarity curve, we detected the peak similarity and the time to the peak using MATLAB “findpeak” function. Instances within each anatomical area were averaged for brain regional inferences. General framework of statistical testing: The datasets in the present study generally have a nested structure such that each condition consists of multiple electrode contacts, each with a brain area and subject identity; therefore, the mixed linear model was used as a general testing framework for model fitting, using the brain area and subject identifiers as nested random factors. Therefore, the comparisons were always made within the same stimulation/recording channels, brain areas and subjects, whenever possible. Practical model setups can vary slightly according to data structures and tested hypotheses. For significance inference, the model comparison approach with a hypothesis testing framework was adopted. Namely, a null model is compared with an alternative model which has the hypothesized effect added on top of the null model. The model with higher evidence is selected, based on a likelihood-ratio test. We additionally provided the Akaike information criterion (AIC) and Bayesian information criterion (BIC) in our supplementary table of post-hoc tests where spurious significant results may arise given the large numbers of tests, even after the false-discovery-rate (FDR) correction. The model fittings and selections were conducted using the lme4 toolbox and the anova function in R. The model specifications were provided along with the significance tables in the supplementary materials. Declarations Acknowledgments: We are grateful to the patients who agreed to participate in our research, and the staff in the EEG lab including clinical fellows and attendings whose help enabled this research. We express our gratitude to our colleagues Riki Matsumoto (Kobe - Japan), Corey Keller (Stanford - USA) and Konrad Körding (Philadelphia – USA) for their valuable feedback on our findings; and to Leland McInnes (Ottawa – Canada) for running a sanity check on our UMAP analysis approach. Funding: This work was supported by research grants from the US National Institute of Neurological Disorders and Stroke (R01NS078396 and R21NS113024), US National Institute of Mental Health (1R01MH109954, and P50MH109429) and US National Science Foundation (BCS1358907) and US National Science Foundation (BCS1358907 and BCS1850938) to JP. Author contributions: DL contributed to the conception, acquisition, analysis, and interpretation of data and creation of new codes and measures used in the work and writing of the manuscript; JRS contributed to the analysis and interpretation of data; ZL contributed to the acquisition of data; JP contributed to conception and design of the work, acquisition and interpretation of data and writing of the manuscript. All authors provided significant feedback throughout the study and manuscript preparation and have approved the submitted version. Competing interests: The authors declare no competing interests. References Sporns, O., Chialvo, D. R., Kaiser, M. & Hilgetag, C. C. Organization, development and function of complex brain networks. Trends Cogn. Sci. 8 , 418-425 (2004). Shine, J. M., Lewis, L. D., Garrett, D. D. & Hwang, K. The impact of the human thalamus on brain-wide information processing. Nature reviews. Neuroscience (2023). https://doi.org:10.1038/s41583-023-00701-0 Matsumoto, R. et al. Functional connectivity in the human language system: a cortico-cortical evoked potential study. Brain 127 , 2316-2330 (2004). Matsumoto, R. et al. Functional connectivity in human cortical motor system: a cortico-cortical evoked potential study. Brain 130 , 181-197 (2007). Wu, T. Q. et al. Multisite thalamic recordings to characterize seizure propagation in the human brain. Brain 146 , 2792-2802 (2023). https://doi.org:10.1093/brain/awad121 Zumsteg, D., Lozano, A. M. & Wennberg, R. A. Depth electrode recorded cerebral responses with deep brain stimulation of the anterior thalamus for epilepsy. Clin Neurophysiol 117 , 1602-1609 (2006). https://doi.org:10.1016/j.clinph.2006.04.008 Ojeda Valencia, G. et al. Signatures of Electrical Stimulation Driven Network Interactions in the Human Limbic System. The Journal of Neuroscience: The Official Journal of the Society for Neuroscience 43 , 6697-6711 (2023). https://doi.org:10.1523/JNEUROSCI.2201-22.2023 Miller, K. J., Muller, K. R. & Hermes, D. Basis profile curve identification to understand electrical stimulation effects in human brain networks. PLoS Comput Biol 17 , e1008710 (2021). https://doi.org:10.1371/journal.pcbi.1008710 Jedynak, M. et al. Variability of Single Pulse Electrical Stimulation Responses Recorded with Intracranial Electroencephalography in Epileptic Patients. Brain Topogr 36 , 119-127 (2023). https://doi.org:10.1007/s10548-022-00928-7 Becht, E. et al. Dimensionality reduction for visualizing single-cell data using UMAP. Nat Biotechnol (2018). https://doi.org:10.1038/nbt.4314 Stieger, J. R. et al. Cross regional coordination of neural activity in the human brain during autobiographical self-referential processing. PNAS In Press (2023). Kunieda, T., Yamao, Y., Kikuchi, T. & Matsumoto, R. New Approach for Exploring Cerebral Functional Connectivity: Review of Cortico-cortical Evoked Potential. Neurologia medico-chirurgica 55 , 374-382 (2015). https://doi.org:10.2176/nmc.ra.2014-0388 Keller, C. J. et al. Intrinsic functional architecture predicts electrically evoked responses in the human brain. Proceedings of the National Academy of Sciences 108 , 10308-10313 (2011). https://doi.org:10.1073/pnas.1019750108 Keller, C. J. et al. Mapping human brain networks with cortico-cortical evoked potentials. Philos Trans R Soc Lond B Biol Sci 369 , 20130528 (2014). https://doi.org:10.1098/rstb.2013.0528 Togo, M. et al. Distinct connectivity patterns in human medial parietal cortices: Evidence from standardized connectivity map using cortico-cortical evoked potential. NeuroImage 263 , 119639 (2022). https://doi.org:10.1016/j.neuroimage.2022.119639 Fogerson, P. M. & Huguenard, J. R. Tapping the Brakes: Cellular and Synaptic Mechanisms that Regulate Thalamic Oscillations. Neuron 92 , 687-704 (2016). https://doi.org:10.1016/j.neuron.2016.10.024 Buzsaki, G. Rhythms of the Brain . (Oxford University Press, 2006). Groppe, D. M. et al. Dominant frequencies of resting human brain activity as measured by the electrocorticogram. Neuroimage 79 , 223-233 (2013). https://doi.org:10.1016/j.neuroimage.2013.04.044 Hacker, C. D., Snyder, A. Z., Pahwa, M., Corbetta, M. & Leuthardt, E. C. Frequency-specific electrophysiologic correlates of resting state fMRI networks. Neuroimage 149 , 446-457 (2017). https://doi.org:10.1016/j.neuroimage.2017.01.054 Fisher, R. et al. Electrical stimulation of the anterior nucleus of thalamus for treatment of refractory epilepsy. Epilepsia 51 , 899-908 (2010). https://doi.org:10.1111/j.1528-1167.2010.02536.x Engel, J., Jr. & Pitkanen, A. Biomarkers for epileptogenesis and its treatment. Neuropharmacology 167 , 107735 (2020). https://doi.org:10.1016/j.neuropharm.2019.107735 Caciagli, L., Bernhardt, B. C., Hong, S.-J., Bernasconi, A. & Bernasconi, N. Functional network alterations and their structural substrate in drug-resistant epilepsy. Front Neurosci 8 , 411 (2014). https://doi.org:10.3389/fnins.2014.00411 Fleury, M. et al. Episodic memory network connectivity in temporal lobe epilepsy. Epilepsia 63 , 2597-2622 (2022). https://doi.org:10.1111/epi.17370 Pittau, F., Grova, C., Moeller, F., Dubeau, F. & Gotman, J. Patterns of altered functional connectivity in mesial temporal lobe epilepsy. Epilepsia 53 , 1013-1023 (2012). https://doi.org:10.1111/j.1528-1167.2012.03464.x Roger, E. et al. Hubs disruption in mesial temporal lobe epilepsy. A resting-state fMRI study on a language-and-memory network. Human Brain Mapping 41 , 779-796 (2020). https://doi.org:10.1002/hbm.24839 Li, L. et al. Topographical reorganization of brain functional connectivity during an early period of epileptogenesis. Epilepsia 62 , 1231-1243 (2021). https://doi.org:10.1111/epi.16863 Mazrooyisebdani, M. et al. Graph Theory Analysis of Functional Connectivity Combined with Machine Learning Approaches Demonstrates Widespread Network Differences and Predicts Clinical Variables in Temporal Lobe Epilepsy. Brain Connect 10 , 39-50 (2020). https://doi.org:10.1089/brain.2019.0702 Ofer, I. et al. Association between seizure freedom and default mode network reorganization in patients with unilateral temporal lobe epilepsy. Epilepsy Behav 90 , 238-246 (2019). https://doi.org:10.1016/j.yebeh.2018.10.025 Liao, W. et al. Default mode network abnormalities in mesial temporal lobe epilepsy: A study combining fMRI and DTI. Human Brain Mapping 32 , 883-895 (2011). https://doi.org:10.1002/hbm.21076 Widjaja, E., Zamyadi, M., Raybaud, C., Snead, O. C. & Smith, M. L. Abnormal Functional Network Connectivity among Resting-State Networks in Children with Frontal Lobe Epilepsy. AJNR Am J Neuroradiol 34 , 2386-2392 (2013). https://doi.org:10.3174/ajnr.A3608 Zhang, Z. et al. Impaired attention network in temporal lobe epilepsy: A resting FMRI study. Neuroscience Letters 458 , 97-101 (2009). https://doi.org:10.1016/j.neulet.2009.04.040 Keller, C. J. et al. Mapping human brain networks with cortico-cortical evoked potentials. Philos Trans R Soc Lond B Biol Sci 369 (2014). https://doi.org:10.1098/rstb.2013.0528 Veit, M. J. et al. Temporal order of signal propagation within and across intrinsic brain networks. Proc Natl Acad Sci U S A 118 (2021). https://doi.org:10.1073/pnas.2105031118 Seguin, C. et al. Communication dynamics in the human connectome shape the cortex-wide propagation of direct electrical stimulation. Neuron 0 (2023). https://doi.org:10.1016/j.neuron.2023.01.027 Guo, Z. H. et al. Epileptogenic network of focal epilepsies mapped with cortico-cortical evoked potentials. Clin Neurophysiol 131 , 2657-2666 (2020). https://doi.org:10.1016/j.clinph.2020.08.012 Arthuis, M. et al. Impaired consciousness during temporal lobe seizures is related to increased long-distance cortical-subcortical synchronization. Brain 132 , 2091-2101 (2009). https://doi.org:10.1093/brain/awp086 Guye, M. et al. The role of corticothalamic coupling in human temporal lobe epilepsy. Brain 129 , 1917-1928 (2006). https://doi.org:10.1093/brain/awl151 Pizzo, F. et al. The Ictal Signature of Thalamus and Basal Ganglia in Focal Epilepsy: A SEEG Study. Neurology 96 , e280-e293 (2021). https://doi.org:10.1212/WNL.0000000000011003 Filipescu, C. et al. The effect of medial pulvinar stimulation on temporal lobe seizures. Epilepsia 60 , e25-e30 (2019). https://doi.org:10.1111/epi.14677 Evangelista, E. et al. Does the Thalamo-Cortical Synchrony Play a Role in Seizure Termination? Front Neurol 6 , 192 (2015). https://doi.org:10.3389/fneur.2015.00192 Gadot, R., Korst, G., Shofty, B., Gavvala, J. R. & Sheth, S. A. Thalamic stereoelectroencephalography in epilepsy surgery: a scoping literature review. J Neurosurg , 1-16 (2022). https://doi.org:10.3171/2022.1.JNS212613 Ilyas, A., Tandon, N. & Lhatoo, S. D. Thalamic neuromodulation for epilepsy: A clinical perspective. Epilepsy Res 183 , 106942 (2022). https://doi.org:10.1016/j.eplepsyres.2022.106942 Chaitanya, G. et al. Robot-assisted stereoelectroencephalography exploration of the limbic thalamus in human focal epilepsy: implantation technique and complications in the first 24 patients. Neurosurg Focus 48 , E2 (2020). https://doi.org:10.3171/2020.1.FOCUS19887 Romeo, A. et al. Early ictal recruitment of midline thalamus in mesial temporal lobe epilepsy. Ann Clin Transl Neurol 6 , 1552-1558 (2019). https://doi.org:10.1002/acn3.50835 McKhann, G. M. Editorial. Dulling the double-edged sword of human SEEG research. Neurosurg Focus 48 , E3 (2020). https://doi.org:10.3171/2020.1.FOCUS2069 Groppe, D. M. et al. iELVis: An open source MATLAB toolbox for localizing and visualizing human intracranial electrode data. J Neurosci Methods 281 , 40-48 (2017). https://doi.org:10.1016/j.jneumeth.2017.01.022 Fischl, B. FreeSurfer. NeuroImage 62 , 774-781 (2012). Jenkinson, M., Beckmann, C. F., Behrens, T. E., Woolrich, M. W. & Smith, S. M. Fsl. Neuroimage 62 , 782-790 (2012). https://doi.org:10.1016/j.neuroimage.2011.09.015 Jenkinson, M. & Smith, S. A global optimisation method for robust affine registration of brain images. Medical image analysis 5 , 143-156 (2001). Greve, D. N. & Fischl, B. Accurate and robust brain image alignment using boundary-based registration. Neuroimage 48 , 63-72 (2009). https://doi.org:10.1016/j.neuroimage.2009.06.060 Papademetris, X. et al. BioImage Suite: An integrated medical image analysis suite: An update. The insight journal 2006 , 209 (2006). Su, J. H. et al. Thalamus Optimized Multi Atlas Segmentation (THOMAS): fast, fully automated segmentation of thalamic nuclei from structural MRI. Neuroimage 194 , 272-282 (2019). https://doi.org:10.1016/j.neuroimage.2019.03.021 Liu, S. & Parvizi, J. Cognitive refractory state caused by spontaneous epileptic high-frequency oscillations in the human brain. Sci Transl Med 11 (2019). https://doi.org:10.1126/scitranslmed.aax7830 McInnes, L., Healy, J. & Melville, J. (arXiv, 2020). Campello, R. J. G. B., Moulavi, D. & Sander, J. (eds Jian Pei et al. ) 160-172 (Springer). Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software 67 , 1-48 (2015). https://doi.org:10.18637/jss.v067.i01 Additional Declarations There is NO Competing Interest. Supplementary Files SupplementaryMaterials.docx Cite Share Download PDF Status: Published Journal Publication published 15 Jul, 2025 Read the published version in Nature Neuroscience → Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-4366486","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":306310867,"identity":"ff9724a6-7867-42ef-be36-90eb19d4f37a","order_by":0,"name":"Josef Parvizi","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA60lEQVRIiWNgGAWjYBACxgYQacDAww/m2gDxAQYGCaK0SDYwAxlpRGiBA4MDxGphnt387DNPQZ2M8fnzBz8XJNjk8x1gPnibB5/D5hwzns1jwMZjdiOZWXpGQprlzANsydZ4tcxIMGacYcAD1MLMIM3747CBwQEeM2n8WtI/A7VI8Bj3H2b+zZMA0sL/jYCWHGOGDwYGPAYMyWzSEC08bPi1zDlTDNSSwCNxI9nMmichzUDyMJux5Rw8Wgxnt29mSPhTZ8/ff/DxbZ4EGwO+480Pb7zBp2UGhhAzHuUgIE9cRI+CUTAKRsGIBgCDKkREyuLACAAAAABJRU5ErkJggg==","orcid":"https://orcid.org/0000-0001-8520-7948","institution":"Stanford University","correspondingAuthor":true,"submittingAuthor":false,"prefix":"","firstName":"Josef","middleName":"","lastName":"Parvizi","suffix":""},{"id":306310868,"identity":"7997dfbc-c84d-4186-9cac-7be31ce118c6","order_by":1,"name":"Dian Lyu","email":"","orcid":"","institution":"Stanford University","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Dian","middleName":"","lastName":"Lyu","suffix":""},{"id":306310869,"identity":"680ee684-dcdf-40ba-a160-0278b2b580f2","order_by":2,"name":"James Stieger","email":"","orcid":"","institution":"Stanford University","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"James","middleName":"","lastName":"Stieger","suffix":""},{"id":306310870,"identity":"263e2568-c3c2-4128-8186-d7f7375a103f","order_by":3,"name":"Zoe Lusk","email":"","orcid":"https://orcid.org/0000-0002-4728-2688","institution":"Stanford University","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Zoe","middleName":"","lastName":"Lusk","suffix":""},{"id":306310871,"identity":"34cd6ae4-1751-4adc-bdcc-5a4d23e6884b","order_by":4,"name":"Vivek Buch","email":"","orcid":"https://orcid.org/0000-0003-1151-6128","institution":"Stanford University","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Vivek","middleName":"","lastName":"Buch","suffix":""}],"badges":[],"createdAt":"2024-05-04 03:25:37","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-4366486/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-4366486/v1","draftVersion":[],"editorialEvents":[{"content":"https://doi.org/10.1038/s41593-025-02009-x","type":"published","date":"2025-07-15T04:00:00+00:00"}],"editorialNote":"","failedWorkflow":false,"files":[{"id":57379907,"identity":"5c69b0a8-458e-42e5-b9e0-ff39f7707d11","added_by":"auto","created_at":"2024-05-30 01:12:32","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":196616,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eElectrophysiological responses (local field potential) evoked by the stimulation of a given brain area\u003c/strong\u003e\u003cem\u003e\u003cstrong\u003e \u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003eare\u003c/strong\u003e\u003cem\u003e\u003cstrong\u003e complex \u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003eand variable\u003c/strong\u003e\u003cem\u003e\u003cstrong\u003e. \u003c/strong\u003e\u003c/em\u003e\u003cem\u003eHere we show randomly selected examples of the complex waveforms generated by electrical stimulation. A large portion of evoked responses either does not reach an arbitrary threshold, e.g., z-score = 5 that is usually used in the literature (first column, with five separate examples) or does not conform to a typical shape of evoked responses with an N1 and N2 components (second column) or neither have a large amplitude nor a canonical N1 and N2 shape (third column). By using univariate measures, one will only map a fraction of true connections in the brain. Plotted signals are the evoked responses stimulated and recorded from a pair of bipolar sites (i.e., from one stimulation bipolar site to one recording bipolar site), of which the anatomical information is provided in the upper corner. The central line is the trials-averaged signal, baseline-corrected and z-scored, while the shaded area depicts the standard error over the 45 trials.\u003c/em\u003e\u003c/p\u003e","description":"","filename":"1.png","url":"https://assets-eu.researchsquare.com/files/rs-4366486/v1/6f72e69bd093f0ded9095505.png"},{"id":57379908,"identity":"47bdf54c-5731-4ee6-ae6f-57982ff29dc6","added_by":"auto","created_at":"2024-05-30 01:12:32","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":529211,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cem\u003e\u003cstrong\u003eIllustration of data processing pipeline. \u003c/strong\u003e\u003c/em\u003e\u003cem\u003e(a) Electrode coverage at the group level. The electrode localization was manually labelled by the experienced neurologist on the team, based on each subject’s brain morphology in their own high-resolution T1 scan (L: left, R: right, ant: anterior, pst: posterior, PMC: posteromedial cortex, SM: sensorimotor, SPL: superior parietal lobule, ACC: anterior cingulate cortex, CLT: claustrum, TH: thalamus, IPL: inferior parietal lobule, MFC: medial frontal cortex, OFC: orbital frontal cortex, STG: superior temporal gyrus, LFC: lateral frontal cortex, INS: insula, FG: fusiform gyrus, HPC: hippocampus, MTG: middle temporal gyrus, PHG: parahippocampal gyrus, AMY: amygdala, TP: temporal pole, MCC: midcingulate cortex, ITG: inferior temporal gyrus). For visualization purposes, each subject’s brain images were normalized to a standard brain surface space (FS_LR: https://osf.io/k89fh/wiki/Surface/) with FreeSurfer (https://surfer.nmr.mgh.harvard.edu/) and Connectome-Workbench (\u003c/em\u003e\u003ca href=\"https://www.humanconnectome.org/software/connectome-workbench\"\u003e\u003cem\u003ehttps://www.humanconnectome.org/software/connectome-workbench\u003c/em\u003e\u003c/a\u003e\u003cem\u003e), and their electrodes were projected to the standard surface (sites out of grey matter are excluded). \u003c/em\u003e\u003cem\u003e\u003cstrong\u003e(b)\u003c/strong\u003e\u003c/em\u003e\u003cem\u003e The localization of electrodes in the thalamus, divided into three thalamic divisions: anterior (purple), mid (blue), and posterior (green) subregions of the thalamus. \u003c/em\u003e\u003cem\u003e\u003cstrong\u003e(c)\u003c/strong\u003e\u003c/em\u003e\u003cem\u003e Roadmap of data analysis. \u003c/em\u003e\u003cem\u003e\u003cstrong\u003e(d) \u003c/strong\u003e\u003c/em\u003e\u003cem\u003eDimensionality reduction of the evoked spectral information of all pair-wise evoked potentials in the activated cluster (i.e., pairs with underlying connections). The spectral information is the concatenated power and inter-trial phase coherence (ITPC) spectrograms of the evoked potentials, which are down-sampled to balance the number of datapoints for earlier and later neural responses (see Methods). While this figure is not intended to show distinct clusters, we have colored the dots (a posteriori) to illustrate that data from ipsilateral thalamic stimulations (red dots) are already distinguishable upon a visual inspection. Notably, the color patterns follow different stimulation sites but for recording sites (also see SI \u003c/em\u003e\u003cem\u003e\u003cstrong\u003eFig. S3 a,c,d\u003c/strong\u003e\u003c/em\u003e\u003cem\u003e). They are not biased by specific subjects (SI \u003c/em\u003e\u003cem\u003e\u003cstrong\u003eFig. S4\u003c/strong\u003e\u003c/em\u003e\u003cem\u003e). \u003c/em\u003e\u003cem\u003e\u003cstrong\u003e(e)\u003c/strong\u003e\u003c/em\u003e\u003cem\u003eNeural feature encoding involved two main steps: UMAP supervised learning and cluster-based permutation testing. The input data was the evoked power and ITPC spectrograms of the group-level whole-brain evoked potentials, while the dependent variable (i.e., the labels) were the ipsilateral cortical (COR-ipsi), ipsilateral thalamic (THAL-ipsi), contralateral cortical (COR-contr), and contralateral thalamic (THAL-contr) pairs. The algorithm successfully characterized these evoked spectrograms into the four categories. HDBscan was used to formally define the clusters in the embedding space, dissociating them from noisy channels. Original spectrograms of the pairs clustered in the four categories were then used to perform cluster-based permutation testing, by which we identified three time-frequency clusters within the spectrograms that were specific to each category (i.e., Neural Feature). \u003c/em\u003eBefore statistical testing, the evoked power spectrograms of the four stimulation sites already show distinct features that can be distinguishable by visual inspection alone (marked by white circles and numbered).\u003c/p\u003e","description":"","filename":"2.png","url":"https://assets-eu.researchsquare.com/files/rs-4366486/v1/7f78a01700626fa456dff3c6.png"},{"id":57380369,"identity":"68d892de-aebc-46d6-8b09-17781150d835","added_by":"auto","created_at":"2024-05-30 01:20:32","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":249408,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cem\u003e\u003cstrong\u003eSpectrograms of stimulation-evoked power within and between anatomical categories\u003c/strong\u003e\u003c/em\u003e\u003cem\u003e (COR: cortex, THAL: thalamus, ipsi: ipsilateral, contr: contralateral, à: causal influence direction). Highlighted contours on the spectrograms indicate significant clusters (n-permutation = 5000, initial cluster forming threshold = 6 t-scores, P\u003c/em\u003e\u003csub\u003e\u003cem\u003ecluster\u003c/em\u003e\u003c/sub\u003e\u003cem\u003e \u0026lt; 0.01). Dashed lines on the spectrograms denotes the segments of conventional frequency bands: [0.5, 5] Hz (delta), [5, 8] Hz (theta), [8, 15] Hz (alpha), [15, 30] Hz (beta), [30, 70] Hz (gamma) and \u0026gt; 70 Hz (high). Same tests were conducted on ITPC spectrograms, which generated similar results (\u003c/em\u003e\u003cem\u003e\u003cstrong\u003eFig. S6\u003c/strong\u003e\u003c/em\u003e\u003cem\u003e). The frequency axis is in a logarithmic (log) scale. The time axis is unevenly sampled to balance the varied lengths for different clusters (see Method “training data preparation”). (a) Thalamic and cortical evoked spectrograms and significant clusters in either ipsilateral or contralateral hemispheres. Significance testing conducted using within-subject one-sample t-test (in a mixed-model design): individual-level t-statistics input into a group-level significance testing. The color bars show group-level t-scores. The color bars show group-level t-scores; significant clusters of the one-sided test is marked by white contours. These spectrograms, without significant markers, have also been used in \u003c/em\u003e\u003cem\u003e\u003cstrong\u003eFigure 2e\u003c/strong\u003e\u003c/em\u003e\u003cem\u003e for illustration. (b) Thalamic vs. cortical evoked spectrogram comparison using within-subject two-sample t-test (two-tailed) on the power spectrograms, significance inference performed with cluster-based permutation testing. Yellow and blue contours respectively highlight significant clusters of the contrast indicated on the subtitle and its reversed contrast. \u0026nbsp;Acronym example “THAL-ipsi” on the subtitle means ipsilateral pairs stimulating from thalamus. (c) Intra-thalamic, thalamocortical, corticothalamic and cortico-cortical evoked spectrogram comparisons. Acronym example “THAL-ipsi” means ipsilateral pairs stimulating from the thalamus. Significant clusters of the one-sided tests are marked by white contours. For THAL-\u0026gt;THAL (contra) connections, no significant cluster was generated here as there were not enough subjects (n\u0026lt;10) to provide sufficient statistical power for the mixed-model testing. However, a significance test can be done at the level of electrode contacts (i.e., without considering the grouping factor of “subject”). This alternative analysis showed a consistent pattern with the subject-averaged spectrogram, and resulted in two significant clusters, one in the gamma band before 60ms, another in the alpha/theta band around 100-275 ms (\u003c/em\u003e\u003cem\u003e\u003cstrong\u003eFig. S8\u003c/strong\u003e\u003c/em\u003e\u003cem\u003e). \u0026nbsp;Acronym example “THAL-THAL (ipsi)” on the subtitle means ipsilateral pairs stimulating from the thalamus and recording from the thalamus.\u003c/em\u003e\u003c/p\u003e","description":"","filename":"3.png","url":"https://assets-eu.researchsquare.com/files/rs-4366486/v1/f06acbddded1ded455dd6942.png"},{"id":57379913,"identity":"d9cc2655-b12e-4dda-8f2c-606d1c487e8f","added_by":"auto","created_at":"2024-05-30 01:12:33","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":518518,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cem\u003e\u003cstrong\u003eElectrophysiological neural features indicating different types of connectivity and illustration of the decoding process.\u003c/strong\u003e\u003c/em\u003e\u003cem\u003e (a) Neural features specified by the spectra-temporal information. The time windows of the neural features are circumscribed by the significant clusters distinct among conditions from the previous cluster-based permutation testing. The time-frequency relationships in the three features are represented by the group-averaged spectrograms of the COR-ipsi, COR-contr and THAL-contr stimulations, as these three categories show clearest and non-overlapped significant clusters at the group level. The frequency axis is in log-scale. The time axis is in a natural scale with even samples, different from the time axis in \u003c/em\u003e\u003cem\u003e\u003cstrong\u003eFigure 3\u003c/strong\u003e\u003c/em\u003e\u003cem\u003e. The values of power and ITPC spectrograms are transformed (logarithmized and square-rooted, respectively) to approximate Gaussian distributions, and then zscored (in both time and frequency directions) to be comparable among connections. Line plots to the right of the spectrograms show the (normalized and log-scaled) spectral density of the power and ITPC during the significant time window, with colored lines for all the time points, and black for the time average. Since the evoked spectrograms were baseline-corrected, depicted curves can be seen as the “bumps” on top of the 1/f background noise. Group-averaged evoked potentials at the temporal domain corresponding to each category of the spectrograms are shown below. (b) Cross-correlation with a sliding window approach. To reveal the presence of each feature in individual connections, the feature information (i.e. the short-lasting spectral information characterized by the power and ITPC in the time-frequency window) is correlated to the spectrograms with Person-correlation r. To examine this correspondence at every time point (sampling rate is 200 per second), an over-lapping sliding window approach is used, whereby the “transient” time-frequency information at each time point is examined while the examining window is sliding over the spectrograms. This generates the dynamic appearance of the feature presence over time. The line plot in the middle is an example case randomly selected form the THAL-contr instances. Every point on the curve, with a specific time in x-axis and similarity measure in y-axis, indicates how similar the current spectral information (power and ITPC sustaining in certain amount of time) matches the neural feature in concern. The maximum r (r*) and the time to r* (latency) was taken as indices of feature representation and analyzed in the subsequent analyses. The heatmap to the right shows the correlation between all the THAL-contr instances and the Feature 3 (F3). An overview of feature presence in all categories of all features is presented in \u003c/em\u003e\u003cem\u003e\u003cstrong\u003eFig. S9\u003c/strong\u003e\u003c/em\u003e\u003cem\u003e.\u003c/em\u003e\u003c/p\u003e","description":"","filename":"4.png","url":"https://assets-eu.researchsquare.com/files/rs-4366486/v1/9c105d6c739722462559cc93.png"},{"id":57379912,"identity":"50ee2576-25fa-4a14-95bc-743d7e5c3ef5","added_by":"auto","created_at":"2024-05-30 01:12:33","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":215641,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cem\u003e\u003cstrong\u003eThe timing of feature presence (i.e., latency of maximum representation of the feature). \u003c/strong\u003e\u003c/em\u003e\u003cem\u003e(a) The top panel shows the distribution of the latencies of F1 presence in all types of connections. F1 is identified in individual instances by finding the peaks of the feature curve. Included in the distribution are connections with strong F1 representation (r* \u0026gt; 0.4) which falls in a sensible latency range of (10, 70) ms, based on our results. Density distributions for thalamocortical connections are filled red, corticothalamic connections are filled blue, and thalamo-thalamic/cortico-cortical connections are filled grey. Colors for contralateral connections are lighter than those for ipsilateral connections. The middle panel is a cartoon showing the median latency of the feature presence among the aforementioned connections. The lower panel depicts within-subject post-hoc comparisons (adjusted for individual and regional differences) between the categories that showed evident differences in the latency distributions. Each colored dot marks the mean latency of all the connections from one subject. The black bar through the dot marks the range of ± one standard error of the within-subject data. Grey lines across groups show the within-subject comparisons. All tests were corrected for multiple comparisons. The full list of post-hoc comparisons is presented in SI Table S9. (b) and (c) shows F2 and F3 results with same visualization schemes. The criteria for data going into the F2 distributions are r* \u0026gt; 0.4 and latency in (70, 165) ms. The criteria for the F3 distributions are r* \u0026gt; 0.5 and latency in (200, 400) ms. As there were much fewer contralateral thalamocortical connections with F1 and F2, the data for these two features showed large variabilities.\u003c/em\u003e\u003c/p\u003e","description":"","filename":"5.png","url":"https://assets-eu.researchsquare.com/files/rs-4366486/v1/e57290617d5923794980dc61.png"},{"id":57379914,"identity":"8542c66f-8e32-4872-a891-6275fc60890a","added_by":"auto","created_at":"2024-05-30 01:12:33","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":299383,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eThe anatomical landscapes of feature presence.\u003c/strong\u003e (a) Early F1 connectivity matrices for bilateral cortical and thalamic subregions (l: left, r: right, ant: anterior, pst: posterior, PMC: posteromedial cortex, SM: sensorimotor, SPL: superior parietal lobule, ACC: anterior cingulate cortex, CLT: claustrum, TH: thalamus, IPL: inferior parietal lobule, MFC: medial frontal cortex, OFC: orbital frontal cortex, STG: superior temporal gyrus, LFC: lateral frontal cortex, INS: insula, FG: fusiform gyrus, HPC: hippocampus, MTG: middle temporal gyrus, PHG: parahippocampal gyrus, AMY: amygdala, TP: temporal pole, MCC: midcingulate cortex, ITG: inferior temporal gyrus); The anatomical localization of all electrode contacts was visually inspected and labeled by an experienced neuroanatomist using the position of the electrodes in the individual subject’s native brain space. Each entry of the matrix has a row and column identity corresponding to site of stimulation and recording, respectively, and its value indicates the feature representation (R) for the corresponding feature. Specifically, R = r\u003csup\u003e-\u003c/sup\u003e, was averaged across connections where the feature was present. Feature presence was binarized with a threshold of r*\u0026gt;0.4 for F1,2 and r*\u0026gt;0.5 for F3. Arbitrary thresholding was applied only for visualization purposes; un-thresholded matrices are presented in the SI \u003cstrong\u003eFig. S10\u003c/strong\u003e. From the matrix, three entries were randomly chosen, with graded R values from low to high, to demonstrate their associated evoked responses in the time domain. For plotting the evoked responses, the black line is the site-averaged signal across the evoked potentials stimulated/recorded in the same anatomical regions – it is a group-level average which may involve different subjects whose stimulated/recording sites are in the same brain region. The grey-shaded area around the line is the standard error over sites. (b) Delayed F3 connectivity matrices for bilateral cortical and thalamic subregions with the same visualization scheme as aforementioned. (c) and (e) show brain heatmaps of the contrasted feature representation among the thalamic divisions (antTH vs. pstTH, midth vs. pstTH, midTH vs. pstTH), respectively for the thalamic inflow (i.e., recording in the thalamus) and outflow (i.e., stimulation in the thalamus) pathways. All values were projected to one hemisphere for compact visualization. FS_LR brain surface space (with symmetric left and right hemisphere) is used to minimize the visualizing bias caused by interhemispheric anatomical differences (https://osf.io/k89fh/wiki/Surface/). Grey dots on the brain surface indicate electrode coverage, and the color radius around the dots indicate local R values of the given contact. Due to sparse recording, the color is also projected to the brain surface with a Gaussian function over distance from the source, to approximate a whole-brain level estimation. Coloring intensity on the brain surface has been adjusted for regional density of electrode coverage. The presented brain heatmaps are not thresholded; formal statistical testing results and model details can be found in \u003cstrong\u003eTable S2-5\u003c/strong\u003e.\u003c/p\u003e","description":"","filename":"6.png","url":"https://assets-eu.researchsquare.com/files/rs-4366486/v1/19f5ca04cda84ea1775ba633.png"},{"id":57379910,"identity":"c6cdb3ba-e389-4ad8-91a9-d837326d57d8","added_by":"auto","created_at":"2024-05-30 01:12:32","extension":"png","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":259982,"visible":true,"origin":"","legend":"\u003cp\u003eSee image above for figure legend.\u003c/p\u003e","description":"","filename":"7.png","url":"https://assets-eu.researchsquare.com/files/rs-4366486/v1/2f1882d42faa50d4e4067e8b.png"},{"id":86835933,"identity":"4bff6bec-d873-4ea4-9cc6-fc3d23a5b936","added_by":"auto","created_at":"2025-07-16 07:09:32","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":3719601,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-4366486/v1/a7aa9d22-e910-4596-9cba-5efb62e753c1.pdf"},{"id":57379911,"identity":"17dde2fe-69f8-4c20-a9f5-397449bea529","added_by":"auto","created_at":"2024-05-30 01:12:32","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":8813027,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryMaterials.docx","url":"https://assets-eu.researchsquare.com/files/rs-4366486/v1/f426af97018ea2d3445c8800.docx"}],"financialInterests":"There is \u003cb\u003eNO\u003c/b\u003e Competing Interest.","formattedTitle":"Causal Cortical and Thalamic Connections in the Human Brain","fulltext":[{"header":"Introduction","content":"\u003cp\u003eThe brain\u0026rsquo;s dynamics are shaped by electrophysiological interactions throughout its subregions, encompassing not only cortical but also subcortical regions. Advancements in neuroimaging have provided significant insights into the global architecture of the brain\u0026rsquo;s functional connectivity\u003csup\u003e1\u003c/sup\u003e. However, we know considerably less about the dynamic, fast-paced and causal electrophysiological relationships in the human brain. This knowledge gap is even more pronounced when it comes to understanding the role of subcortical areas such as the thalamus, which is known to play a key role in modulating the global dynamics of the brain\u003csup\u003e2\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eA classic method of studying causal electrophysiological relationships in the brain is by sending repeated single electrical pulses to intracranially implanted electrodes while recording the presence or absence of electrophysiological changes in all other areas of the brain where recording electrodes are present\u003csup\u003e3\u003c/sup\u003e. This method has been traditionally referred to as the study of cortico-cortical evoked potentials\u003csup\u003e4\u003c/sup\u003e (CCEPs) since it has primarily been used to study causal electrophysiological connectivity between pairs of \u003cem\u003ecortical\u003c/em\u003e regions. Similar studies of thalamocortical and corticothalamic connections have been rarely conducted since direct recordings from and stimulations in the human thalamus are extremely rare in human neuroscience research\u003csup\u003e5\u0026ndash;7\u003c/sup\u003e. As a result, it remains unknown if stimulation of the thalamus will have a different, or the same, effect on other regions of the brain. Moreover, with a few exceptions\u003csup\u003e7,8\u003c/sup\u003e, previous studies have been largely reliant on simple univariate measures detecting large peaks or the time-to-peak in the evoked signals, and limit responses to fixed windows of interest. These traditional methods are not able to capture the complex dynamics of physiological responses generated by electrical stimulations (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e); and they can vary significantly depending on the chosen cutoff z-scores for determining an effect, current intensity and the distance between stimulating and recording electrode contacts\u003csup\u003e9\u003c/sup\u003e. These problems become more significant when we lack prior data to inform our hypotheses about evoked potentials caused by human subcortical structures such as the thalamus.\u003c/p\u003e \u003cp\u003eIn the current study, we aim to use a novel approach to provide a map of causal electrophysiological connections in the human brain including the thalamus.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e"},{"header":"Results","content":"\u003cp\u003eOur results are based on intracranial recordings and stimulations in twenty-seven participants (40.7% female; mean age\u0026thinsp;\u0026plusmn;\u0026thinsp;SD: 34.9\u0026thinsp;\u0026plusmn;\u0026thinsp;10.0 years) diagnosed with focal epilepsy. Each participant provided data from an average of 180 (\u0026plusmn;\u0026thinsp;46) electrode sites. The entire dataset encompasses 4864 sites across all subjects, with each adjacent pair stimulated approximately 45 times (with sufficient inter-trial intervals to avoid overlapping effects), while recordings were taken from all other available sites.\u003c/p\u003e \u003cp\u003eWe leveraged a recently developed clinical method that involves The brain's functional architecture is intricately shaped by causal connections between its cortical and subcortical structures. Here, we studied 27 participants with 4864 electrodes implanted across the anterior, mediodorsal, and pulvinar thalamic regions, and the cortex. Using data from electrical stimulation procedures and a data-driven approach informed by neurophysiological standards, we dissociated three unique spectral patterns generated by the perturbation of a given brain area. Among these, a novel waveform emerged, marked by delayed-onset slow oscillations in both ipsilateral and contralateral cortices following thalamic stimulations, suggesting a mechanism by which a thalamic site can influence bilateral cortical activity. Moreover, cortical stimulations evoked earlier signals in the thalamus than in other connected cortical areas suggesting that the thalamus receives a copy of signals before they are exchanged across the cortex. Our causal connectivity data can be used to inform biologically-inspired computational models of the functional architecture of the brain.\u003c/p\u003e \u003cp\u003e \u003cem\u003ethalamic recordings\u003c/em\u003e in neurosurgical patients undergoing invasive stereotactic electroencephalography (sEEG)\u003csup\u003e5\u003c/sup\u003e. In this procedure, intracranial electrodes were implanted in three different sites of the thalamus \u0026mdash; corresponding to anterior, mediodorsal, and pulvinar nuclei \u0026mdash; by extending electrodes initially placed in the cortical areas. This procedure was designed for extended clinical indication without increasing the number of implanted electrodes.\u003c/p\u003e \u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eIdentifying Common Stable Patterns in Stimulation Evoked Potentials\u003c/h2\u003e \u003cp\u003eTo determine true physiological effects upon electrical stimulation, we used a non-linear manifold learning algorithm\u003csup\u003e10\u003c/sup\u003e integrated with human judgement, based on the spectral information of evoked potentials (EPs) across time and frequencies (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea-c). In this approach, we first labeled the data partially and tentatively, i.e., subject to subsequent data-driven corrections, based on two measures of connectivity: (i) trial-averaged \u003cem\u003epower\u003c/em\u003e of EPs as a measure of strength, and (ii) \u003cem\u003einter-trial phase coherence\u003c/em\u003e (ITPC) across repeated stimulations as a measure of consistency for the evoked responses\u003csup\u003e11\u003c/sup\u003e. We applied a semi-supervised machine learning approach to decipher the inherent data structure from this preliminary labelling, and re-classified the data, within each subject, into two separate clusters (\u003cb\u003eFig. S2\u003c/b\u003e): one cluster of data where stimulation of the seed site caused significant EPs in the target site, and another cluster of data where stimulation of the seed site did not cause any significant EPs in the target site. For simplicity, we called these two clusters \u0026ldquo;\u003cem\u003eactivated\u003c/em\u003e\u0026rdquo; and \u0026ldquo;\u003cem\u003enon-activated\u0026rdquo;\u003c/em\u003e, respectively.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e As a result, EPs in the activated cluster were differentiated from the spontaneous activity based on prominent changes in both power and ITPC relative to the baseline (pre-stimulation intrinsic neural activities), and also differentiated from the unnatural signals that we have manually labelled as noise (e.g., stimulation artifact or bad channels) by providing those typical cases. Once all subjects\u0026rsquo; EPs were successfully classified as \u0026ldquo;activated\u0026rdquo; vs. \u0026ldquo;non-activated\u0026rdquo;, we ensured the cluster\u0026rsquo;s meaningful representation of the data by manually inspecting the original and reconstructed data for every subject. With the machine-learning refined labelling (activated vs. non-activated) as a benchmark, the preliminary activation criteria were re-evaluated for further sanity checks (\u003cb\u003eTable S1\u003c/b\u003e).\u003c/p\u003e \u003cp\u003eWe next aimed to characterize the prominent electrophysiological properties of the EPs in the \u0026ldquo;activated\u0026rdquo; cluster. Dimensionality reduction with uniform manifold approximation and projection (UMAP) revealed that the inherent structure of the data was different for stimulation-recording pairs within the same hemisphere (i.e., ipsilateral) compared to those across the two hemispheres (i.e., contralateral connections). Surprisingly, we also discovered that the inherent structure of the data was different for EPs generated by stimulation of the thalamic sites vs. cortical sites (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ed; more label mapping can be found in \u003cb\u003eFig. S3\u003c/b\u003e). The thalamic vs. cortical recording sites, however, do not exhibit such obvious differences (\u003cb\u003eFig. S3c\u003c/b\u003e vs. \u003cb\u003eFig. S3d\u003c/b\u003e). We refer to this anatomical features (1. Stimulation from the cortex or thalamus, 2. Ipsilateral or contralateral connections) as UMAP localizers. We then searched for the significant features that set the data apart in the embedding space. To dissociate the data that contain the features of interest, we applied supervised learning with the UMAP localizers as labels. As expected, it resulted in four clusters in the embedding space with distinct spectral patterns, i.e., clusters of thalamic, cortical, ipsilateral, and contralateral EPs \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ee\u003cb\u003e)\u003c/b\u003e.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec4\" class=\"Section2\"\u003e \u003ch2\u003eDissociating Distinct Neural Features in Evoked Responses\u003c/h2\u003e \u003cp\u003eOnce the four clusters were identified, we used the power and ITPC spectrograms from each cluster for further statistical testing. Cluster-based permutation testing identified the significant time/frequency boundaries of three distinct neural features in the data that we labeled as Feature 1 (F1), Feature 2 (F2), and Feature 3 (F3) as tentatively shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ee by visual identification. Details of statistical testing procedures are presented in the \u003cem\u003eMethods\u003c/em\u003e. We summarize the three features as follows (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e3\u003c/span\u003ea, b): F1 is characterized by significant increases in power of high frequency activity in the gamma range and the strength of ITPC within the first 10-60ms. Of note, F1 was clearly dissociable from the artifact which was showed as a high-frequency (\u0026gt;\u0026thinsp;100 Hz) peak power immediately after the stimulation (\u0026lt;\u0026thinsp;10ms) and lacking signal continuation in both frequency and time dimensions (\u003cb\u003eFig. S5)\u003c/b\u003e. F2 is characterizable by an increased power in the theta to alpha range, without phase locking, happening\u0026thinsp;~\u0026thinsp;120 ms upon stimulation. Later than F2 (~\u0026thinsp;200 ms post-stimulation), F3 is featured by a prolonged power increase in the theta band, with a phase-locking feature which indicates oscillation (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e4\u003c/span\u003ea).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eWithin-subject comparisons showed that the ipsilateral cortical stimulations compared to contralateral cortical stimulations (i.e., COR-ipsi\u0026thinsp;\u0026gt;\u0026thinsp;COR-contr) caused significantly stronger F1 and F2 signals (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e3\u003c/span\u003eb). Similarly, ipsilateral thalamic stimulations compared to contralateral thalamic stimulations (i.e., THAL-ipsi\u0026thinsp;\u0026gt;\u0026thinsp;THAL-contr) caused stronger F1 and F2 signals (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e3\u003c/span\u003eb) i.e., stimulations of the thalamic and cortical sites lead to a stronger F1 and F2 signals in ipsilateral (compared to contralateral) brain sites. Of note, F1 was lacking in the contralateral recordings to the thalamic stimulation (i.e., THAL-contr).\u003c/p\u003e \u003cp\u003eImportantly, the EPs generated from thalamic stimulations were distinct from the cortical stimulations by the presence of a strong F3 signal (i.e., higher theta frequency change starting emerging around 165 ms and persisting for about 250 ms (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e3\u003c/span\u003eb). By checking the signals in the time domain, combining the evidence of the evoked power and ITPC, we confirmed (i) the presence of oscillatory activity embedded in the F3 signal, and (ii) that the F3 signal is not a delayed F2 or continuation of it. In the time domain, F3 indicates a\u0026thinsp;~\u0026thinsp;5 Hz oscillation, but F2 indicates a wider curve with jittered peaks across sites. In the time-frequency domain, both F2 and F3 share similar power spectrum, but their ITPC pattern differs significantly in that F2 does not have a phase-locking feature (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e4\u003c/span\u003ea). Moreover, unlike F2, F3 is lacking in cortical EPs, instead, it is specific to thalamic EPs with both THAL-ipsi and THAL-contr clusters showing this delayed and long-lasting F3 activity.\u003c/p\u003e \u003cp\u003eWe then investigated the spectrograms separately for different recording sites. Profiles of THAL\u0026agrave;COR, COR\u0026agrave;THAL and COR\u0026agrave;COR in both ipsilateral and contralateral pairs reflected the main results of the UMAP localizers, which showed that the stimulation sites, rather than the recording sites, dominate the spectral features.\u003c/p\u003e \u003cp\u003eNotably, multi-site thalamic coverage enabled us to examine intra- and inter-thalamic causal connectivity at the individual level (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e3\u003c/span\u003ec\u003cb\u003e)\u003c/b\u003e. The results provided unprecedent data from the thalamus where we could compare the connectivity of different thalamic subregions with the available brain recording sites and examine the effect of electrical stimulation of a given thalamus upon the activity of sites in other thalamic subregions within the same or across the other hemisphere (\u003cb\u003eFig. S7\u003c/b\u003e).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec5\" class=\"Section2\"\u003e \u003ch2\u003eTiming and Global Patterns of Connectivity Indicated by Feature Presence\u003c/h2\u003e \u003cp\u003eNext, we aimed to understand how each of the three neural features unfolds in time within individual instances of connectivity for the four categories of interest (i.e., COR-ipsi, COR-contr, THAL-ipsi and THAL-contr). The group-averaged spectral information within the significant time windows from the COR-ipsi, COR-contr and THAL-contr categories was used to create the feature template for F1, F2 and F3, respectively (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e4\u003c/span\u003ea). Using a sliding-window correlation approach (see details in \u003cem\u003eMethods and\u003c/em\u003e illustration in Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e4\u003c/span\u003eb), we estimated the strength of feature representations (r*) and their presence in time.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eCorresponding to the previous group-level cluster-based permutation testing on the spectrograms, the decoding results also showed that both F1 and F2 peaks were seen more strongly for ipsilateral than contralateral EPs (Dr* = 0.11, t\u0026thinsp;=\u0026thinsp;48.05, p\u003csub\u003eFDR-corr\u003c/sub\u003e\u0026thinsp;\u0026lt;\u0026thinsp;\u0026lt;\u0026thinsp;0.01, n-connection [conn.]\u0026thinsp;=\u0026thinsp;57881, n-subject [sbj.]\u0026thinsp;=\u0026thinsp;26 for F1, and Dr* = 0.09, t\u0026thinsp;=\u0026thinsp;44.90, p\u003csub\u003eFDR-corr\u003c/sub\u003e\u0026thinsp;\u0026lt;\u0026thinsp;\u0026lt;\u0026thinsp;0.01, n-conn. = 57881, n-sbj. = 26 for F2, using the stimulation category \u0026ldquo;COR/THAL\u0026rdquo; as a covariate, corrected for the nested individual, regional and site effects in a hierarchical linear model [HLM]). Additionally, the decoding method was able to provide more precise time information than the previous test, which showed that both F1 and F2 were detected earlier in the EPs responding to cortical versus thalamic stimulation (Dlatency\u0026thinsp;=\u0026thinsp;5.60 ms, t\u0026thinsp;=\u0026thinsp;27.20, p\u003csub\u003eFDR-corr\u003c/sub\u003e\u0026thinsp;\u0026lt;\u0026thinsp;\u0026lt;\u0026thinsp;0.01, n-conn. = 37444, n-sbj. = 26 for robust F1 [r*\u0026gt;0.4, latency in the 10\u0026ndash;100 ms range], and Dlatency\u0026thinsp;=\u0026thinsp;8.18 ms, t\u0026thinsp;=\u0026thinsp;22.49, p\u003csub\u003eFDR-corr\u003c/sub\u003e\u0026thinsp;\u0026lt;\u0026thinsp;\u0026lt;\u0026thinsp;0.01, n-conn. = 48300, n-sbj. = 26 for robust F2 [r*\u0026gt;0.4, latency in the 70\u0026ndash;200 ms range]; using the cross-hemispheric category \u0026ldquo;ipsi/contr\u0026rdquo; as a covariate in an HLM).\u003c/p\u003e \u003cp\u003eThe decoding results also echoed the previous finding that F3 was seen mainly in evoked potentials pertaining to the thalamic stimulations, with F3 peaking higher in ipsilateral than in contralateral EPs (Dr* = 0.11, t\u0026thinsp;=\u0026thinsp;48.05, p\u003csub\u003eFDR-corr\u003c/sub\u003e = 0.02, n-conn. = 6497, n-sbj. = 25), but without significant difference in the latency between ipsilateral and contralateral EPs (Dr* = 5.03 ms, t\u0026thinsp;=\u0026thinsp;1.86, p\u0026thinsp;=\u0026thinsp;0.06, n-conn. = 3106, n-sbj. = 25 for robust F3 with r* \u0026gt; 0.5, latency\u0026thinsp;\u0026gt;\u0026thinsp;200 ms). Additional post-hoc comparisons among sub-groups are presented in Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e5\u003c/span\u003e and SI \u003cb\u003eTable S9\u003c/b\u003e. The correspondence of the decoding results with the previous group-level spectrogram analysis suggests that our devised feature indices can reliably indicate the presence of the multivariate features with single values.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eWe further examined the whole-brain spatiotemporal patterns indicated by the three neural features. To do so, we constructed connectivity matrix with each cell filled with the connectivity index of R for indicating the signal strength between two brain regions (R = \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\stackrel{-}{{r}^{*}}\\)\u003c/span\u003e\u003c/span\u003e, averaged across the sites where the feature is present in the given regions). By organizing brain areas in the order of anatomical proximity, interesting patterns emerged: (i) the F1 matrix demonstrated a profile of modular connectivity, i.e., higher F1 representation was seen for connections between pairs located within the same hemisphere or lobe and within the ipsilateral thalamus (Fig.\u0026nbsp;\u003cspan refid=\"Fig14\" class=\"InternalRef\"\u003e6\u003c/span\u003ea, b); (ii) the F2 matrix revealed wide-spread representations across regions and hemispheres and demonstrated that the global architecture of connections based on F2 did not differ much from the F1 matrix. This seems to suggest that F2 is a combination of jittered rebounds following the same paths for the F1 signals (\u003cb\u003eFig. S10\u003c/b\u003e); (iii) unlike the F1 and F2 matrices which have a modular architecture, the F3 was sparsely present among the cortical connections, but largely present in the cortical responses (i.e., widespread and bilateral) associated with thalamic stimulation (Fig.\u0026nbsp;\u003cspan refid=\"Fig14\" class=\"InternalRef\"\u003e6\u003c/span\u003ef, g).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eIn a closer examination, we found that F1 corticothalamic and F3 thalamocortical connectivity features had the following unique characteristics in common: (i) F1 responses recorded in the thalamus had widespread and bilateral origins across the brain (i.e., stimulation of many cortical regions bilaterally could evoke responses in the thalamus) (Fig.\u0026nbsp;\u003cspan refid=\"Fig14\" class=\"InternalRef\"\u003e6\u003c/span\u003ea), and vice versa for F3, the stimulation of a thalamic site would cause widespread cortical areas bilaterally to show delayed slow oscillations (Fig.\u0026nbsp;\u003cspan refid=\"Fig14\" class=\"InternalRef\"\u003e6\u003c/span\u003eb); and (ii) both corticothalamic F1 and thalamocortical F3 responses were \u003cem\u003estronger\u003c/em\u003e and \u003cem\u003eearlier\u003c/em\u003e than corticocortical F1 and corticocortical F3 responses, respectively. The statistical details are presented as below.\u003c/p\u003e \u003cp\u003eWe compared the cortical versus thalamic targets responding to the same cortical seed of stimulation in the same participant, and confirmed the F1 responses in the thalamus were \u003cem\u003eearlier\u003c/em\u003e and \u003cem\u003estronger\u003c/em\u003e (Dlatency\u0026thinsp;=\u0026thinsp;2.79 ms, t\u0026thinsp;=\u0026thinsp;12.74, p\u003csub\u003eFDR-corr\u003c/sub\u003e\u0026thinsp;\u0026lt;\u0026thinsp;\u0026lt;\u0026thinsp;0.01, n-conn. = 48446, n-sbj. = 26; Dr* = 0.06, t\u0026thinsp;=\u0026thinsp;18.20, p\u003csub\u003eFDR-corr\u003c/sub\u003e\u0026thinsp;\u0026lt;\u0026thinsp;\u0026lt;\u0026thinsp;0.01, n-conn. = 30965, n-sbj. = 26, with HLM corrected for the covariance effect of cross-hemispheric connections and the nested random effect of 1378 different stimulation sites of different brain areas in all subjects) (Fig.\u0026nbsp;\u003cspan refid=\"Fig14\" class=\"InternalRef\"\u003e6\u003c/span\u003ea for R matrix, SI \u003cb\u003eFig. S10\u003c/b\u003e for latency matrix). In other words, stimulation of a given cortical site changed the activity of the thalamus before it changed the activity of its cortical targets (Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e5\u003c/span\u003ea and SI \u003cb\u003eTable S9\u003c/b\u003e for post-hoc latency comparisons).\u003c/p\u003e \u003cp\u003eMoreover, the F3 matrix revealed that the stimulation of thalamic sites caused delayed low-frequency oscillatory effects (i.e., F3) in a large mantle of the cortex in not only ipsilateral but also contralateral hemispheres (Fig.\u0026nbsp;\u003cspan refid=\"Fig14\" class=\"InternalRef\"\u003e6\u003c/span\u003eb). By contrast, the F3 feature was only sparsely present with the stimulation of non-thalamic seeds. Importantly, the latency of F3 responses caused by thalamic stimulations was significantly earlier than the non-thalamic seeds (36.27 ms earlier, t\u0026thinsp;=\u0026thinsp;21.60, p\u003csub\u003eFDR-corr\u003c/sub\u003e\u0026thinsp;\u0026lt;\u0026thinsp;\u0026lt;\u0026thinsp;0.01, n-conn. = 20823, corrected for the covariance of cross-hemispheric effect, and the nested random effect of 3041 recording sites of different areas in all subjects with an HLM) (Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e5\u003c/span\u003ec).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec6\" class=\"Section2\"\u003e \u003ch2\u003eDifferential connectivity among thalamic subregions\u003c/h2\u003e \u003cp\u003eUp to this point, we have remained agnostic to the specific anatomical location of the thalamic sites where evoked responses were recorded or stimulations were seeded. However, each subject had more than one electrode in the thalamus targeting the anterior, mediodorsal and pulvinar nuceli\u003csup\u003e5\u003c/sup\u003e. Given that the precise anatomical boundaries of each thalamic nucleus within each subject may be variable and difficult to ascertain with neuroimaging data, and that the field of evoked responses or electrical stimulations may cross the nuclear boundaries, we refrained from labeling our thalamic sites according to specific nuclei, and instead, referred to them as anterior (antTH), middle (midTH), and posterior thalamus (pstTH).\u003c/p\u003e \u003cp\u003eLeveraging the anatomical information of recordings and simulation sites - at the individual level - we compared the connectivity profiles of the three different thalamic sites. This analysis revealed unique site-specific profiles of connectivity between the thalamus and the cerebral cortex (Fig.\u0026nbsp;\u003cspan refid=\"Fig14\" class=\"InternalRef\"\u003e6\u003c/span\u003ec, d, e\u003cb\u003e)\u003c/b\u003e. As describing the details of differential connectivity of thalamic sites are outside the scope of this report, we only highlight some of the most salient findings. The full results of pair-wise comparisons among the thalamic subregions for their bi-directional connections with each of the brain regions are provided in the supplementary \u003cb\u003eTable S2-7\u003c/b\u003e.\u003c/p\u003e \u003cp\u003eNext, we compared the delayed thalamocortical F3 signal representation map with the ones for early thalamocortical F1 and early corticothalamic F1 maps, with pair-wise comparisons within the same pair of thalamic vs. non-thalamic sites. Using the measure of normalized feature representation \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\stackrel{-}{zscore\\left({r}^{*}\\right)}\\)\u003c/span\u003e\u003c/span\u003e, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig15\" class=\"InternalRef\"\u003e7\u003c/span\u003e and detailed in \u003cb\u003eTables S6 and S7\u003c/b\u003e, between a subregion of the thalamus and a given brain region of interest, we found significant asymmetry between thalamocortical F3 representations and thalamocortical F1 representations. Many regions had stronger F3 than F1 representations suggesting that the delayed thalamic (F3) outflow is represented across a more widespread cortical areas compared to the fast thalamic (F1) outflow signal (see brain regions with more brown colors in Fig.\u0026nbsp;\u003cspan refid=\"Fig15\" class=\"InternalRef\"\u003e7\u003c/span\u003ea and positive estimate numbers in \u003cb\u003eTable S6\u003c/b\u003e). The same applies to the comparison between thalamocortical F3 and corticothalamic F1 representations. In this comparison, the hippocampus (HPC) was an exception among all other examined brain regions: representation of F1 signals evoked by HPC and recorded in anterior and posterior thalamic sites was stronger than the representation of F3 signals recorded in HPC and evoked by the thalamus, with pair-wise comparisons for the same sites in the regions of interest (Dz\u0026thinsp;=\u0026thinsp;0.30, t\u0026thinsp;=\u0026thinsp;2.51, p\u003csub\u003eFDR-corr\u003c/sub\u003e =0.04, n-conn. = 190, n-sbj. = 21 for antTH; Dz\u0026thinsp;=\u0026thinsp;0.62, t\u0026thinsp;=\u0026thinsp;4.72, p\u003csub\u003eFDR-corr\u003c/sub\u003e\u0026thinsp;\u0026lt;\u0026thinsp;\u0026lt;\u0026thinsp;0.01, n-conn. = 194, n-sbj. = 18 for pstTH) (\u003cb\u003eTable S7\u003c/b\u003e).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e"},{"header":"Discussion","content":"\u003cp\u003eIn our study, we addressed the limitations of past studies by leveraging a novel approach to characterize the complex features of electrophysiological responses that are evoked by not only cortical but also thalamic stimulations in the human brain. This is important because, as we documented in Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e, reliance on the magnitude of recorded neural responses or expecting typical N1 and N2 profiles in the evoked responses may greatly bias our views of causal connections in the brain. The magnitude of evoked responses can be affected by several confounding factors such as distance from the source of stimulation, volume conduction, and electrode impedance\u003csup\u003e12\u003c/sup\u003e and the time-to-peak may be elusive when there is no clear N1 or N2 peaks. While typical N1 and N2 waveforms reported in the CCEP literature are abundant, it is important to note that the past reports were based on group level analysis of data from cortical stimulations\u003csup\u003e13,14\u003c/sup\u003e. As noted in Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e, individual evoked signals may be far more complex than simple N1 and N2 waveforms and the magnitude of evoked responses may not always pass an arbitrary threshold we set for significance. This issue is even more pressing when waveforms from subcortical structures such as the thalamus are included. Therefore, in line with recent attempts by other investigators\u003csup\u003e7,8\u003c/sup\u003e, we felt it was imperative to employ a multivariate and unbiased approach for sorting the distinct profiles of connectivity across cortical and thalamic sites.\u003c/p\u003e \u003cp\u003eIn our encoding method, we relied on the machine learning algorithms of UMAP to identify the inherent variabilities among the datasets. We emphasize that our data-driven approach was not blind to biologically plausible features. Instead, it was directly guided by the 2/3 of the data that had already been labelled based on established neurophysiological criteria for significant evoked responses (i.e., significant amplitude above the noise level, significant inter-trial phase coherence concerning reliability of responses across trials, and more importantly, labeling what constitutes a stimulation artifact). Our post-hoc analysis also confirmed that the machine-learning refined labeling largely corresponded to the activation detection results based on the preset biological criteria, among which, the criterion based on significant ITPC (i.e., consistency in the evoked responses over trials) had the highest correspondence (0.80). In contrast, the traditional activation detection criterion of large peak detection (\u0026gt;\u0026thinsp;5 z-score) had a high precision (0.99) but low sensitivity (0.53), which is in line with our observation: many thalamic evoked responses (especially the contralaterally recorded ones) would be missed if based solely on the amplitude criterion as those responses usually have lower amplitudes on the voltage traces, though not necessarily possessing less robust evoked features.\u003c/p\u003e \u003cp\u003eOur new approach yielded important results suggesting that the complex electrophysiological signals evoked by repeated single pulse electrical stimulations have three distinct features (i.e., F1, F2, and F3). In the context of past literature, unitary F1 and F2 features may correspond to clear N1 and N2 peaks, but we emphasize that changes in the neural activity of a given brain site evoked by the stimulation of another site cannot always be depicted as a uniform waveform consisting of N1 and N2, but rather, a complex signal with different components of the evoked waveform which may in turn reflect distinct physiological substrates with distinct sources. While we postulate that F1 and F2 may be related to the well-known N1 and N2 components of the EP signal described in the literature \u003csup\u003e12,14,15\u003c/sup\u003e, our data suggests that F3 is a unique feature of its own and not a prolonged F2 \u0026ndash; even though F2 and F3 features have seemingly similar power spectrogram peaking at a low frequency range. Our rationale is as follows: (i) F3 has a peak oscillation at about 5 Hz, but F2 has wider frequency range with jittered peaks across trials. (ii) F3\u0026rsquo;s ITPC spectrogram differs significantly from the F2\u0026rsquo;s, as F3 has a phase-locking feature at the theta frequency but F2 does not (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e4\u003c/span\u003ea). The oscillatory feature of F3 was also observed at the same time-frequency cluster as long-lasting theta oscillations (\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;2 cycles). The oscillatory feature of F3, unlike F2, was observed in the time domain: it was strongly time-locked and exhibited at least two cycles of oscillation - even when averaging the EPs recorded across all recording areas and across all thalamic stimulation sites. (iii) The temporal window of significant time-frequency clusters in F2 and F3 do not overlap (i.e., F3 always happens after F2, peaking\u0026thinsp;\u0026gt;\u0026thinsp;200 ms post-stimulation). (iv) Lastly, F3 feature is present specifically when the thalamus is stimulated (including both THAL-ipsi and THAL-contr clusters), while F2 does not have this anatomical specificity.\u003c/p\u003e \u003cp\u003eWe documented that the F3 feature was detected with the stimulation of \u003cem\u003eall three thalamic sites\u003c/em\u003e. This suggests that the F3 signal is not produced by the stimulation of only one specific thalamic nucleus. However, it's important to acknowledge a limitation of our study, which stemmed from the absence of data from a number of other thalamic nuclei. For instance, we have no information about the profile of electrophysiological responses evoked by the stimulation of thalamic nuclei such as the geniculate bodies. This limitation was primarily due to the clinical constraints and ethical considerations inherent in studies involving human subjects.\u003c/p\u003e \u003cp\u003eWe also note that the F3 oscillatory evoked responses were also seen with the stimulation of a few non-thalamic brain regions - especially the hippocampus and amygdala \u0026ndash; both of which are non-cortical structures [Figure \u003cspan refid=\"Fig14\" class=\"InternalRef\"\u003e6\u003c/span\u003eb]). interestingly, the latency of F3 signal evoked by the stimulation of these non-thalamic sites was significantly delayed compared to the latency of F3 signal evoked by the stimulation of the thalamic sites, and they were mostly present in ipsilateral recording sites while the thalamic stimulations evoked F3 in far more widespread areas involving bilateral hemispheres.\u003c/p\u003e \u003cp\u003eAs the thalamus is considered to play a key role in cortical rhythms\u003csup\u003e16\u003c/sup\u003e, and as the hippocampus is known for its role in generating theta rhythms\u003csup\u003e17\u003c/sup\u003e, we hope our findings will motivate future systematic studies to determine how theta oscillations observed throughout a broad surface of higher association areas in the cerebral cortex during cognitive activities tasks or rest\u003csup\u003e18,19\u003c/sup\u003e are related to the delayed modulatory effects exerted by the thalamic (or medial temporal lobe) structures. It also remains to be determined if the bilateral presence of the F3 signal evoked by thalamic stimulations is functionally related to, or is crucial for, the integration of information across the two hemispheres.\u003c/p\u003e \u003cp\u003eOur approach with signal decoding of the stimulation and recording data revealed that stimulation of a given cortical region evoked \u003cem\u003eearlier\u003c/em\u003e and \u003cem\u003estronger\u003c/em\u003e responses in the thalamus compared to the cortical targets. A hypothetical interpretation of this finding is that the thalamus acts as a universal receiver of information being exchanged across cortical regions (i.e., a \u0026ldquo;listener of corticocortical dialogue\u0026rdquo;) \u0026ndash; which is a necessary for it to function as a key structure for information integration as detailed recently\u003csup\u003e2\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eAdditionally, the multisite thalamic stimulation and recording within the same individual offered us an unprecedent opportunity to examine the intrathalamic and interthalamic causal connections. We demonstrate that the stimulation of a given thalamic site led to fast high frequency modulation of activity in other thalamic sites ipsilaterally as well as delayed low frequency modulation of thalamic sites contralaterally. These findings provide a plausible electrophysiological mechanism for intra-thalamic and inter-thalamic connectivity, which may have important implications and clinical relevance: Stimulating a specific thalamic nucleus with DBS (e.g., for treatment of epilepsy\u003csup\u003e20\u003c/sup\u003e) is likely to impact the activity of other thalamic nuclei (beyond the one directly targeted) and the cortical networks they are associated with. Thus, the therapeutic benefits of DBS might be related to a broader network of brain regions that are modulated rather than a focal network targeted. However, this does not imply that the effect of stimulation of each targeted thalamic site is diffuse and all-encompassing, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e5\u003c/span\u003e, stimulations in three subregions of the thalamus do not elicit anatomically identical responses, but instead show relative anatomical preference.\u003c/p\u003e \u003cp\u003eWe are mindful that our data was acquired from patients with a neurological disease (epilepsy) and one may ask to what extent they are generalizable to normal human brains. While we acknowledge that molecular and cellular changes accompany local circuit rewiring in patients with chronic epilepsy\u003csup\u003e21\u003c/sup\u003e, there is no firm evidence suggesting changes in the fundamental topology of connections at the system level in these patients\u003csup\u003e22\u003c/sup\u003e. In the extant literature, unless studying the circuit level and local changes in a region-of-interest type of analysis\u003csup\u003e23\u0026ndash;25\u003c/sup\u003e or graph theoretic measures for deriving high-level estimation of network capacity\u003csup\u003e26,27\u003c/sup\u003e, the global network topology of the brain has been shown not to be significantly different between epileptic cohort and the normal control\u003csup\u003e28\u003c/sup\u003e. Several studies focusing on epileptic brain reorganization have shown that the same sets of canonical brain networks can be discovered by the data-driven methods, just as in normal controls\u003csup\u003e29\u0026ndash;31\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eWe are also mindful that electrode coverage in our subjects was limited and sparse. This is unfortunately a limitation of the intracranial approach in human participants that we, as researchers, have little control over since the implantation of electrodes in all participants was clinically driven. We made the best effort to reduce the confounding effects by maximizing the number of subjects, and diversifying the cohort with different epilepsies so as to avoid a system bias of encountering pathological tissues concentrated in specific areas.\u003c/p\u003e \u003cp\u003eDespite the limitations of our study, we hope that our findings can serve as the steppingstone towards a wider understanding of the architectural topology of interactions across the cerebral cortex and between those and subcortical structures such as the thalamus.\u003c/p\u003e \u003cp\u003eIn closing, we acknowledge that the causal connectivity maps presented in our current study provide only one aspect of the functional architecture of the brain, and as such, they serve as rough approximations for the dynamic causal interactions occurring in real-life scenarios and during cognitive processing and are complimentary to those connectivity measures generated by structural diffusion MRIs or resting state or task specific imaging studies\u003csup\u003e3,4,32\u0026ndash;34\u003c/sup\u003e. To highlight the functional significance of RSEPS-based connectivity maps, we recently confirmed that the strength of RSEPS connectivity between sites correlates with the strength of their co-activations during a cognitive task, i.e., if a specific pairs of neuronal populations across the cortex and thalamus had high RSEP-based connectivity strength they had higher correlation of high-gamma activations during a memory task\u003csup\u003e11\u003c/sup\u003e. Additionally, it has been shown that seizures propagate from cortical seizure onset zones to thalamic sites\u003csup\u003e5\u003c/sup\u003e or other cortical areas\u003csup\u003e35\u003c/sup\u003e that show faster and stronger RSEPS connectivity measures with the seizure onset zone.\u003c/p\u003e"},{"header":"Online Methods","content":"\u003cdiv id=\"Sec9\" class=\"Section2\"\u003e\n \u003ch2\u003eResource Availability\u003c/h2\u003e\n \u003cp\u003e\u003cstrong\u003eMaterials availability\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003eRaw data, including electrophysiology collected in this study will be shared publicly after publication.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec10\" class=\"Section2\"\u003e\n \u003ch2\u003eData and code availability:\u003c/h2\u003e\n \u003cul\u003e\n \u003cli\u003e\n \u003cp\u003eThe electrophysiological data have been deposited at Mendeley and are publicly available as of the date of publication. The DOI will be listed in the key resources table.\u003c/p\u003e\n \u003c/li\u003e\n \u003cli\u003e\n \u003cp\u003eAll original codes have been deposited at Zenodo and is publicly available as of the date of publication. DOIs will be listed in the key resources table.\u003c/p\u003e\n \u003c/li\u003e\n \u003cli\u003e\n \u003cp\u003eAny additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.\u003c/p\u003e\n \u003c/li\u003e\n \u003c/ul\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec11\" class=\"Section2\"\u003e\n \u003ch2\u003eParticipants\u003c/h2\u003e\n \u003cp\u003eIn this study, twenty-seven participants (40.7% female; mean age\u0026thinsp;\u0026plusmn;\u0026thinsp;SD: 34.9\u0026thinsp;\u0026plusmn;\u0026thinsp;10.0 years) diagnosed with focal epilepsy were recruited (see Table \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003e). Each participant had 180\u0026thinsp;\u0026plusmn;\u0026thinsp;46 (mean\u0026thinsp;\u0026plusmn;\u0026thinsp;standard deviation) SEEG contacts implanted. The total number of electrode contacts being studied is 4864. All participants underwent invasive electrophysiological monitoring at our medical center as part of their treatment for refractory epilepsy. The placement of electrodes was determined exclusively based on clinical requirements. Prior to their involvement, all participants provided written informed consent, and the study protocol was approved by the Institutional Review Board for human experimentation.\u0026nbsp;\u003c/p\u003e\n \u003ctable id=\"Tab1\" border=\"1\"\u003e\n \u003ccaption language=\"En\"\u003e\n \u003cdiv class=\"CaptionNumber\"\u003eTable 1\u003c/div\u003e\n \u003cdiv class=\"CaptionContent\"\u003e\n \u003cp\u003eParticipant demographics\u003c/p\u003e\n \u003c/div\u003e\n \u003c/caption\u003e\n \u003cthead\u003e\n \u003ctr\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eParticipants\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eAge\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eGender\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eRace\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eHemisphere coverage\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eNumber of electrodes\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eSeizure onset zone\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eEpilepsy (year)\u003c/p\u003e\n \u003c/th\u003e\n \u003c/tr\u003e\n \u003c/thead\u003e\n \u003ctbody\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS01_166\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e38\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eF\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eWhite\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e176\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eright hippocampus\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e2\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS02_167\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e23\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eF\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eMixed\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eR\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e168\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eright posterior insula\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e4\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS03_169\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e23\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eM\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eWhite\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e186\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003einferior occipital gyrus\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e19\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS04_170\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e40\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eM\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eAsian\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eL\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e116\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eleft posterior inferior temporal lobe\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e24\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS05_171\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e52\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eM\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eAsian\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e252\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eleft lateral temporal and right orbitofrontal region\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e35\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS06_172\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e47\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eM\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eAsian\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e174\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003emesial temporal and insula\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e14\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS07_176\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e35\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eF\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eWhite\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eL\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e156\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eleft hippocampus\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e2\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS08_177\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e19\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eM\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eWhite\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e168\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003ehippocampus\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e4\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS09_178\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e28\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eM\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eWhite\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e252\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eright anterior temporal lobe\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e9\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS10_181\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e23\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eF\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eWhite\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e248\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eright occipital lobe\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e5\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS11_182\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e37\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eF\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eWhite\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e260\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003ehippocampus\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e27\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS12_183\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e33\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eM\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eWhite\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e174\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eright anterior temporal lobe\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e5\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS13_185\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e36\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eF\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eWhite\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e142\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eleft mesial temporal lobe\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e26\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS14_188\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e31\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eM\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eWhite\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e122\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003ehippocampus\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e10\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS15_189\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e27\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eM\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eHispanic\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e184\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eleft posterior insula and left hippocampus\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e19\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS16_190\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e50\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eF\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eWhite\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eR\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e116\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eright primary somatosensory cortex\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e4\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS17_192\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e20\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eM\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eHispanic\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eL\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e88\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eleft hippocampus\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e5\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS18_193\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e40\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eM\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eHispanic\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e182\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eright periventricular nodular heterotopia, right occipital region\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e30\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS19_194\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e37\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eM\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eHispanic\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e256\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eleft medial hippocampus, right posterior insula and pulvinar\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e33\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS20_195\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e47\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eF\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eBlack\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e138\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003emesial temporal lobe and insula\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e17\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS21_196\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e40\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eF\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eWhite\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e158\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eleft orbital frontal and left mesial temporal lobe\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e20\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS22_197\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e20\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eM\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eBlack\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e214\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eright supramarginal area\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e15\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS23_198\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e38\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eM\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eBlack\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e204\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003emesial temporal lobe\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e4\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS24_199\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e36\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eM\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eWhite\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e178\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eright superior temporal gyrus\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e20\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS25_201\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e27\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eM\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eHispanic\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e174\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003emesial temporal lobe\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e7\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS26_202\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e48\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eF\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eAsian\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e194\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eleft mesial occipital area\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e16\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS27_205\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e48\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eF\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eWhite\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eB\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e184\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eright anterior hippocampus\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e8\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tbody\u003e\n \u003c/table\u003e\n \u003cp\u003eM, male; F, female; R, right; L, left; B, both.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec12\" class=\"Section2\"\u003e\n \u003ch2\u003ePatient Safety \u0026amp; Ethics\u003c/h2\u003e\n \u003cp\u003eAs noted, we did not implant extra electrodes for thalamic recordings. We only extended the electrodes that were clinically planned for cortical monitoring to reach the thalamus. Our procedure was based on patients\u0026rsquo; informed consent and all our research related activities were IRB approved. To minimize injury, we used a reduced diameter obturating stylet and reduced diameter electrodes with 0.86mm diameter (AdTech, Inc). With this surgical approach, we have observed no complications, no thalamic hemorrhage or edema, and no neurological symptoms post operatively in any of the patients examined. We note that the thalamic recording in epilepsy patients has been widely practiced in some centers in Europe for two decades yielding clinically important information\u003csup\u003e36\u0026ndash;40\u003c/sup\u003e. The method has recently gained ground in the U.S.\u003csup\u003e41\u0026ndash;44\u003c/sup\u003e, and a recent neurosurgery editorial recommended and encouraged more centers in the U.S. to do the same\u003csup\u003e45\u003c/sup\u003e. Additionally, iEEG implantation strategy is decided by a group of clinicians who meet weekly to decide the clinical needs of the patient, and not by this research project. Also, as noted, all patients provide informed consent prior to being enrolled in any research studies.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec13\" class=\"Section2\"\u003e\n \u003ch2\u003eElectrode Implantation\u003c/h2\u003e\n \u003cp\u003eMRI sequences for imaging the thalamus and its subnuclei included BRAVO, fGATIR and MP2RAGE. We co-registered the post-implant CT with a pre-implant MRI. High resolution T1, fast gray matter acquisition T1 inversion recovery (FGATIR), and T1 post-contrast imaging were used for planning. Trajectories were planned to traverse in an orthogonal plane to capture the cortical frontal or temporal operculum, insula, and to be extended into specific sites of interest in the thalamus. The priority for all trajectories was safety and avoidance of blood vessels. To achieve implementation of our multisite sampling approach, the cortical electrode trajectories were extended to reach three different thalamic sites: anterior, mid, and posterior thalamus as detailed in our recent publication\u003csup\u003e5\u003c/sup\u003e. Approximate locations and number of electrodes along with their trajectories were all planned in a multidisciplinary surgical epilepsy conference with detailed review of presurgical data leading to a clinical hypothesis of the most likely seizure onset zones. Electrode placements were at the discretion of the clinical care team, who were not specifically involved with this research project.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec14\" class=\"Section2\"\u003e\n \u003ch2\u003eCo-localization of electrodes\u003c/h2\u003e\n \u003cp\u003eA thin cut head-CT, obtained after electrode implantation, was co-registered to pre-operative 3Tesla MRI data for verification of the trajectory. Precise electrode positioning in the FreeSurfer surface space, voxel space and MNI space were automatically extracted by the iElVis toolbox\u003csup\u003e46\u003c/sup\u003e. The T1-weighted MRI scan was used to generate 3D cortical volume and subcortical segmentation using recon-all command of Freesurfer v7.3.2\u003csup\u003e47\u003c/sup\u003e. The post-implant CT scan were co-registered to the pre-implant MRI using the \u003cem\u003eflirt\u003c/em\u003e from the Oxford Centre for Functional MRI of the Brain Software Library\u003csup\u003e48,49\u003c/sup\u003e or using \u003cem\u003ebbregister\u003c/em\u003e from Freesurfer\u003csup\u003e50\u003c/sup\u003e to get the best results. We manually labelled each electrode on the T1-registered CT image using BioImageSuite\u003csup\u003e51\u003c/sup\u003e. The electrode coordinates in the native anatomical space were inspected and manually labeled by the experienced neurologist J.P., based on the individual brain\u0026rsquo;s morphology and landmarks.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec15\" class=\"Section2\"\u003e\n \u003ch2\u003eThalamic parcellation\u003c/h2\u003e\n \u003cp\u003eTo check the concordance between manual and automatic labeling of thalamic electrode sites, individual contact centers of electrode masses were defined in native T1 space for each subject. The center of mass of each contact was then converted into an X,Y,Z coordinate in the MNI space. A recently developed and widely used MNI thalamic atlas\u003csup\u003e52\u003c/sup\u003e was utilized to parcellate the nuclear location of each contact within the thalamus. For this, a 1mm cubic voxel region was created around each contact center of mass (contact neighborhood). For each contact neighborhood, the fraction of voxels that overlapped with each thalamic nucleus in the THOMAS atlas was calculated. Contact neighborhoods that fell solely within a single THOMAS nuclear mask would have a value of 1 for that specific nucleus, while contact neighborhoods that had no overlapping voxels with a given nuclear mask would have a value assigned to 0 for that specific nucleus. Thus, for every contact neighborhood anatomically inside the thalamus, the fractional overlap with each THOMAS nuclear mask was calculable. For each site, this within-subject fractional overlap was summed across subjects and across all thalamic contacts to generate overall contact neighborhood fractional overlap values.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec16\" class=\"Section2\"\u003e\n \u003ch2\u003eIdentifying Regions of Interest\u003c/h2\u003e\n \u003cp\u003eAll electrodes contact locations were manually labelled with brain regions of interest (Fig.\u0026nbsp;\u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003ea) by the experienced neurologist J.P. on the team. This process involves carefully inspecting the CT-MRI fused images in axial, sagittal, and coronal planes for the implanted electrodes, and their surrounding anatomical landmarks in the brain morphology (e.g., basal ganglia including putamen, and pallidum etc.). Contacts entirely or partially located in the white matter (e.g., internal capsule, or cingulum etc.) were labelled as white matter. We emphasize that these regions of interest were not identified in the standard atlas space, but in each participant\u0026rsquo;s own high-resolution T1 scan to ensure anatomical precision of electrode localization.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec17\" class=\"Section2\"\u003e\n \u003ch2\u003eRepeated single electrical pulse stimulation (RSEPS)\u003c/h2\u003e\n \u003cp\u003eStimulations were performed in a bi-polar manner stimulating between a pair of electrodes. Each RSEPS trial entailed an instantaneous (pulse duration\u0026thinsp;=\u0026thinsp;0.2 ms) 6 mA biphasic square-wave pulse. For electrode contacts near the seizure onset zones (identified by the clinical team), we stimulated at 4 mA. Repeated pulses (n\u0026thinsp;=\u0026thinsp;42\u0026thinsp;\u0026plusmn;\u0026thinsp;2 trials) were delivered in a row with an interval of 2 seconds. We excluded data from (i) bipolar channels with both contacts located in the white matter (see above for details of anatomical localization), (ii) bipolar contacts not within the same anatomical region (e.g., one contact in insula and the other in orbitofrontal cortex), (iii) recording bipolar sites too close to the stimulation site (Euclidean distance between their middle points\u0026thinsp;\u0026lt;\u0026thinsp;=\u0026thinsp;5 mm), and (iv) channels largely contaminated by stimulation artifact, i.e., recording large signals (\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;7 SD of the baseline) only in the artifact window (\u0026lt;\u0026thinsp;10 ms), and (v) bad channels based on more specific criterion as detailed below.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec18\" class=\"Section2\"\u003e\n \u003ch2\u003eIntracranial EEG data acquisition and preprocessing\u003c/h2\u003e\n \u003cp\u003eThe data was collected using the Nihon Kohden recording system with a sampling rate of 1000 Hz. We used an in-house programmed preprocessing pipeline to denoise the intracranial EEG data before any statistical testing. The pipeline included notch filtering at 60, 120 and 180 Hz, data exclusion, and re-referencing. Excluded data includes pathological, noisy channels, and bad trials, the identification criteria of which are described as follows. We used a recently developed algorithm to identify pathological channels by the presence of pathological high-frequency oscillations (HFOs) \u003csup\u003e53\u003c/sup\u003e. Noisy channels were identified by the presence of extremely high raw amplitude (\u0026gt;\u0026thinsp;5 standard deviations [SD] across all channels) and/or the prevalence of spikes (\u0026gt;\u0026thinsp;3 median of the distribution across all channels), i.e. jumps between consecutive data points larger than 80 \u0026micro;V. Trials with extreme variations in the signal, i.e. having extreme mean (\u0026gt;\u0026thinsp;4 SD across all trials) of the absolute magnitude of the signal values, and extreme standard deviation (\u0026gt;\u0026thinsp;3.5 SD) of the variance of the signal values along the sampled timepoints, were excluded for being trial-averaged to create the time-domain evoked potential. Following data screening, we applied bi-polar re-referencing instead of common averaging, given that the stimulations were performed in a bi-polar manner, i.e., injected between a pair of electrodes. Common average reference of all non-noisy channels for the intracranial EEG recording analysis has also been experimented and compared with the bi-polar referencing scheme to make sure that the bi-polar subtraction did not remove signals. After manual inspection, we chose to use the bi-polar referencing for the further analyses.\u003c/p\u003e\n \u003cp\u003eTo extract temporal-frequency information, the continuous bipolar voltage traces were decomposed using \u003cem\u003eMorlet\u003c/em\u003e wavelet filtering (at log-spaced frequencies between 1 and 256 Hz [59 total frequencies]; each wavelet having a width of 5 cycles). The output instantaneous power timeseries has a sampling rate of 200 at each frequency. The timeseries was then epoched, baseline corrected and averaged across good trials. Values in the power spectrogram were then log-transformed to approximate normal distribution and z-scored across time and frequency. To measure how EPs are consistent across trials, we calculated the inter-trial phase coherence (ITPC) by taking the complex phase component of the wavelet decomposition, averaging across trials, and then taking the absolute. Values in the ITPC spectrogram were then square-root transformed to approximate normal distribution and then z-scored across time and frequencies. The remaining text of the paper will refer to the evoked potentials (EP), the evoked power and ITPC spectrograms as their preprocessed version by default.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec19\" class=\"Section2\"\u003e\n \u003ch2\u003eActivation detection using semi-supervised Uniform Manifold Approximation and Projection (UMAP)\u003c/h2\u003e\n \u003cdiv id=\"Sec20\" class=\"Section3\"\u003e\n \u003ch2\u003ePreliminary activation labelling\u003c/h2\u003e\n \u003cp\u003eThe key to determining a casual connectivity lies on correctly detecting a true physiological activation among the volume conductions of an electrical stimulation. We employed a semi-supervised learning approach utilizing the Uniform Manifold Approximation and Projection (UMAP) algorithm with the python toolbox \u003cem\u003eumap-learn\u003c/em\u003e (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://umap-learn.readthedocs.io/en/latest/\u003c/span\u003e\u003c/span\u003e)\u003csup\u003e54\u003c/sup\u003e. Unlike traditional dimensionality reduction techniques, semi-supervised UMAP allows for integration of labeled and unlabeled data, using our partial knowledge to guide the manifold structure of a true activation.\u003c/p\u003e\n \u003cp\u003eBefore training, we produced a preliminary activation label based on the \u0026ldquo;consensus\u0026rdquo; of the pre-defined 3 criteria: (1) \u0026ldquo;Uniform-Cutoff\u0026rdquo;: a hard threshold of EP z-score\u0026thinsp;\u0026gt;\u0026thinsp;7 at any time point in the [10, 600] ms time window. (2) \u0026ldquo;Channel-Adaptive\u0026rdquo;: stimulation-channel specific thresholds (SD\u0026thinsp;\u0026gt;\u0026thinsp;2) for peak heights and peak dominances of the EPs generated from the same bipolar stimulation. EP peaks were detected using the MATLAB \u0026ldquo;findpeaks\u0026rdquo; function for all peaks (positive or negative), and the peak and prominence distributions were constructed from all the detected peaks, upon which we applied the thresholds. (3) \u0026ldquo;Phase-locking\u0026rdquo;: individual specific thresholds with a confidence level of 95% in the distribution of connectivity index scores across all stimulating-recording pairs of the same subject. The connectivity index is based on ITPC, recently created to measure causal connectivity as discussed in a recent study\u003csup\u003e11\u003c/sup\u003e. Apart from these criteria reported from the literature, we also imposed some criteria that should be met to be physiologically alike: the voltage trace of the EP should have smooth and fast rises in the early phase (10\u0026thinsp;~\u0026thinsp;60 ms), later peaks (after 250 ms) can have lower peak heights (z-score\u0026thinsp;\u0026gt;\u0026thinsp;5), and the signal should approximate to baseline after 800 ms. When all the criteria agree, we labelled a stimulating-recording instance as either \u0026ldquo;1\u0026rdquo; (activated) or \u0026ldquo;0\u0026rdquo; (non-activated); otherwise (~\u0026thinsp;1/3 of the cases), we labelled it as \u0026ldquo;-1\u0026rdquo; (unlabeled), which is left for the algorithm to apply learned rules and decide for us.\u003c/p\u003e\n \u003c/div\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec21\" class=\"Section2\"\u003e\n \u003ch2\u003eTraining data preparation\u003c/h2\u003e\n \u003cp\u003eWe applied a semi-supervised UMAP with power and ITPC spectrograms as input data, as the spectrograms explicitly quantified the rich waveform information (components and consistency) of the EPs. We prepared the training data with the following steps. We down sampled and smoothed the spectrograms by averaging the nearest neighbors to reduce the data redundancy, save working memory for the intensive computation, and to denoise and highlight the features of interest. We used characteristic time-frequency kernels to smooth (averaging the nearest neighbors) and down-sampled the spectrograms, based on Leland McInnes k-means clustering from all stimulation-evoked potential peaks, and our prior knowledge that the first 10 ms upon a stimulation onset is mainly artifact (\u003cstrong\u003eFig. S5\u003c/strong\u003e). Therefore, we down sampled the temporal domain into 5, 25, 20, and 10 points to cover the windows of [-30, 10] ms (artifact window), [11, 117] ms, [118, 283] ms, and [284, 800] ms. The artifact window was kept in the training data for the algorithm to learn to recognize. We used more data points for representing the [11, 283] ms because most evoked potentials happen in this time window. On the frequency domain, they were down sampled into 4 datapoints to cover the [0.5, 5] Hz (delta), [5, 8] Hz (theta), [8, 15] Hz (alpha), [15, 30] Hz (beta), [30, 70] Hz (gamma) and \u0026gt;\u0026thinsp;70 Hz (high). The down sampled power and ITPC spectrograms were then vectorized and concatenated as a single vector of multiple features (number of features\u0026thinsp;=\u0026thinsp;timepoints*frequencies of power\u0026thinsp;+\u0026thinsp;timepoints*frequencies of ITPC) for each stimulating-recording instance. Each stimulating-recording vector was then concatenated to be the input training data. The input data was curated using the \u003cem\u003esklearn\u003c/em\u003e toolbox in python (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://scikit-learn.org/stable/\u003c/span\u003e\u003c/span\u003e), including missing data imputation with mean statistics (negligible cases, about 0.089% of the datapoints), and quantile transformation to distribute the data normally. When both transformed and un-transformed data were experimented, the results did not change but showed a clearer clustering effect using the transformed data; hence, the results using transformed data are reported in this paper.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec22\" class=\"Section2\"\u003e\n \u003ch2\u003eSubject-level semi-supervised UMAP learning\u003c/h2\u003e\n \u003cp\u003eThe number of latent dimensions for the UMAP mapping was set to 2, the size of the local neighborhood was set to 5, the minimum distance of the points in the latent space was set to 0.01. The rest of the parameters were kept at default. As a result, each participant\u0026rsquo;s EP data was classified into activation and non-activation clusters that were grouped at the far ends in latent space, while the few floating data points that were not clustered were considered outliers. Each person\u0026rsquo;s clustering was manually checked by (1) visualizing the spectrograms of the centroid and boundary points of each cluster, and (2) visualizing the inverse projection of the original space from the latent space. 26 out of 27 participants\u0026rsquo; data were grouped into two meaningful clusters, while one participant\u0026rsquo;s data (S21_196) with extensive amount of noisy data failed to cluster and was excluded from the following analyses.\u003c/p\u003e\n \u003cdiv id=\"Sec23\" class=\"Section3\"\u003e\n \u003ch2\u003eNeural feature extraction for activated signals\u003c/h2\u003e\n \u003cdiv id=\"Sec24\" class=\"Section4\"\u003e\n \u003ch2\u003eGroup-level UMAP learning\u003c/h2\u003e\n \u003cp\u003eTo gain an overview of the data structure of the activated spectrograms, we applied an unsupervised UMAP with input data of denoised power and ITPC spectrograms of the activated channels (number of latent dimensions\u0026thinsp;=\u0026thinsp;4). To minimize the influence of electrical artifacts, we created a non-activated spectral information (power and ITPC) template for each stimulating channel by averaging all the non-activated pairs stimulating from this channel. This non-activated spectral template contains mostly the artifact of volume conduction (\u003cstrong\u003eFig. S5\u003c/strong\u003e), as the physiological events would be averaged out when not synchronized by an activation. We used the non-activated spectral template as a control and subtracted it from all the activated spectrograms from the same stimulation channel within the same participant. When the non-activated cases for a stimulation channel are too few (\u0026lt;\u0026thinsp;10) and thus unable to generate a stable mean as a contrast, this happened only once, we excluded the case from the subsequent analyses. The remaining traces of the activated spectral information were considered denoised and fed into the group-level analysis.\u003c/p\u003e\n \u003cp\u003eWe applied a supervised UMAP with input data of denoised power and ITPC spectrograms of the activated channels, to map to the anatomical labeling of interest. The provided anatomical labels were \u0026ldquo;THAL-ipsi\u0026rdquo;, \u0026ldquo;THAL-contr\u0026rdquo;, \u0026ldquo;COR-ipsi\u0026rdquo; and \u0026ldquo;COR-contr\u0026rdquo;, corresponding to \u003cem\u003estimulating\u003c/em\u003e from thalamus/cortex and recording from the ipsilateral/contralateral hemisphere.\u003c/p\u003e\n \u003cp\u003eThe same training data preparation as mentioned before was applied. The number of latent dimensions for the UMAP mapping was set to 2, the size of the local neighborhood was set to 15, the minimum distance of the points in the latent space was set to 0.1. These parameters varied from the within-subject semi-supervised UMAP learning because the data variability within-subject is assumed to be smaller than the group data, while we aim to find commonalities that override individual differences. Therefore, the local neighborhood and the latent space are relaxed to bigger numbers (As a general principle, we experimented with different sets of parameters [within a reasonable range] to achieve the best clustering results from visual inspection). The rest of the parameters were kept at default. We further employed a hierarchical-density clustering algorithm to define the boundaries of the embedding clusters against the background noise. Finally, we investigated the centroid and averaged the presentation of the data in each cluster to verify that the clustering preserved meaningful and distinct features within the data.\u003c/p\u003e\n \u003cp\u003eTo ensure the data-driven clusters are generalizable and meaningful, we fed the resulting embedding to a series of feature evaluation procedures, including using a support vector classifier to test their generalizability to the \u0026frac14; unseen testing data. We showed that the embedding for separating the cortical and thalamus evoked responses has 75.9% predictability to the unseen data, much higher than chance rate (25%). Especially, the clusters of thalamus and cortex stimulated data are mostly well preserved, suggesting these two types of connectivity show most distinctive features as encapsulated by their power and ITPC spectrograms. We further employed a hierarchical-density clustering algorithm\u003csup\u003e55\u003c/sup\u003e to define the boundaries of the embedding clusters against the background noise, using the \u003cem\u003eHDBscan\u003c/em\u003e toolbox (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://hdbscan.readthedocs.io/en/latest/\u003c/span\u003e\u003c/span\u003e). We then investigated the averaged spectrograms (i.e. the cluster centroid) for each cluster, for a sanity check of the clustering results.\u003c/p\u003e\n \u003c/div\u003e\n \u003c/div\u003e\n \u003cdiv id=\"Sec25\" class=\"Section3\"\u003e\n \u003ch2\u003eCluster-based permutation test\u003c/h2\u003e\n \u003cp\u003eCategory-specific spectral features were determined by cluster-based permutation significance testing: (1) one-sample T-test for determining the significant clusters for each category, and (2) paired T-test for determining the significant differences between categories (with a 2 [THAL/COR]-by-2 [ipsi/contr] design). Specifically, two-sample t-tests were conducted within subjects and the individual statistics (t-scores) were fed into a group-level one-sample T-test for final inference. The testing was performed using the toolbox MNE-python (v1.5.0) (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://mne.tools/stable/index.html\u003c/span\u003e\u003c/span\u003e). The number of permutations was set to 5000, the initial cluster forming threshold was bigger than 6 t-scores, and the cluster-wise permutation confidence interval was set to 99% (i.e., \u003cem\u003eP\u003c/em\u003e\u003csub\u003ecluster\u003c/sub\u003e \u0026lt; 0.01).\u003c/p\u003e\n \u003c/div\u003e\n \u003cdiv id=\"Sec26\" class=\"Section3\"\u003e\n \u003ch2\u003eNeural feature decoding\u003c/h2\u003e\n \u003cp\u003eBased on the cluster-based permutation test, we identified three distinct neural features (i.e., significant clusters) that are specific to each anatomical category of interest (i.e., UMAP localizers). They are represented by the power and ITPC patterns in the window of ([10, 50] ms of the COR-ipsi category (averaged across all instances in the corresponding category), in the window of [70, 160] ms of the COR-contr category, and the [165, 280] in the THAL-contr category. The temporal-frequency relationships in these time ranges and anatomical categories were considered as neural features. They were decoded in each instance of stimulated/record pairs. To do this, we utilized the spectral information from the significant cluster in the group-level spectrograms as a template. We employed a sliding-window cross-correlation to match this template with the specific instances of the spectrograms, point by point in time. Different similarity (or distance) measures were experimented and turned out to produce similar results; therefore, we used the Pearson correlation \u003cem\u003er\u003c/em\u003e for its easy interpretability. The resulting timeseries of \u003cem\u003er\u003c/em\u003e has a samping rate of 200 per second, which is the same as the spectral dataset.\u003c/p\u003e\n \u003cp\u003eA sliding-window cross-correlation with the 3 identified neural features generated three temporal similarity curves for each stimulating-recoding instance, each corresponding to the neural feature\u0026rsquo;s temporal evolvement in this specific pair of stimulated/recorded sites. With the temporal similarity curve, we detected the peak similarity and the time to the peak using MATLAB \u0026ldquo;findpeak\u0026rdquo; function. Instances within each anatomical area were averaged for brain regional inferences.\u003c/p\u003e\n \u003c/div\u003e\n \u003cdiv id=\"Sec27\" class=\"Section3\"\u003e\n \u003ch2\u003eGeneral framework of statistical testing:\u003c/h2\u003e\n \u003cp\u003eThe datasets in the present study generally have a nested structure such that each condition consists of multiple electrode contacts, each with a brain area and subject identity; therefore, the mixed linear model was used as a general testing framework for model fitting, using the brain area and subject identifiers as nested random factors. Therefore, the comparisons were always made within the same stimulation/recording channels, brain areas and subjects, whenever possible. Practical model setups can vary slightly according to data structures and tested hypotheses. For significance inference, the model comparison approach with a hypothesis testing framework was adopted. Namely, a null model is compared with an alternative model which has the hypothesized effect added on top of the null model. The model with higher evidence is selected, based on a likelihood-ratio test. We additionally provided the Akaike information criterion (AIC) and Bayesian information criterion (BIC) in our supplementary table of post-hoc tests where spurious significant results may arise given the large numbers of tests, even after the false-discovery-rate (FDR) correction. The model fittings and selections were conducted using the \u003cem\u003elme4\u003c/em\u003e toolbox and the \u003cem\u003eanova\u003c/em\u003e function in R. The model specifications were provided along with the significance tables in the supplementary materials.\u003c/p\u003e\n \u003c/div\u003e\n\u003c/div\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eAcknowledgments:\u003c/strong\u003e We are grateful to the patients who agreed to participate in our research, and the staff in the EEG lab including clinical fellows and attendings whose help enabled this research. We express our gratitude to our colleagues Riki Matsumoto (Kobe - Japan), Corey Keller (Stanford - USA) and Konrad K\u0026ouml;rding (Philadelphia \u0026ndash; USA) for their valuable feedback on our findings; and to Leland McInnes (Ottawa \u0026ndash; Canada) for running a sanity check on our UMAP analysis approach.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunding:\u003c/strong\u003e This work was supported by research grants from the US National Institute of Neurological Disorders and Stroke (R01NS078396 and R21NS113024), US National Institute of Mental Health (1R01MH109954, and P50MH109429) and US National Science Foundation (BCS1358907) and US National Science Foundation (BCS1358907 and BCS1850938) to JP.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthor contributions:\u003c/strong\u003e DL contributed to the conception, acquisition, analysis, and interpretation of data and creation of new codes and measures used in the work and writing of the manuscript; JRS contributed to the analysis and interpretation of data; ZL contributed to the acquisition of data; JP contributed to conception and design of the work, acquisition and interpretation of data and writing of the manuscript. All authors provided significant feedback throughout the study and manuscript preparation and have approved the submitted version.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompeting interests:\u003c/strong\u003e The authors declare no competing interests.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eSporns, O., Chialvo, D. R., Kaiser, M. \u0026amp; Hilgetag, C. C. Organization, development and function of complex brain networks. \u003cem\u003eTrends Cogn. Sci.\u003c/em\u003e \u003cstrong\u003e8\u003c/strong\u003e, 418-425 (2004). \u003c/li\u003e\n\u003cli\u003eShine, J. M., Lewis, L. D., Garrett, D. D. \u0026amp; Hwang, K. The impact of the human thalamus on brain-wide information processing. \u003cem\u003eNature reviews. Neuroscience\u003c/em\u003e (2023). https://doi.org:10.1038/s41583-023-00701-0\u003c/li\u003e\n\u003cli\u003eMatsumoto, R.\u003cem\u003e et al.\u003c/em\u003e Functional connectivity in the human language system: a cortico-cortical evoked potential study. \u003cem\u003eBrain\u003c/em\u003e \u003cstrong\u003e127\u003c/strong\u003e, 2316-2330 (2004). \u003c/li\u003e\n\u003cli\u003eMatsumoto, R.\u003cem\u003e et al.\u003c/em\u003e Functional connectivity in human cortical motor system: a cortico-cortical evoked potential study. \u003cem\u003eBrain\u003c/em\u003e \u003cstrong\u003e130\u003c/strong\u003e, 181-197 (2007). \u003c/li\u003e\n\u003cli\u003eWu, T. Q.\u003cem\u003e et al.\u003c/em\u003e Multisite thalamic recordings to characterize seizure propagation in the human brain. \u003cem\u003eBrain\u003c/em\u003e \u003cstrong\u003e146\u003c/strong\u003e, 2792-2802 (2023). https://doi.org:10.1093/brain/awad121\u003c/li\u003e\n\u003cli\u003eZumsteg, D., Lozano, A. M. \u0026amp; Wennberg, R. A. Depth electrode recorded cerebral responses with deep brain stimulation of the anterior thalamus for epilepsy. \u003cem\u003eClin Neurophysiol\u003c/em\u003e \u003cstrong\u003e117\u003c/strong\u003e, 1602-1609 (2006). https://doi.org:10.1016/j.clinph.2006.04.008\u003c/li\u003e\n\u003cli\u003eOjeda Valencia, G.\u003cem\u003e et al.\u003c/em\u003e Signatures of Electrical Stimulation Driven Network Interactions in the Human Limbic System. \u003cem\u003eThe Journal of Neuroscience: The Official Journal of the Society for Neuroscience\u003c/em\u003e \u003cstrong\u003e43\u003c/strong\u003e, 6697-6711 (2023). https://doi.org:10.1523/JNEUROSCI.2201-22.2023\u003c/li\u003e\n\u003cli\u003eMiller, K. J., Muller, K. R. \u0026amp; Hermes, D. Basis profile curve identification to understand electrical stimulation effects in human brain networks. \u003cem\u003ePLoS Comput Biol\u003c/em\u003e \u003cstrong\u003e17\u003c/strong\u003e, e1008710 (2021). https://doi.org:10.1371/journal.pcbi.1008710\u003c/li\u003e\n\u003cli\u003eJedynak, M.\u003cem\u003e et al.\u003c/em\u003e Variability of Single Pulse Electrical Stimulation Responses Recorded with Intracranial Electroencephalography in Epileptic Patients. \u003cem\u003eBrain Topogr\u003c/em\u003e \u003cstrong\u003e36\u003c/strong\u003e, 119-127 (2023). https://doi.org:10.1007/s10548-022-00928-7\u003c/li\u003e\n\u003cli\u003eBecht, E.\u003cem\u003e et al.\u003c/em\u003e Dimensionality reduction for visualizing single-cell data using UMAP. \u003cem\u003eNat Biotechnol\u003c/em\u003e (2018). https://doi.org:10.1038/nbt.4314\u003c/li\u003e\n\u003cli\u003eStieger, J. R.\u003cem\u003e et al.\u003c/em\u003e Cross regional coordination of neural activity in the human brain during autobiographical self-referential processing. \u003cem\u003ePNAS\u003c/em\u003e \u003cstrong\u003eIn Press\u003c/strong\u003e (2023). \u003c/li\u003e\n\u003cli\u003eKunieda, T., Yamao, Y., Kikuchi, T. \u0026amp; Matsumoto, R. New Approach for Exploring Cerebral Functional Connectivity: Review of Cortico-cortical Evoked Potential. \u003cem\u003eNeurologia medico-chirurgica\u003c/em\u003e \u003cstrong\u003e55\u003c/strong\u003e, 374-382 (2015). https://doi.org:10.2176/nmc.ra.2014-0388\u003c/li\u003e\n\u003cli\u003eKeller, C. J.\u003cem\u003e et al.\u003c/em\u003e Intrinsic functional architecture predicts electrically evoked responses in the human brain. \u003cem\u003eProceedings of the National Academy of Sciences\u003c/em\u003e \u003cstrong\u003e108\u003c/strong\u003e, 10308-10313 (2011). https://doi.org:10.1073/pnas.1019750108\u003c/li\u003e\n\u003cli\u003eKeller, C. J.\u003cem\u003e et al.\u003c/em\u003e Mapping human brain networks with cortico-cortical evoked potentials. \u003cem\u003ePhilos Trans R Soc Lond B Biol Sci\u003c/em\u003e \u003cstrong\u003e369\u003c/strong\u003e, 20130528 (2014). https://doi.org:10.1098/rstb.2013.0528\u003c/li\u003e\n\u003cli\u003eTogo, M.\u003cem\u003e et al.\u003c/em\u003e Distinct connectivity patterns in human medial parietal cortices: Evidence from standardized connectivity map using cortico-cortical evoked potential. \u003cem\u003eNeuroImage\u003c/em\u003e \u003cstrong\u003e263\u003c/strong\u003e, 119639 (2022). https://doi.org:10.1016/j.neuroimage.2022.119639\u003c/li\u003e\n\u003cli\u003eFogerson, P. M. \u0026amp; Huguenard, J. R. Tapping the Brakes: Cellular and Synaptic Mechanisms that Regulate Thalamic Oscillations. \u003cem\u003eNeuron\u003c/em\u003e \u003cstrong\u003e92\u003c/strong\u003e, 687-704 (2016). https://doi.org:10.1016/j.neuron.2016.10.024\u003c/li\u003e\n\u003cli\u003eBuzsaki, G. \u003cem\u003eRhythms of the Brain\u003c/em\u003e. (Oxford University Press, 2006).\u003c/li\u003e\n\u003cli\u003eGroppe, D. M.\u003cem\u003e et al.\u003c/em\u003e Dominant frequencies of resting human brain activity as measured by the electrocorticogram. \u003cem\u003eNeuroimage\u003c/em\u003e \u003cstrong\u003e79\u003c/strong\u003e, 223-233 (2013). https://doi.org:10.1016/j.neuroimage.2013.04.044\u003c/li\u003e\n\u003cli\u003eHacker, C. D., Snyder, A. Z., Pahwa, M., Corbetta, M. \u0026amp; Leuthardt, E. C. Frequency-specific electrophysiologic correlates of resting state fMRI networks. \u003cem\u003eNeuroimage\u003c/em\u003e \u003cstrong\u003e149\u003c/strong\u003e, 446-457 (2017). https://doi.org:10.1016/j.neuroimage.2017.01.054\u003c/li\u003e\n\u003cli\u003eFisher, R.\u003cem\u003e et al.\u003c/em\u003e Electrical stimulation of the anterior nucleus of thalamus for treatment of refractory epilepsy. \u003cem\u003eEpilepsia\u003c/em\u003e \u003cstrong\u003e51\u003c/strong\u003e, 899-908 (2010). https://doi.org:10.1111/j.1528-1167.2010.02536.x\u003c/li\u003e\n\u003cli\u003eEngel, J., Jr. \u0026amp; Pitkanen, A. Biomarkers for epileptogenesis and its treatment. \u003cem\u003eNeuropharmacology\u003c/em\u003e \u003cstrong\u003e167\u003c/strong\u003e, 107735 (2020). https://doi.org:10.1016/j.neuropharm.2019.107735\u003c/li\u003e\n\u003cli\u003eCaciagli, L., Bernhardt, B. C., Hong, S.-J., Bernasconi, A. \u0026amp; Bernasconi, N. Functional network alterations and their structural substrate in drug-resistant epilepsy. \u003cem\u003eFront Neurosci\u003c/em\u003e \u003cstrong\u003e8\u003c/strong\u003e, 411 (2014). https://doi.org:10.3389/fnins.2014.00411\u003c/li\u003e\n\u003cli\u003eFleury, M.\u003cem\u003e et al.\u003c/em\u003e Episodic memory network connectivity in temporal lobe epilepsy. \u003cem\u003eEpilepsia\u003c/em\u003e \u003cstrong\u003e63\u003c/strong\u003e, 2597-2622 (2022). https://doi.org:10.1111/epi.17370\u003c/li\u003e\n\u003cli\u003ePittau, F., Grova, C., Moeller, F., Dubeau, F. \u0026amp; Gotman, J. Patterns of altered functional connectivity in mesial temporal lobe epilepsy. \u003cem\u003eEpilepsia\u003c/em\u003e \u003cstrong\u003e53\u003c/strong\u003e, 1013-1023 (2012). https://doi.org:10.1111/j.1528-1167.2012.03464.x\u003c/li\u003e\n\u003cli\u003eRoger, E.\u003cem\u003e et al.\u003c/em\u003e Hubs disruption in mesial temporal lobe epilepsy. A resting-state fMRI study on a language-and-memory network. \u003cem\u003eHuman Brain Mapping\u003c/em\u003e \u003cstrong\u003e41\u003c/strong\u003e, 779-796 (2020). https://doi.org:10.1002/hbm.24839\u003c/li\u003e\n\u003cli\u003eLi, L.\u003cem\u003e et al.\u003c/em\u003e Topographical reorganization of brain functional connectivity during an early period of epileptogenesis. \u003cem\u003eEpilepsia\u003c/em\u003e \u003cstrong\u003e62\u003c/strong\u003e, 1231-1243 (2021). https://doi.org:10.1111/epi.16863\u003c/li\u003e\n\u003cli\u003eMazrooyisebdani, M.\u003cem\u003e et al.\u003c/em\u003e Graph Theory Analysis of Functional Connectivity Combined with Machine Learning Approaches Demonstrates Widespread Network Differences and Predicts Clinical Variables in Temporal Lobe Epilepsy. \u003cem\u003eBrain Connect\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, 39-50 (2020). https://doi.org:10.1089/brain.2019.0702\u003c/li\u003e\n\u003cli\u003eOfer, I.\u003cem\u003e et al.\u003c/em\u003e Association between seizure freedom and default mode network reorganization in patients with unilateral temporal lobe epilepsy. \u003cem\u003eEpilepsy Behav\u003c/em\u003e \u003cstrong\u003e90\u003c/strong\u003e, 238-246 (2019). https://doi.org:10.1016/j.yebeh.2018.10.025\u003c/li\u003e\n\u003cli\u003eLiao, W.\u003cem\u003e et al.\u003c/em\u003e Default mode network abnormalities in mesial temporal lobe epilepsy: A study combining fMRI and DTI. \u003cem\u003eHuman Brain Mapping\u003c/em\u003e \u003cstrong\u003e32\u003c/strong\u003e, 883-895 (2011). https://doi.org:10.1002/hbm.21076\u003c/li\u003e\n\u003cli\u003eWidjaja, E., Zamyadi, M., Raybaud, C., Snead, O. C. \u0026amp; Smith, M. L. Abnormal Functional Network Connectivity among Resting-State Networks in Children with Frontal Lobe Epilepsy. \u003cem\u003eAJNR Am J Neuroradiol\u003c/em\u003e \u003cstrong\u003e34\u003c/strong\u003e, 2386-2392 (2013). https://doi.org:10.3174/ajnr.A3608\u003c/li\u003e\n\u003cli\u003eZhang, Z.\u003cem\u003e et al.\u003c/em\u003e Impaired attention network in temporal lobe epilepsy: A resting FMRI study. \u003cem\u003eNeuroscience Letters\u003c/em\u003e \u003cstrong\u003e458\u003c/strong\u003e, 97-101 (2009). https://doi.org:10.1016/j.neulet.2009.04.040\u003c/li\u003e\n\u003cli\u003eKeller, C. J.\u003cem\u003e et al.\u003c/em\u003e Mapping human brain networks with cortico-cortical evoked potentials. \u003cem\u003ePhilos Trans R Soc Lond B Biol Sci\u003c/em\u003e \u003cstrong\u003e369\u003c/strong\u003e (2014). https://doi.org:10.1098/rstb.2013.0528\u003c/li\u003e\n\u003cli\u003eVeit, M. J.\u003cem\u003e et al.\u003c/em\u003e Temporal order of signal propagation within and across intrinsic brain networks. \u003cem\u003eProc Natl Acad Sci U S A\u003c/em\u003e \u003cstrong\u003e118\u003c/strong\u003e (2021). https://doi.org:10.1073/pnas.2105031118\u003c/li\u003e\n\u003cli\u003eSeguin, C.\u003cem\u003e et al.\u003c/em\u003e Communication dynamics in the human connectome shape the cortex-wide propagation of direct electrical stimulation. \u003cem\u003eNeuron\u003c/em\u003e \u003cstrong\u003e0\u003c/strong\u003e (2023). https://doi.org:10.1016/j.neuron.2023.01.027\u003c/li\u003e\n\u003cli\u003eGuo, Z. H.\u003cem\u003e et al.\u003c/em\u003e Epileptogenic network of focal epilepsies mapped with cortico-cortical evoked potentials. \u003cem\u003eClin Neurophysiol\u003c/em\u003e \u003cstrong\u003e131\u003c/strong\u003e, 2657-2666 (2020). https://doi.org:10.1016/j.clinph.2020.08.012\u003c/li\u003e\n\u003cli\u003eArthuis, M.\u003cem\u003e et al.\u003c/em\u003e Impaired consciousness during temporal lobe seizures is related to increased long-distance cortical-subcortical synchronization. \u003cem\u003eBrain\u003c/em\u003e \u003cstrong\u003e132\u003c/strong\u003e, 2091-2101 (2009). https://doi.org:10.1093/brain/awp086\u003c/li\u003e\n\u003cli\u003eGuye, M.\u003cem\u003e et al.\u003c/em\u003e The role of corticothalamic coupling in human temporal lobe epilepsy. \u003cem\u003eBrain\u003c/em\u003e \u003cstrong\u003e129\u003c/strong\u003e, 1917-1928 (2006). https://doi.org:10.1093/brain/awl151\u003c/li\u003e\n\u003cli\u003ePizzo, F.\u003cem\u003e et al.\u003c/em\u003e The Ictal Signature of Thalamus and Basal Ganglia in Focal Epilepsy: A SEEG Study. \u003cem\u003eNeurology\u003c/em\u003e \u003cstrong\u003e96\u003c/strong\u003e, e280-e293 (2021). https://doi.org:10.1212/WNL.0000000000011003\u003c/li\u003e\n\u003cli\u003eFilipescu, C.\u003cem\u003e et al.\u003c/em\u003e The effect of medial pulvinar stimulation on temporal lobe seizures. \u003cem\u003eEpilepsia\u003c/em\u003e \u003cstrong\u003e60\u003c/strong\u003e, e25-e30 (2019). https://doi.org:10.1111/epi.14677\u003c/li\u003e\n\u003cli\u003eEvangelista, E.\u003cem\u003e et al.\u003c/em\u003e Does the Thalamo-Cortical Synchrony Play a Role in Seizure Termination? \u003cem\u003eFront Neurol\u003c/em\u003e \u003cstrong\u003e6\u003c/strong\u003e, 192 (2015). https://doi.org:10.3389/fneur.2015.00192\u003c/li\u003e\n\u003cli\u003eGadot, R., Korst, G., Shofty, B., Gavvala, J. R. \u0026amp; Sheth, S. A. Thalamic stereoelectroencephalography in epilepsy surgery: a scoping literature review. \u003cem\u003eJ Neurosurg\u003c/em\u003e, 1-16 (2022). https://doi.org:10.3171/2022.1.JNS212613\u003c/li\u003e\n\u003cli\u003eIlyas, A., Tandon, N. \u0026amp; Lhatoo, S. D. Thalamic neuromodulation for epilepsy: A clinical perspective. \u003cem\u003eEpilepsy Res\u003c/em\u003e \u003cstrong\u003e183\u003c/strong\u003e, 106942 (2022). https://doi.org:10.1016/j.eplepsyres.2022.106942\u003c/li\u003e\n\u003cli\u003eChaitanya, G.\u003cem\u003e et al.\u003c/em\u003e Robot-assisted stereoelectroencephalography exploration of the limbic thalamus in human focal epilepsy: implantation technique and complications in the first 24 patients. \u003cem\u003eNeurosurg Focus\u003c/em\u003e \u003cstrong\u003e48\u003c/strong\u003e, E2 (2020). https://doi.org:10.3171/2020.1.FOCUS19887\u003c/li\u003e\n\u003cli\u003eRomeo, A.\u003cem\u003e et al.\u003c/em\u003e Early ictal recruitment of midline thalamus in mesial temporal lobe epilepsy. \u003cem\u003eAnn Clin Transl Neurol\u003c/em\u003e \u003cstrong\u003e6\u003c/strong\u003e, 1552-1558 (2019). https://doi.org:10.1002/acn3.50835\u003c/li\u003e\n\u003cli\u003eMcKhann, G. M. Editorial. Dulling the double-edged sword of human SEEG research. \u003cem\u003eNeurosurg Focus\u003c/em\u003e \u003cstrong\u003e48\u003c/strong\u003e, E3 (2020). https://doi.org:10.3171/2020.1.FOCUS2069\u003c/li\u003e\n\u003cli\u003eGroppe, D. M.\u003cem\u003e et al.\u003c/em\u003e iELVis: An open source MATLAB toolbox for localizing and visualizing human intracranial electrode data. \u003cem\u003eJ Neurosci Methods\u003c/em\u003e \u003cstrong\u003e281\u003c/strong\u003e, 40-48 (2017). https://doi.org:10.1016/j.jneumeth.2017.01.022\u003c/li\u003e\n\u003cli\u003eFischl, B. FreeSurfer. \u003cem\u003eNeuroImage\u003c/em\u003e \u003cstrong\u003e62\u003c/strong\u003e, 774-781 (2012). \u003c/li\u003e\n\u003cli\u003eJenkinson, M., Beckmann, C. F., Behrens, T. E., Woolrich, M. W. \u0026amp; Smith, S. M. Fsl. \u003cem\u003eNeuroimage\u003c/em\u003e \u003cstrong\u003e62\u003c/strong\u003e, 782-790 (2012). https://doi.org:10.1016/j.neuroimage.2011.09.015\u003c/li\u003e\n\u003cli\u003eJenkinson, M. \u0026amp; Smith, S. A global optimisation method for robust affine registration of brain images. \u003cem\u003eMedical image analysis\u003c/em\u003e \u003cstrong\u003e5\u003c/strong\u003e, 143-156 (2001). \u003c/li\u003e\n\u003cli\u003eGreve, D. N. \u0026amp; Fischl, B. Accurate and robust brain image alignment using boundary-based registration. \u003cem\u003eNeuroimage\u003c/em\u003e \u003cstrong\u003e48\u003c/strong\u003e, 63-72 (2009). https://doi.org:10.1016/j.neuroimage.2009.06.060\u003c/li\u003e\n\u003cli\u003ePapademetris, X.\u003cem\u003e et al.\u003c/em\u003e BioImage Suite: An integrated medical image analysis suite: An update. \u003cem\u003eThe insight journal\u003c/em\u003e \u003cstrong\u003e2006\u003c/strong\u003e, 209 (2006). \u003c/li\u003e\n\u003cli\u003eSu, J. H.\u003cem\u003e et al.\u003c/em\u003e Thalamus Optimized Multi Atlas Segmentation (THOMAS): fast, fully automated segmentation of thalamic nuclei from structural MRI. \u003cem\u003eNeuroimage\u003c/em\u003e \u003cstrong\u003e194\u003c/strong\u003e, 272-282 (2019). https://doi.org:10.1016/j.neuroimage.2019.03.021\u003c/li\u003e\n\u003cli\u003eLiu, S. \u0026amp; Parvizi, J. Cognitive refractory state caused by spontaneous epileptic high-frequency oscillations in the human brain. \u003cem\u003eSci Transl Med\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e (2019). https://doi.org:10.1126/scitranslmed.aax7830\u003c/li\u003e\n\u003cli\u003eMcInnes, L., Healy, J. \u0026amp; Melville, J. (arXiv, 2020).\u003c/li\u003e\n\u003cli\u003eCampello, R. J. G. B., Moulavi, D. \u0026amp; Sander, J. (eds Jian Pei\u003cem\u003e et al.\u003c/em\u003e) 160-172 (Springer).\u003c/li\u003e\n\u003cli\u003eBates, D., M\u0026auml;chler, M., Bolker, B. \u0026amp; Walker, S. Fitting Linear Mixed-Effects Models Using lme4. \u003cem\u003eJournal of Statistical Software\u003c/em\u003e \u003cstrong\u003e67\u003c/strong\u003e, 1-48 (2015). https://doi.org:10.18637/jss.v067.i01\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"nature-portfolio","isNatureJournal":true,"hasQc":false,"allowDirectSubmit":false,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"","title":"Nature Portfolio","twitterHandle":"","acdcEnabled":false,"dfaEnabled":false,"editorialSystem":"ejp","reportingPortfolio":"","inReviewEnabled":true,"inReviewRevisionsEnabled":false},"keywords":"","lastPublishedDoi":"10.21203/rs.3.rs-4366486/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-4366486/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eThe brain's functional architecture is intricately shaped by causal connections between its cortical and subcortical structures. Here, we studied 27 participants with 4864 electrodes implanted across the anterior, mediodorsal, and pulvinar thalamic regions, and the cortex. Using data from electrical stimulation procedures and a data-driven approach informed by neurophysiological standards, we dissociated three unique spectral patterns generated by the perturbation of a given brain area. Among these, a novel waveform emerged, marked by delayed-onset slow oscillations in both ipsilateral and contralateral cortices following thalamic stimulations, suggesting a mechanism by which a thalamic site can influence bilateral cortical activity. Moreover, cortical stimulations evoked earlier signals in the thalamus than in other connected cortical areas suggesting that the thalamus receives a copy of signals before they are exchanged across the cortex. Our causal connectivity data can be used to inform biologically-inspired computational models of the functional architecture of the brain.\u003c/p\u003e","manuscriptTitle":"Causal Cortical and Thalamic Connections in the Human Brain","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2024-05-30 01:12:28","doi":"10.21203/rs.3.rs-4366486/v1","editorialEvents":[],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"nature-neuroscience","isNatureJournal":true,"hasQc":false,"allowDirectSubmit":false,"externalIdentity":"neuro","sideBox":"Learn more about [Nature Neuroscience](http://www.nature.com/neuro/)","snPcode":"","submissionUrl":"","title":"Nature Neuroscience","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"ejp","reportingPortfolio":"Nature Research","inReviewEnabled":true,"inReviewRevisionsEnabled":false}}],"origin":"","ownerIdentity":"10572f1a-c73a-4008-953a-885cbf4c617b","owner":[],"postedDate":"May 30th, 2024","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"published-in-journal","subjectAreas":[{"id":32353680,"name":"Biological sciences/Neuroscience/Computational neuroscience/Dynamical systems"},{"id":32353681,"name":"Biological sciences/Neuroscience/Diseases of the nervous system/Epilepsy"},{"id":32353682,"name":"Biological sciences/Neuroscience/Cognitive neuroscience/Consciousness"}],"tags":[],"updatedAt":"2025-07-16T07:09:25+00:00","versionOfRecord":{"articleIdentity":"rs-4366486","link":"https://doi.org/10.1038/s41593-025-02009-x","journal":{"identity":"nature-neuroscience","isVorOnly":false,"title":"Nature Neuroscience"},"publishedOn":"2025-07-15 04:00:00","publishedOnDateReadable":"July 15th, 2025"},"versionCreatedAt":"2024-05-30 01:12:28","video":"","vorDoi":"10.1038/s41593-025-02009-x","vorDoiUrl":"https://doi.org/10.1038/s41593-025-02009-x","workflowStages":[]},"version":"v1","identity":"rs-4366486","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-4366486","identity":"rs-4366486","version":["v1"]},"buildId":"-HB7Z8yhvgn0wM9Nzuekk","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}