Computational modeling of seizure onset patterns to underpin their underlying mechanisms ======================================================================================== * Leila Abrishami Shokooh * Frédéric Lesage * Dang Khoa Nguyen ## Abstract In the study of epilepsy, it is of crucial importance to understand the transition from interictal into ictal activities (ictogenesis). Different mechanisms have been suggested for the generation of ictal activity; yet, it remains unclear whether different physiological mechanisms underly different seizure onset patterns. Herein, by implementing a computational model that takes into account some of the most relevant physiological events (e.g., depolarization block, collapse, and recovery of inhibitory activities) and different scenarios of imbalanced excitatory-inhibitory activities, we explored if seizures with different onset patterns stem from different underlying mechanisms. Our model revealed that depending on the excitation level, seizures could be generated due to both enhancement and collapse of inhibition for specific range of parameters. Successfully reproducing some of the commonly observed seizure onset patterns, our findings indicated that different onset patterns can arise from different underlying mechanisms. **Significance Statement** Various seizure onset patterns have been reported; however, it yet remains unknown whether seizures with distinct onset patterns originate from different underlying mechanisms. The common belief on seizure generation focuses on the imbalance between synaptic excitation and inhibition which has led to the identification of distinct and, in some cases, even contradictory mechanisms for seizure initiation. In this study, by incorporating some of these various physiological mechanisms in a unified framework, we reproduced some commonly observed seizure onset patterns. Our results suggest the existence of different mechanisms responsible for the generation of seizures with distinct onset patterns which can enhance our understanding of seizure generation mechanisms with significant implications on developing therapeutic measures in seizure control. ## Introduction Epilepsy, second only to stroke, is one of the most common chronic neurological disorders. It affects 50 million people worldwide [1] and is characterized by recurrent spontaneous seizures, i.e., sudden transient behavioral disturbances due to abnormal excessive neural discharges [2]. About 30% of patients have drug-resistant epilepsy [3–5] which could be due to the diverse network and cellular mechanisms responsible for seizure generation [6]. Variable electrophysiological features have been reported for seizures in intracranial EEG studies which has resulted in the definition of specific patterns for seizure onset [7,8]. Some onset patterns have been associated with the activity of specific brain regions [9–11] and synaptic mechanisms [12,13]. For instance, using an in vitro optogenetic experimental setting, it has been reported that seizures with Low Voltage Fast Activity (LVFA) onset are generated due to the involvement of hyperpolarizing inhibitory events while the generation of seizures with hypersynchronous onset depends on enhanced excitation [12,13]. These findings raise the question as to whether distinct onset patterns of ictal activity arise from fundamentally different mechanisms. The mechanisms underlying seizure generation involve the interaction of many complex cellular and synaptic components [14,15] which could govern the emergence of distinct iEEG ictal patterns [16–18]. Studying the imbalance between synaptic excitation and inhibition has been the main focus of the research on ictogenesis [13,19–24]. It was commonly believed that abnormally enhanced glutamatergic excitation and collapse of GABAergic inhibition were pro-epileptic and the main driver for seizure generation, whereas the increased GABAergic activity was viewed to protect against seizure generation [25–28]. However, this notion might require reconsideration in light of more recent clinical and experimental data that have offered evidence for the critical role of the GABAergic network in focal seizures. More specifically, some studies have reported enhancement of inhibitory activity before and at the very onset of ictal activity [29–33], casting doubt on the suggested role for the collapse of inhibition in seizure generation. In summary, distinct and even contradictory mechanisms have been suggested to underly seizure initiation. Incorporating these distinct mechanisms in a unified framework would enable us to compare the impact of these different mechanisms in seizure generation and would aid in providing an insight into whether different underlying mechanisms could lead to seizures with different onset patterns. To investigate the underlying mechanisms for generation of some of the most common patterns of ictal activity, we constructed a neural mass model and, by updating, expanding, and exploring a well-established framework [34,35], we generated different patterns of seizure generation. More specifically, we considered different scenarios of disturbed balance between excitatory and inhibitory activities and identified a range of parameters producing distinct onset patterns providing interpretation in terms of underlying physiological mechanisms. The model successfully reproduced a number of commonly observed onset patterns including Low Voltage Fast Activity (LVFA), spindles of alpha/beta, rhythmic Spike and Wave (rSW/rPSW), burst suppression, High Amplitude Fast Activity (HAFA) and rhythmic spikes/sharp waves. The results, collectively, point toward more complicated scenarios for seizure generation than simple excessive excitation and insufficient inhibition and suggest that different seizure onset patterns can arise from different cellular and synaptic mechanisms. ## Materials and Methods ### Intracranial electroencephalogram In a previous work [36], approved by the Centre Hospitalier de l’Université de Montréal (CHUM) Ethics Committee, we reviewed 103 iEEG ictal recordings using intracranial EEG electrodes (128- channel recording system, Harmonie, Stellate, Montreal, Canada, 2000Hz sampling rate) from 20 patients with drug-resistant focal epilepsy who were candidates for resective surgery and, identified nine common seizure onset patterns [7,8]. In the current study, we used some of those seizure onset classes to validate and classify our model simulations. ### Computational model In the study of epilepsy, the computational approach has been widely utilized in order to provide a mechanistic insight into mechanisms responsible for seizure generation [37,38]. In general, there are two main approaches in constructing mathematical models of ictal activity: biophysically realistic models and phenomenological models. Biophysically realistic modeling approaches rely on detailed cellular and synaptic mechanisms which lead to highly parameterized and complicated models [39–43]. On the other hand, phenomenological models describe the dynamic without reference to mechanisms and are generally more tractable and computationally less expensive [44,45]. These models, while being straightforward in implementation, do not have a physiological equivalent in their elements. Therefore, using phenomenological models, findings are disconnected from physiological phenomena, rendering interpretation more limited when designing therapeutic measures. More recently a hybrid approach has been introduced that combines both biophysically realistic and phenomenological models and is aimed at capturing the benefits of these two approaches through the incorporation of the relevant biophysical features into a simple framework [46]. In this study, we adopted the hybrid modeling approach. We used a well-established coupled neuronal population model, the Jansen model [47], and expanded it with some relevant biophysical features. We incorporated the slow and fast inhibitory mechanisms as suggested by Wendling et al. [35] to respectively account for the dendrite and soma targeting inhibitory interneurons. As presented in Fig.1a, the model consists of four components: the pyramidal cells (PYR), the excitatory interneurons (E), the slow dendrite targeting inhibitory interneurons (ID) and the fast soma targeting inhibitory interneurons (IS). The postsynaptic activity of the pyramidal cells (the main excitatory component) is considered as the model output to simulate the ictal iEEG recordings. ![Figure 1.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/17/2023.02.16.23286033/F1.medium.gif) [Figure 1.](http://medrxiv.org/content/early/2023/02/17/2023.02.16.23286033/F1) Figure 1. Excitatory inhibitory neuronal populational model and the activation functions. (a) The neuronal population model which is composed of excitatory and inhibitory components: pyramidal cells (PYR), excitatory interneurons (E), dendrite targeting inhibitory interneurons (ID) and soma targeting interneurons (IS). (b) The nonlinear activation function of a neural population that relates the average postsynaptic potential to an average firing rate for that population. The dotted black curve represents the commonly used sigmoid function and the solid black curve shows the activation function we used which considers depolarization block. (c) Decreasing the depolarization block threshold for the inhibitory populations. In this figure, the red activation function has a lower depolarization block threshold (θ = 0) compared to the black activation function (θ = 4). In this approach the average postsynaptic potential of each neural population is modeled by two main parts: 1) an activation function that describes a nonlinear relationship between the neural population’s average postsynaptic potentials and its output firing rate and 2) a linear transfer function, ![Graphic][1], that transforms the input firing rate into average postsynaptic potential. In this equation, W and τw respectively represent the average synaptic gain and average synaptic time delays for each neural population. The model can be expressed as a system of coupled second-order differential equations [35]: ![Formula][2] ![Formula][3] where *A, B* and *G* represent the synaptic gains of respectively the excitatory, dendrite targeting inhibitory interneurons and soma targeting inhibitory interneurons. P(t) describes inputs from other regions which are modeled as Gaussian white noise similar to the original model [47]. Parameters *C**1* to *C**7* are the connectivity constants which represent the average number of synaptic contacts (Table 1). In Eq.2, Fj where j∈ {e,s,d}, represents the activation function for the excitatory population when *j* = *e* and for dendrite and soma targeting inhibitory interneurons when respectively *j* = *d* and *j* = *s*. The variable *θ**j* represents the depolarization block threshold for dendrite targeting and soma targeting inhibitory interneurons respectively when *j* = *d* and *j* = *s*. View this table: [Table 1.](http://medrxiv.org/content/early/2023/02/17/2023.02.16.23286033/T1) Table 1. The model parameters used in Eq.1 and Eq.2. Usually, the activation function is assumed to follow a sigmoid-shaped function; however, more recent studies suggest a bell-shaped activation function which considers depolarization block [48–50]. An activity-driven depolarization block is suggested to be a rate-limiting mechanism to protect cells from excessive firing and can be due to synaptic resources or energy depletion during high levels of activity [49]. Accordingly, in this study, we replaced the commonly used sigmoid function of the original model with the firing rate function in Eq.2 to account for the effect of depolarization block in our model (Fig.1b) [50]. A second mechanism that was incorporated in the model was the collapse and recovery of the inhibitory activity. Some studies have reported seizure initiation due to the failure of inhibition with an important role for the depolarization block of inhibitory cells [28,51–54]. On the other hand, there is experimental evidence for the enhancement of inhibition before and at the very onset of seizure, contradicting the widely accepted paradigm on inhibitory failure mechanisms [29–33,55,56]. To explore both scenarios, we further expanded the model by considering variable depolarization block threshold for the activation function of the inhibitory populations (i.e., *θ**s*, *θ**d* in Eq.2) (Fig.1c). In all simulations, the depolarization block threshold of the excitatory populations was set to 15 (i.e., *θ**e*=15) which was chosen to adjust the excitatory activation function to have the same slope at half activation as the sigmoid function used by Jansen et al. [47] in their original model (Table 1). Inhibitory neurons due to their limited size become activated for smaller inputs and reach their depolarization block with a faster dynamic compared to the large pyramidal cells [48]. In order to consider these differences between the excitatory and inhibitory populations in our model, we modified parameters *r**i* and *v**i* as presented in Table 1 and chose *θ**d* and *θ**s* value to be 4 so that it would have the same slope at half activation as a sigmoid function with parameters *r**i* and *v**i*. However, wherever we intended to explore the effect of the collapse of an inhibitory population, we reduced the depolarization block threshold of that inhibitory population to zero (i.e., *θ**d*, *θ**s*). Fig.1c depicts the activation function of the inhibitory populations for both depolarization block threshold values equal to 0 and 4 with respectively the red and black curves. To account for different scenarios of disturbed balance between excitation and inhibition, we used a stochastic approach. In short, the model parameters for excitatory and inhibit synaptic gains (*A, B* and *G*) were divided into one-unit intervals and randomly sampled within each interval. Then, for each interval of the parameters, the most probable pattern of activity was identified. In the next step, starting from the range of *A, B* and *G* parameters, which generated background activity, we examined if the collapse and recovery of the inhibitory activity, which respectively were modeled through decreasing and increasing the parameters *θ**d* and *θ**s*, could result in the emergence of different seizure onset patterns. For the parameter *A*, we considered a range of values that varied from very low excitation levels to values inducing depolarization block in the excitatory population and explored the effect of different inhibitory modulations in transforming the background activity into ictal patterns; we did not explore “A” values beyond this level. ## Results We first validated our model by producing some of the commonly observed patterns of ictal activity. Using Eq.1 and Eq.2 for fixed values of *θ**e*, *θ**d* and *θ**s* and only by varying the synaptic gains *A, B*, and *G*, we could generate six ictal patterns observed in iEEG recordings as described below (Fig.2): 1. Low Voltage Fast Activity (LVFA): the initial amplitude is lower than 30mV with the frequency range of γ activity (Fig.2(b)), 2. Rhythmic spike/sharp waves: rhythmic spikes with the duration of 20ms to 70ms or sharp waves (duration of 70-200ms); because of the discernable similarity between these two patterns, we merged these two patterns in a single group to simplify our classification (Fig.2(c)), 3. Rhythmic Spike and Wave (rSW): spike and slow wave complexes with the frequency ranging between 2 and 4 Hz (Fig.2(d)), 4. High Amplitude Fast Activity (HAFA): the frequency is greater than 13 Hz and amplitude higher 30mV (Fig.2(e)), 5. Rhythmic α/β: rhythmic activity in the frequency range of α/β activity (Fig.2(f)), 6. Burst suppression: burst of rhythmic spikes interrupted with brief periods of voltage suppression/attenuation (Fig.3(g)). ![Figure 2.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/17/2023.02.16.23286033/F2.medium.gif) [Figure 2.](http://medrxiv.org/content/early/2023/02/17/2023.02.16.23286033/F2) Figure 2. Model simulations vs. real iEEG ictal recordings. The signals in blue are the simulated signals produced by our computational model and the signals in grey are the real depth iEEG recordings for (a) background, (b) Low Voltage Fast Activity (LVFA), (c) rhythmic slow spikes/sharp waves, (d) rhythmic Spike and Wave (rSW), (e) High Amplitude Fast Activity (HAFA), (f) rhythmic α/β and (g) burst suppression ictal activities. ![Figure 3.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/17/2023.02.16.23286033/F3.medium.gif) [Figure 3.](http://medrxiv.org/content/early/2023/02/17/2023.02.16.23286033/F3) Figure 3. The results from the Monte Carlo simulation. The colors represent the most probable pattern of ictal activity for each one-unit interval of the parameters A, B, and G with θd = θs = 4 for the left column and θd = θs = 0 for the right column. To identify the range of parameters producing each activity pattern, we first divided the range of synaptic gains *A, B* and *G* into one-unit intervals. Then, taking the Monte Carlo approach, we drew 100 random values for *A, B* and *G* within each interval and ran the simulations for both low and high values of depolarization block threshold of inhibitory populations (*θ**d*, *θ**s*). In the next step, based on the time and frequency domain features of each simulated signal we classified them into six distinct activity patterns that the model could generate as displayed in Fig.2. We finally identified the range of parameters *A, B* and *G* that generated each pattern of activity for both low and high values of *θ**d* and *θ**s* (Fig.3) The result of our Monte Carlo simulations is presented in Fig.3 and shows that the background activity (shown in blue) was generated for both low and high excitation levels (i.e., *A* values). For the mid-range excitation values, non-background patterns of activity emerge (i.e., ictal activities). However, comparing the left and right columns we can see that for high values of *θ**d* and *θ**s*, due to higher inhibitory effect, the values of *A* that generate ictal activities (i.e., non-background activity) are less restricted compared to *A* values that produced ictal patterns when low values of *θ**d* and *θ**s* are in place. To explore conditions in which the background activity can turn into an ictal pattern, we started from the background activity (i.e., the blue region in Fig.3) with both low and high *A* values (i.e., low and high levels of excitation) and tested as to whether collapse or recovery of the inhibitory activities (modeled by altering *θ**d* and *θ**s*) could result in seizure initiation and generation of different seizure onset patterns. ### Seizure generation with low preictal excitation For values of *A* lower than 3, both low and high values of *θ**d* and *θ**s* resulted in background activity –only for 2