Single unit-derived connectivity networks in tuberous sclerosis complex reveal propensity for network hypersynchrony driven by tuber-tuber interactions

This prospective study was approved by the London-Brent Research Ethics Committee on behalf of the UK Health Research Authority (IRAS ID 255823). In line with UK law, the parents of all included patients provided informed consent for study participation. All research was conducted in accordance with relevant guidelines/regulations and the Declaration of Helsinki.

Participants

Three patients with TSC were included in the study (Table 1). All were selected to undergo stereoelectroencephalography (SEEG) as part of their clinical work-up for drug resistant epilepsy. Clinical decision making, including decision to undergo SEEG and the number and location of electrodes were not affected by involvement in the study. All underwent predominantly frontal explorations with some electrodes extending into the insula and parietal cortices.

Clinical workflows for robot-assisted SEEG implantation were identical to previously published reports and post-operative CT scans were registered to the pre-operative CT scan using in-house scripts¹⁴. Electrode localisation was performed by members of the clinical team using the SEEG assistant tool in 3D slicer¹⁵.

Microelectrode recording and data processing

Enrolled patients underwent insertion of three (NNAIS05), two (NNAIS06) and 5 (NNAIS07) hybrid micro-macro electrodes (AdTech MM16A-SP05X-000), each with 10 microelectrode contacts along the shaft (Fig. 2a). Hybrid electrodes replaced existing planned macroelectrode trajectories and were chosen based on the length of the planned electrode trajectories; the hybrid micro-macro electrodes had a fixed active length of 28 mm and it was therefore not possible to specifically target the purported epileptogenic tubers. Post hoc, it was revealed that none of the isolated neurons were from putative epileptogenic tubers that were subsequently resected or ablated. Macro channels were connected to the clinical system (Natus Quantum) to record at 2 kHz whilst micro channels were connected to a research system with pre-amplification (Ripple Summit) to record at 30 kHz. The same white matter macro contacts were chosen as ground and reference for both systems and connected to each other to reduce noise.

Microelectrode recordings were conducted for up to 30 min per day for each recording day possible. This was in the awake state during the day whilst the child was sitting in bed, either watching the television or playing puzzle/activity games. Recordings were performed before the first recorded seizure in all patients.

Microelectrode recordings (Fig. 2b) were processed using Matlab R2020b, FieldTrip & WaveClus3^16,17. On channels with acceptable impedances on impedance testing (50–400 kOhms) and local field potential (LFP) visible at low frequencies, a discrete fourier transform filter (at 50,100 and 150 Hz) was applied to reduce line noise. Spikes were extracted using positive and negative detection at a 4 standard deviation above the noise threshold and subsequently sorted using WaveClus3 and manually selected based on a combination of the temperature plots and firing rate. Putative units were excluded if > 2% of the inter spike intervals (ISIs) were < 2ms and their firing rates were less than 0.5 Hz. As we wanted to sample all interictal connectivity, interictal epileptiform discharges (IEDs) were not sought out and continuous data was used¹⁸.

Each spike train was assigned a location based on the fusion of the post-operative CT to MRI as either within the tuber, in the transmantle tail or in radiologically unaffected cortex (Fig. 2c). This was possible as the predominant SEEG implantation strategy at our institute is to sample putative tubers, the depths reaching the white matter tails, with occasional electrodes in normal areas of cortex guided by the non-invasive presurgical evaluation data. Spike properties were assessed using the frequency, the spike half width and coefficient of variance of the ISIs (standard deviation of ISIs / mean ISI)¹⁹.

Generalised linear models were used to generate post-spike filters (PSFs) and coupling filters (CFs) for each neuron within each recording session (Fig. 1). These filters are statistically robust models of auto- and cross-correlation properties of the spike trains, allowing an assessment of the statistical influence of past firing (PSF) and the past firing of other neurons (CF) on future firing^11,12,20. The GLM has a baseline firing rate parameter controlling the average propensity of a neuron to fire and all other effects are encoded as a time-varying gain functions called filters that multiplicatively modulate the firing rate. The PSF encodes the modulation of a neuron’s firing rate given it has just fired, whilst the CF encodes the modulation of a neuron’s firing rate given that another neuron has just fired (Fig. 1). The CF therefore acts as a directed connectivity measure between the influencing neuron and target neuron. Following normalisation of the filters to enable shape comparisons, two main analyses were performed as follows.

First, to quantify variation in auto- or cross-correlation structures defined by the PSF and CFs, we decomposed the filter time series using principal component analysis (PCA). Differences in the first 2 principal components were compared between histological regions (3 groups for PSFs and 9 groups for CFs).

Second, we developed novel, network-based measures to assess autonomous synchrogenicity and network synchrogenicity (Fig. 1). Because healthy brains must constantly switch between synchronous and asynchronous states, we hypothesised that the mechanisms of hypersynchrony are due to abnormal responses to synchronous inputs. When a region of brain tissue receives inputs that drive neurons to fire at the same time, the dynamical response to those inputs will produce outputs that can either be synchronous or asynchronous. During interictal, asynchronous activity, we cannot directly observe the drive toward synchrony, but we can estimate this response using the fitted GLMs from above.

Autonomous synchrogenicity

Given two neurons with normalised PSFs, and, we define the dot product of their PSFs as.

where T is duration of the PSFs. \(\:{W}_{ij}\) measures the similarity of the predicted future firing of neurons i and j given that they fired simultaneously. For n neurons, this n x n matrix of pairwise dot products \(\:W = \left[ {{W_{ij}}} \right]\)can be viewed as a weighted adjacency matrix for an abstract network among neurons from which we computed the following three common network measures using the Brain Connectivity Toolbox²¹:

• Strength (i.e., weighted degree): Node strength is the sum of weights of edges connected to the node. High strength implies a densely connected node.

• Eigenvector Centrality: Eigenvector centrality is a self-referential measure of centrality; nodes have high eigenvector centrality if they connect to other nodes that have high eigenvector centrality. A high eigenvector centrality implies a high level of influence on the network.

• Clustering Coefficient: The clustering coefficient is the weight of triangles around a node and is equivalent to the weight of node’s neighbours that are neighbours of each other. A high clustering coefficient implies that a node’s neighbours are all densely connected to each other, while a low clustering coefficient indicates that a node’s neighbours are spread out across the network.

Each abstract network measure separately encapsulates the idea that some nodes are more tightly interconnected within the abstract network. Functionally, neurons with high values for strength and eigenvector centrality and low values for clustering coefficient are predicted by the GLM to produce time-correlated outputs with many other neurons in the ensemble, provided their inputs cause them to spike at the same time. This could happen, for example, if a subset of neurons had highly stereotyped oscillatory PSFs. Thus, these network measures capture our notion of autonomous synchrogenicity and comparing between histological regions facilitates an assessment of which regions drive the synchrogenicity of the network.

Network synchrogenicity

Given two pairs of neurons with normalised CFs, and, we define the dot product of their CFs as.

where T is duration of the CFs. \(\:{W}_{(i,j)(k,l)}\) measures the similarity of predicted future firing of neurons j and k given that neurons i and l fired simultaneously. The complete matrix of pairwise dot products \(\:W = \left[ {{W_{(i,j)(k,l)}}} \right]\)can be viewed as a weighted adjacency matrix for an abstract network. Note that the neurons i, j, k, and l need not be distinct.

For n neurons, this n² x n² matrix of pairwise dot products \(\:W = \left[ {{W_{(i,j)(k,l)}}} \right]\)can be viewed as a weighted adjacency matrix for an abstract network among pairs of neurons, from which we can compute strength, clustering coefficient, and eigenvector centrality, as above. Analogous to autonomous synchrogencity above, neuron pairs with high values for strength and eigenvector centrality and low values for clustering coefficient are predicted by the GLM to produce time-correlated outputs with many other neurons in the ensemble, provided their inputs cause them to spike at the same time, this time based on their cross-correlation patterns. The important distinction in this case is that, unlike autonomous synchrogenicity, the inputs and outputs can be different neurons, so synchronous inputs in one part of the network can propagate to generate synchronous outputs at another part of the network. Thus, these network measures capture our notion of network synchrogenicity. Again, comparing between histological regions facilitates an assessment of which regions drive the synchrogenicity of the network.

To compare network measures between the different histological regions, we used standard linear models. To adjust for networks of different sizes, we compared the ranks of the network measures within each subject and each recording session, normalised to a scale between 0 and 1. Statistical analyses were performed using parametric tests. P-values < 0.05 were considered statistically significant.

Code and data availability

Code used for this study is available at https://github.com/aswinchari/SingleUnitConnectivity/ and data (the single unit spike trains after spike sorting, subsequent filters and networks) are available at https://osf.io/ab39m/.