<figure><img src="image_1.png" /></figure>

HAL
open science

# A Model for the Propagation of Seizure Activity in Normal Brain Tissue

Damien Depannemaecker, Mallory Carlu, Jules Bouté, Alain Destexhe

► To cite this version:

Damien Depannemaecker, Mallory Carlu, Jules Bouté, Alain Destexhe. A Model for the Propagation of Seizure Activity in Normal Brain Tissue. eNeuro, 2022, 9 (6), pp.ENEURO.0234-21.2022. 10.1523/ENEURO.0234-21.2022. hal-04174848

HAL Id: hal-04174848

https://hal.science/hal-04174848

Submitted on 4 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Model for the Propagation of Seizure Activity in Normal Brain Tissue

Damien Depannemaecker,*, Mallory Carlu,† Jules Bouté, and Alain Destexhe

https://doi.org/10.1523/ENEURO.0234-21.2022

French National Centre for Scientific Research (CNRS), Paris-Saclay Institute of Neuroscience (NeuroPSI), Paris-Saclay University, 91198 Gif sur Yvette, France

## Abstract

Epilepsies are characterized by paroxysmal electrophysiological events and seizures, which can propagate across the brain. One of the main unsolved questions in epilepsy is how epileptic activity can invade normal tissue and thus propagate across the brain. To investigate this question, we consider three computational models at the neural network scale to study the underlying dynamics of seizure propagation, understand which specific features play a role, and relate them to clinical or experimental observations. We consider both the internal connectivity structure between neurons and the input properties in our characterization. We show that a paroxysmal input is sometimes controlled by the network while in other instances, it can lead the network activity to itself produce paroxysmal activity, and thus will further propagate to efferent networks. We further show how the details of the network architecture are essential to determine this switch to a seizure-like regime. We investigated the nature of the instability involved and in particular found a central role for the inhibitory connectivity. We propose a probabilistic approach to the propagative/non-propagative scenarios, which may serve as a guide to control the seizure by using appropriate stimuli.

**Key words:** computational modeling; epilepsy; inhibitory population; seizure control; seizure propagation; spiking network model

## Significance Statement

Our computational study investigates how epileptic activity invades normal brain tissue, and shows the specific role of the inhibitory population, and its dynamical and structural aspects, using three different neuronal networks. We find that both structural and dynamic aspects are important to determine whether seizure activity invades the network. We show the existence of a specific time window favorable to the reversal of the seizure propagation by appropriate stimuli.

## Introduction

Epilepsy is one of the most common neurologic disorders (Beghi, 2020) and can take numerous forms. It is associated with the presence of paroxysmal electrophysiological events and seizures, usually recorded in humans using the electroencephalogram (EEG). However, EEG recordings do not allow us to probe the activity of single neurons within the network. More recently, the recording conducted with microelectrode arrays made it possible to obtain spike information on the order of a hundred neurons in human epileptic patients

(Peyrache et al., 2012; Dehghani et al., 2016; Paulk et al., 2022).

Such microelectrode recordings showed that neuronal activity during seizures does not necessarily correspond to synchronized spikes over the whole neuron population, as previously modeled (Soltesz and Staley, 2008), including models at different scales from cellular to whole-brain levels (Depannemaecker et al., 2021, 2022). In fact, it turns out that the dynamics of neural networks during

<footnote>Received May 19, 2021; accepted September 29, 2022; First published November 2, 2022.
The authors declare no competing financial interests.</footnote>

<footnote>Author contributions: D.D., M.C., and A.D. designed research; D.D., M.C., and J.B. performed research; D.D. and M.C. contributed unpublished reagents/analytic tools; D.D., M.C., and J.B. analyzed data; D.D. and M.C. wrote the paper.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Examples of inhibitory recruitment during seizures. **a**, Raster plot of three different seizures from the same patient, 92 neurons were identified, 24 putative inhibitory cells (red) and 68 putative excitatory cells (green). **b**, Corresponding firing rate of the putative inhibitory population (red) and the putative excitatory population (green). A plateau of high activity of the putative inhibitory cells can be observed during the seizure (highlighted in dashed purple oval). This was done with data from the study by Dehghani et al. (2016).</figcaption></figure>

seizures are more complex (Jiruska et al., 2013) and still poorly understood. In particular, it is not known how the paroxysmal activity of the seizure propagates, driving other networks into seizure activity.

Here, we investigate this problem using computational models. As a starting point, we consider examples of seizures where the inhibitory network is strongly recruited, while the firing of excitatory cells is diminished. Figure 1 shows three seizures from a patient that were recorded, using Utah-arrays, before resection surgery in a case of intractable epilepsy. From these intracranial recordings, 92 neurons have been identified and isolated and were classified into the following two groups: fast-spiking (FS) neurons and regular-spiking (RS) neurons, based on spike shape, autocorrelograms, firing rates and cell-to-cell interactions (Peyrache et al., 2012). Remarkably, direct cell-to-cell functional interactions were observed, which demonstrated that some of the FS cells are inhibitory while some of the RS cells are excitatory (for details, see Peyrache et al., 2012). The three seizures shown in Figure 1 were taken from the analysis of Dehghani et al. (2016; see this article for details) and are shown with the firing rate of each population of cells. During the seizure, we can observe a plateau of high activity of FS cells, and a strongly reduced activity of RS cells. This phenomenon of unbalanced dynamics between RS and FS cells was only seen during seizures in this patient (Dehghani et al.,

2016). It shows that, in these three examples, the seizure was manifested by a strong “control” by the inhibitory FS cells, which almost silenced excitatory RS cells. Interestingly, a very different conclusion would have been reached if no discrimination between RS and FS cells were performed, which underlies the importance of discriminating RS and FS cells for a correct interpretation of the dynamics during seizures. Based on such measurements, we built computational models based on a larger number of cells to consider network effects that are not directly accessible with the recordings. We were interested in how seizure activity propagates or not, and in the determinants of such propagation.

The region of the brain where the seizure starts is called the seizure focus, although in certain patients it can be distributed over several foci (Nadler and Spencer, 2014), then the seizure spreads to other regions of the brain. When another such region is reached, it can in turn be driven into seizure activity, in which case the seizure activity propagates. It can also control it (as in Fig. 1), in which case the seizure would remain confined to a more restricted brain region.

To gain an understanding of the dynamics underlying these two scenarios, we study the response of networks using three different neuron models (adaptive exponential integrate and fire (AdEx), conductance-based adaptive exponential integrate-and-fire (CAdEx), and Hodgkin-Huxley (HH) models), interacting through conductance-based synapses, to an incoming paroxysmal (seizure-like) perturbation. We observe two types of behavior, which we represent in Figure 2: one where the incoming perturbation successfully increases the activity of the excitatory population, thus making it stronger than the input; and the other where only the inhibitory population strongly increases its activity, thus controlling the perturbation. In the first case, where the excitatory population discharges very strongly, it is therefore likely to transmit, or even amplify, the perturbation transmitted to the next cortical column. We have therefore called this situation the

This work was funded from the European Union's Horizon 2020 Framework Programme for Research and Innovation under the Specific Grant Agreement No. 945539 (Human Brain Project SGA3) and the Centre National de la Recherche Scientifique (France).

*D.D. and M.C. contributed equally to this work.

Correspondence should be addressed to Damien Depannemaecker at damien.depannemaecker@univ-amu.fr or Alain Destexhe at alain destexhe@cnrs.fr.
https://doi.org/10.1523/ENEURO.0234-21.2022

Copyright © 2022 Depannemaecker et al.

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.
<figure><img src="image_4.png" /><figcaption>Figure 2. Cartoon of the modeled scenarios.</figcaption></figure>

propagative scenario. In the opposite case, where the firing rate of the excitatory population remains much lower than the perturbation, the seizure-like event will not spread to the neighboring region, we therefore call this situation the non-propagative scenario. We then propose a more precise approach, based on the AdEx network, that mixes structural and dynamic ingredients to unravel key aspects of the mechanisms at play. Focusing on the different input connectivity profiles for each node in the network, we are able to build separate groups of neurons that display significantly different dynamics with respect to the perturbation. Finally, we study the possibility of a proactive approach, based on the application of an extra stimulus with the aim of reversing the propagative behavior, thus controlling the spread of the seizure.

# Materials and Methods

## Computational models

We use for this study a mathematical model of electro-physiological activity based on ordinary differential equations, describing the dynamics of the membrane potential of the neurons through their interactions.

Each neuron model in the network is described by Equations 1 and 2, the AdEx model (Brette and Gerstner, 2005; Naud et al., 2008), as follows:

$$ 
\begin{aligned}
C \frac{dV}{dt} &= g_L(E_L - V) + g_L \Delta_T \exp\left(\frac{V - V_T}{\Delta_T}\right) - w + I_{\text{syn}} \\
\tau_w \frac{dw}{dt} &= a(V - E_L) - w.
\end{aligned}
\quad (1) $$

Where $I_{\text{syn}}$ is the synaptic input. When the membrane potential crosses a threshold, a spike is emitted, and the system is reset as follows:

$$ \text{if } V \geq V_D \text{ then } \begin{cases} V \to V_R \\ w \to w + b \end{cases}. \quad (2) $$

Parameters used for the excitatory (RS) and inhibitory (FS) populations are respectively as follows: $V_t = -50$ mV and $V_t = -48$ mV, $\Delta_t = 2$ mV and $\Delta_t = 0.5$ mV, $b = 100$ pA and $b = 0$ pA, and $\tau_w = 1000$ ms for RS. For both populations: $C = 200$ pF, $g_L = 10$ nS, $E_L = -65$ mV, $a = 0$ nS, $V_{\text{reset}} = -65$ mV, $t_{\text{refractory}} = 5$ ms.

To compare some of the results obtained with the AdEx model, we used two other models of neuronal activity. First, the CAdEx model, which solves some of the limitations of the AdEx model (Górski et al., 2021). The equations read:

$$ 
\begin{aligned}
C \frac{dV}{dt} &= g_L(E_L - V) + g_L \Delta_T \exp\left(\frac{V - V_T}{\Delta_T}\right) + g_A(E_A - V) + I_{\text{syn}} \\
\tau_A \frac{dg_A}{dt} &= \frac{g_A}{1 + \exp\left(\frac{V_A - V}{\Delta_A}\right)} - g_A.
\end{aligned}
\quad (3) $$

When the membrane potential crosses a threshold, a spike is emitted, and the system is reset as in:

$$ \text{if } V \geq V_D \text{ then } \begin{cases} V \to V_R \\ g_A \to g_A + \delta g_A \end{cases}. \quad (4) $$

Parameters used for inhibitory (FS) populations are as follows: $g_L = 10$ nS, $E_L = -65$ mV, $V_T = -50$ mV, $g_A = 0.1$ nS, $E_A = -70$ mV, $\delta g_A = 0$ nS, $C = 200$ pF, $\Delta_T = 0.5$ ms, $V_A = -45$ mV, $I_{\text{syn}} = 0.0$ nA, $t_{\text{refractory}} = 5$ ms, $V_{\text{reset}} = -65$ mV, $\tau_A = 0.01$ ms, $\Delta_A = 0.5$ mV and for the excitatory (RS) $g_L = 10$ nS, $E_L = -65$ mV, $V_T = -50$ mV, $\delta g_A = 1$ nS, $E_A = -65$ mV, $\delta g_A = 1$ nS, $C = 200$ pF, $\Delta_T = 2$ mV, $V_A = -30$ mV, $I_{\text{syn}} = 0.0$ nA, $t_{\text{refractory}} = 5$ ms, $V_{\text{reset}} = -65$ mV, $\tau_A = 1.0$ s, $\Delta_A = 1$ mV Then we use the HH model (Hodgkin and Huxley, 1952), with the following equations:

$$ C_m \frac{dV}{dt} = -g_L(E_L - V) - g_K n^4 (V - E_K) - g_{Na} m^3 h (V - E_{Na}) + I_{\text{syn}}, \quad (5) $$

with gating variables (in ms):

$$ 
\begin{aligned}
\frac{dn}{dt} &= \frac{0.032(15. - V + V_T)}{\left(\exp\left(\frac{15. - V + V_T}{5.}\right) - 1\right)} \\
&\quad \times (1. - n) - 0.5 \exp\left(\frac{10. - V + V_T}{40.}\right) n \\
\frac{dh}{dt} &= 0.128 \exp\left(\frac{17. - V + V_T}{18.}\right) \\
&\quad \times (1. - h) - \frac{4.}{1 + \exp\left(\frac{40. - V + V_T}{5.}\right)} h \\
\frac{dm}{dt} &= \frac{0.32(13. - V + V_T)}{\left(\exp\left(\frac{13. - V + V_T}{4.}\right) - 1\right)} \\
&\quad \times (1 - m) - \frac{0.28(V - V_T - 40.)}{\left(\exp\left(\frac{V - V_T - 40.}{5.}\right) - 1\right)} m,
\end{aligned}
\quad (6) $$

with $C_m = 200$ pF, $E_L = -65$ mV, $E_{Na} = 60$ mV, $E_K = -90$ mV, $g_L = 10$ nS, $g_{Na} = 20$ nS, $g_K = 6$ nS, $V_{\text{Texc}} = -50$ mV, $V_{\text{Tinh}} = -52$ mV.

For all types of neuron models, the parameters have been chosen in biophysical range (Hodgkin and Huxley,
1952; Hille, 1992; Naud et al., 2008; Górski et al., 2021) to keep the basal asynchronous irregular (AI) activities (Brunel, 2000) into a range of firing rates coherent with experimental observations (El Boustani et al., 2007; Destexhe, 2009; Zerlaut et al., 2018).

The network is built according to a sparse and random (Erdos–Renyi type) architecture where a fixed probability of connection between each neurons is set to 5%. We consider a network model of 10,000 neurons, built according to specific properties of the cortex. This network is made of an inhibitory (FS) and an excitatory (RS) population, respectively, representing 20% and 80% of the total size of the system as previously done in the study by Carlu et al. (2020). The communication between neurons occurs through conductance-based synapses. The synaptic current is described by the following:

$$I_{\text{syn}} = g_{E}(E_{E} - V) + g_{I}(E_{I} - V), \quad (7)$$

where $E_E = 0$ mV is the reversal potential of excitatory synapses and $E_I = -80$ mV is the reversal potential of inhibitory synapses. $g_E$ and $g_I$ are, respectively, the excitatory and inhibitory conductances, which increase by quantity $Q_E = 1.5$ nS and $Q_I = 5$ nS for each incoming spike. The increment of conductance is followed by an exponential decrease according to the following:

$$\frac{dg_{E/I}}{dt} = -\frac{g_{E/I}}{\tau_{\text{syn}}}, \quad (8)$$

with the following:

$$\tau_{\text{syn}} = 5 \text{ ms.}$$

The network thus formed receives an external input, based on the activity of a third population (excitatory) of the same size as the excitatory population. Each of its neurons is connected to the rest of the network according to the same rule as mentioned earlier (fixed probability, 5% for each connection). This external population produces spikes with a Poissonian distribution at a given tunable rate. The external perturbation that mimics the incoming seizure occurs through the augmentation of this firing rate.

The shape of the latter is described by the following:

$$\nu_{\text{pert}}(t) = \beta + \alpha \times \left( \frac{\exp(-(t - T_1)^2}{2 \times \tau_{\text{on}}^2}) \times H(-(t - T_1)) + H(-(t - T_2)) \times H(t - T_1) + \exp(-(t - T_2)^2}{2 \times \tau_{\text{off}}^2}) \times H(t - T_2) \right), \quad (9)$$

where $H$ is the Heaviside function and $\beta = 6$ Hz is the basal constant input. This function takes the general form of a high plateau, where $T_1$ and $T_2$ are the times when the perturbation reaches its beginning and end, respectively, and $\alpha$ defines its maximal height. $\tau_{\text{on}}$ and $\tau_{\text{off}}$ are, respectively, time constants associated with the exponential rise and decay of the perturbation.

For all three types of networks, it is possible to have different connectivities (i.e., different set of random connectivities) and realizations of a Poisson drive (i.e., the generator of the Poisson noise can vary). It is also possible

to fix the seed of the noise either for the connectivity or for the Poisson drive (or both) to analyze specific conditions.

We create network connectivities by allowing a 5% chance of connection between any two neurons, which will indeed lead to an average of 5% of connection, but with some variation. Some neurons can have more after connections from inhibitory neurons than others, which will make them more inhibited, and the same goes with excitatory connections, creating a variation between neurons because of the random nature of the network.

## Coarse graining and continuous analysis

To analyze in details what mechanisms are at play in the network during a seizure-like event, we resort to a combination of two methods: a so-called “structural coarse graining”; that is, we gather neuron models in n groups according to their inhibitory in-degree (the number of inhibitory connections they receive, as introduced before, and we study their time evolution through statistics of their membrane potential (mean and alignment) over these groups. In other words, at each integration time step, we will obtain n values of mean membrane potentials, one for each group, as well as n values of the Kuramoto order parameter (measuring alignment in groups).

To obtain the Kuramoto order parameter, we first transform the single-neuron membrane potentials into phase variables by applying a linear mapping $v_j \in [V_R, V_D] \to \theta_j^v[0, \pi]$. Then the Kuramoto order parameter is computed through the following equation:

$$R \exp(i\Psi) = \frac{1}{N} \sum_{j} \exp(i\theta)_{j}^{v}. \quad (10)$$

$R \in [0, 1]$ gives the degree of “alignment” (if it persists in time, one would say synchronization): $R=1$ implies full alignment, while $R=0$ implies no alignment whatsoever. $\Psi \in [0, \pi]$ tells us the mean phase of the transformed variables (directly related to the mean membrane potential).

Let us mention one caveat here. The membrane potentials are not mapped on the full circle to avoid artificial periodicity of the obtained angles: having $V = V_R$ is not the same as having $V = V_D$. One may thus ask why such a measure is used instead of the usual measures of dispersion such as the Standard deviation (SD). We use the Kuramoto order parameter because it gives a naturally normalized quantity, thus allowing a direct comparison of what is happening at each time step.

## Data availability

The code/software described in the paper is freely available online at https://github.com/HumanBrainProject/PropNoProp.

## Results

We start by showing how, in networks of various neuron models, a paroxysmal external stimulation can trigger a
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 3.** *a–f*, Firing rate of the network populations in response to a perturbation (blue, the incoming perturbation; green, excitatory and in red inhibitory populations): propagative and non-propagative scenarios (respectively, left and right columns) for AdEx model (*a*, *b*), with amplitude of perturbation *α* = 80 Hz and *τ*<sub>on/off</sub> = 100 ms; CAdEx model (*c*, *d*) with *α* = 70 Hz and *τ*<sub>on/off</sub> = 80 ms; HH model with *α* = 60 Hz and (*e*), and with *α* = 140 Hz and *τ*<sub>on/off</sub> = 60 ms (*f*). For each model, the networks are the same in the propagative or non-propagative scenarios; the only difference comes from the incoming input with different realizations.</figcaption></figure>

seizure, depending on various parameters. We show how
the dynamics can differ from model to model and in their
common features. Then, we propose a structural analysis
based on the mean firing rates of individual neuron mod-
els to guide a particular coarse graining, which we use as
a filter to observe the dynamics and gain further under-
standing, from both qualitative and quantitative per-
spectives. Finally, we show how this study can guide a
proactive approach to reduce the chances of seizure
propagation.

**Propagative and non-propagative scenarios**

Throughout this study, we assume that the networks
depicted in the previous section represent a small cortical
area receiving connections from an epileptic focus.
Specifically, the arrival of the seizure is modeled by a sud-
den rise in the firing rate of the external (afferent) Poisson
input from the region where the seizure originates. In
other words, we are not concerned with how seizures
*originate* (epileptogenesis), but how they can propagate.
Therefore, we will frame our analysis into two main sce-
narios: propagative (i.e., the network develops an excita-
tory firing rate greater than the input), which makes it able

to propagate the seizure to efferent regions; and non-
propagative, where the excitatory firing rate is lower
than the input, thus attenuating the incoming signal. As
described in the Materials and Methods section, the
perturbation starts with an exponential growth followed
by a plateau and ends with an exponential decrease,
going back to the basal level (Fig. 3, blue curves). We
show in Figure 3 the response of the various networks
to this type of perturbation.

Here we can distinguish between two classes of
macroscopic differences between propagative and non-
propagative scenarios.

In the first class (for AdEx and CAdEx), the difference is
binary, which means that the network either features a
very strong increase in the firing rate of the inhibitory and
excitatory populations, or that the sharp increase in the
firing rate only concerns the inhibitory population, thus
strongly limiting the activity of the excitatory population
(consequently preventing the seizure from spreading to
the next region). From this perspective, the propagative
scenario can be understood as a loss of balance be-
tween excitatory and inhibitory firing rates, which the
network struggles to find once the excitatory population
has reached very high firing rates. Interestingly, these
<figure><img src="image_11.png" /><figcaption>**Figure 4.** Grid search on the amplitude and slope of the incoming perturbation for each network. **a, b,** The percentage of realizations that propagate (Prop.), respectively, for AdEx (a) and CAdEx (b) networks. **c, d,** For HH networks, the means and SDs (over realizations) of the difference in firing rates between excitatory (c) and Poisson (d) populations ($\Delta$ firing rate = $v_e - v_{Pois}$), averaged over the length of the plateau.</figcaption></figure>

two scenarios can occur for the same global shape of the perturbation but changing only the noise and network realizations. It must be noted that the 200 Hz maximum frequencies measured here are the results of the temporal binning of the global spiking dynamics, taken as $T=5$ ms, which corresponds to the refractory time of the single neurons in AdEx and CAdEx. Upon choosing a shorter binning (e.g., $T=1$ ms), higher-frequency peaks are observed, going up to 800 Hz, thus hinting at overall faster dynamics.

In the second class (HH), there is a rather continuous difference between propagative and non-propagative scenarios, as can be seen in Figure 4c, depending on the amplitude of the perturbation.

## Influence of the perturbation shape

To study how the shape of the perturbation affects the networks response, in Figure 4 we screened different time constants of the exponential growth rates and maximum amplitudes of the plateau with 100 seeds (for both network and noise realizations for each couple of values) and probed, in the case of AdEx (Fig. 4a) and CAdEx (Fig. 4b), the number of realizations, which yield propagative behavior, as it shows binary possible scenarios. Meanwhile, in the HH case, the perspective is a little different: we chose to show two figures, displaying means and SDs over realizations of the difference in firing rate between excitatory and Poisson populations (averaged over the plateau), $\Delta_{firing rate} = v_e - v_{Pois}$ (Fig. 4c,d, respectively). As can be expected, for all networks (AdEx, CAdEx, and HH) the amplitude of the perturbation plays a determinant role in the type of scenario we eventually find (propagative or not), although, in opposite directions and of different natures. Indeed, for both AdEx and CAdEx, increasing the amplitude increases the chance of having a propagative scenario for a fixed slope, in a binary fashion, while in the case of HH the contrary is observed, and in a continuous fashion. Also, we observe a slight coupling effect between

slope and amplitude: for higher amplitudes, the propagation range extends to slower perturbations. On the contrary, in the HH network, it seems that the slope does not play any major role, hinting at a much less dynamic effect: the difference manifests itself as local equilibria of the networks under consideration that were reached no matter the time course. Moreover, the SDs, in addition to showing no clear dependence on either amplitude or slope, are very small compared with the means, thus providing evidence that noise also does not play any significant role here. These observations highlight once again the deep differences between the two types of networks and their respective phenomenology.

Interestingly, in the case of AdEx and CAdEx, there exists a limited, bistable region here, ~80 Hz, where the perturbation may or may not propagate in the network, depending on the noise realization. Thus, the scenario does not trivially depend on the amplitude and time constants of the perturbation in this region, which makes the latter a perfect test bed to study more deeply the internal mechanisms at play, and will thus be the main focus of the remainder of this article.

## Influence of structural aspects on the dynamics

In the following, we turn our attention to the bistable region of AdEx networks, where the two scenarios are present, and investigate the origin of the source of the divergence. There are two main differences between the simulations under consideration: the realization of the network connectivity and the realization of the external input, as both rely on random number generators. We have therefore successively fixed each of them and observed that the two behaviors were still present. Also, the global scenarios were indistinguishable from those showed so far. First, this allows us to fix the network connectivity (which will become determinant in this part) without losing the richness of the phenomenology. Second, this tells us that what shapes the
<figure><img src="image_12.png" /><figcaption>**Figure 5.** Influence of connectivity on single neurons firing rates. **a**, Influence of Poissonian ($N_{\text{inp}}^{\text{Pois}}$), excitatory ($N_{\text{inp}}^{\text{Exc}}$), and inhibitory ($N_{\text{inp}}^{\text{Inh}}$) in-degree on the firing rates of excitatory neurons ($\nu_{E}^{\text{NP}}$), and inhibitory neurons ($\nu_{I}^{\text{NP}}$) in the non-propagative scenario of the AdEx network. The standard Pearson correlation coefficient $\rho$ is estimated. **b**, Time-averaged single-neuron firing rates and differences in propagative versus non-propagative regimes, as a function of both inhibitory and Poissonian in-degrees.</figcaption></figure>

distinction between the two phenomena is more complex than a simple question of structure, or realization of the input. Another perspective is then needed to explore the internal dynamics of the network in both scenarios. As the models into consideration have very large number of dimensions, as well as quite intricate structures, brute force analytical approaches are simply not conceivable.

Let us then take a step back and investigate the relationship between the firing rate of each neuron and its number of afferent input connections for the three kinds of input (inp): excitatory ($N_{\text{inp}}^{\text{Exc}}$), where Exc is excitatory; inhibitory ($N_{\text{inp}}^{\text{Inh}}$), where Inh is inhibitory; and Poissonian ($N_{\text{inp}}^{\text{Pois}}$), where Pois is Poissonian. Figure 5a shows the average (firing rates ($\nu_{E}^{\text{NP}}$ and $\nu_{I}^{\text{NP}}$, where NP is non-propagative) measured over the whole non-propagative scenario
<figure><img src="image_13.png" /><figcaption>**Figure 6.** Grid search on the in-degree inhibitory probability of connection for the AdEx network. **a**, **b**, Percentage of propagation (Prop.) with parameters, as follows: α = 70 Hz and τ = 70 ms (**a**); and α = 95 Hz and τ = 70 ms (**b**), where for both figures Pii is the probability of connection from inhibitory neurons to in-hibitory neurons. Decreasing the probability of connection from inhibitory to excitatory neurons or increasing the probability of con-nection from inhibitory to inhibitory neurons tends to decrease the overall inhibition in the network and thus facilitates propagative behavior.</figcaption></figure>

for each neuron in the AdEx network (simply defined as the total number of spikes divided by the total integration time, after having discarded a transient), plotted as a function of the three different connectivities.

Note here that averaging over simulations for the sake of robustness might be a delicate matter, as we might lose constitutive differences in the process. As we are dealing with highly variable situations, we have to make compromises between generalizability and relevance. Therefore, we start with a single realization to then guide larger and more systematic investigations.

Interestingly, we see a much stronger influence coming from the inhibitory in-degree than from the Poissonian and excitatory ones. Counterintuitively, it even seems that excitatory in-degree has almost no effect at all on total measured firing rates. Indeed, from the point of view of Pearson's correlation, inhibitory in-degree is much more (anti)-correlated with the firing rate than the excitatory in-degree (almost no correlation) or the Poissonian in-degree (little correlation). Note that we observe the same structure for propagative scenarios (results not shown). Based on these results, we can analyze whether the most salient in-degrees (inhibitory and Poissonian) has any influence on the difference between propagative and non-propagative scenarios (Fig. 5b). Here, we see that the global dependency of the average single-neuron firing rates on inhibitory and Poissonian connectivity does not qualitatively change between propagative and non-propagative scenarios. However, the differences $\nu^P - \nu^{NP}$ display an inverse dependency on both variables: despite maintaining a qualitative similarity between first and second columns, the seizure tends to compensate the initial disparity of firing rates. In other words, the neurons that are initially less firing, because of their structural properties, are the most impacted by the seizure. Furthermore, it must be noted that, although there is no correlation between inhibitory and Poissonian in-degrees (as can be expected from random connectivities), the third column highlights that they both play a role in the single-neuron long-term dynamics.

To further understand the effect of the inhibitory connectivity, we chose two points from Figure 4a, one known to be always non-propagative, with $\tau = 70$ ms and $\alpha =$

70 Hz, and the other to be always propagative, with $\tau = 70$ ms and $\alpha = 95$ Hz. In both situations, we varied the probability of connection from the inhibitory population between $p = 0.04$ and $p = 0.06$, as shown in Figure 6. Note that the influence of the incoming inhibitory connectivity shifts the boundary between propagative and non-propagative behaviors. This is an important influencing factor in relation to the shape of the perturbation and in particular its amplitude.

It is worth pointing that these results establish a clear link between structure and dynamics, but structure is by itself not a sufficient criterion to understand the underlying mechanisms. We therefore focus on the temporal evolution of the propagating activity.

Beforehand, we take a step back and probe whether the differences in the individual mean firing rates are associated with specific roles in the dynamics. To achieve so, we start classifying, for the AdEx network, the indices of neurons in the raster plot according to the total number of spikes they emit during the whole simulation. We chose for this purpose a representative propagative scenario. The sequence of propagation of the perturbation then appears visually in Figure 7a. We observe, in the case of propagation, a fast cascade (consistent with the experimental observations; Neumann et al., 2017): some neurons are quickly driven into a sequence at the onset of the seizure. In addition, there is no perfect synchronization of the action potentials of all neurons. This is an interesting result, coherent with experimental observations on epilepsy (Jiruska et al., 2013).

Second, we examine the same situation, but sorting neuron indices as a function of the number of inhibitory inputs they receive, as shown Figure 7b, as it is the most influential structural feature we observed in our model. Here too, the cascade phenomenon is clearly visible, indicating that the inhibitory input connectivity has a central influence on the dynamics at play during the perturbation in the propagative scenario.

Figure 8 shows the same pictures for CAdEx and HH networks. We see here that the behaviors of the CAdEx network are very similar to AdEx: sorting with firing rate or inhibitory in-degree gives very similar structures, and we can also distinguish here the cascading effect at the onset
<figure><img src="image_14.png" /><figcaption>**Figure 7.** Dynamics in the propagative scenario (AdEx). **a**, In a raster plot of a simulation with propagative behavior, neuron indices are sorted according to the number of spikes during the simulation. A “cascade” phenomenon can be observed when zooming on the onset of the perturbation propagation in the excitatory population. **b**, The same cascade phenomenon is observed when neuron indices are sorted as a function of the number of inhibitory inputs received. Note that the absence of excitatory activity after the perturbation is because of a strong adaptation current (Eqs. 1, 2).</figcaption></figure>

of the perturbation, following the indices. HH networks show quite a different phenomenology. First, the two sorting do not show the same structures, which hints at a more subtle mapping between inhibitory in-degree and long-term single-neuron model dynamics. In the firing rate sorting, we can still distinguish three blocs of distinct activity, and thus of populations, corresponding to the three key periods of the simulation: before stimulation, at the onset, and during the stimulation. Interestingly, it seems that before and during the stimulation different populations of neuron models are distinctly mobilized. While before the stimulation, the central neurons (with respect to their indices) are active, a double cascade contaminates the rest of them (toward higher and lower indices) at the onset, ending in a general surge of activity. This must be contrasted with the in-degree sorting panel, where the cascade is more unidirectional, as the main activity slides from low-connectivity indices (less connected) to the higher ones, until all neurons fire. This emphasizes the importance of the perspective chosen to analyze complex behavior: none of these perspectives alone completely explains the intricate interplay among structure, long-term, and short-term dynamics.

Altogether, these results show the relevance of adopting a perspective based on the inhibitory in-degree: it gives an operational method to rank single neurons, and this ranking is clearly associated with specific dynamic features, hence allowing us to study the role of the internal organization of the network before and during the paroxysmal event. As the cascade phenomenon is similarly visible in all types of networks, in the next section we focus

on the AdEx network. We push further this analysis by comparing propagative and non-propagative scenarios and make use of the continuous measures introduced in Materials and Methods.

## Continuous measures on subgroups of neurons

Focusing on the AdEx network, we first consider groups of neurons defined by their inhibitory in-degree. Note that these are somewhat artificial, as they are only statistical reflections of topological aspects of the network (i.e., there is no reason to think a priori that all neurons having *n* inhibitory inputs would have more privileged links among themselves than with those having a different number). However, they allow in principle a variable degree of categorization, based on the sampling of the inhibitory in-degree distribution, which eventually leads to different levels of (nonlinear) coarse graining (although we will consider only one such sampling here). Second, we switch our analysis to continuous variables, which allow a finer and more systematic analysis of the dynamics, as they do not depend on spike times. Indeed, although spike timings are the most accessible collective measures in real-life systems, which make them the most fitted candidates for “transferable” studies, we want here to take advantage of the virtues of mathematical modeling to probe the underlying mechanisms in these simulations, to then be able to draw conclusions on more accessible observables. We focus here uniquely on membrane potentials, as they are the closest proxy of the firing dynamics in the network and chose to use two main measures based on them: the mean $\mu_V$ and a modified Kuramoto
<figure><img src="image_15.png" /><figcaption>**Figure 8.** Dynamics in the propagative scenario (CAdEx and HH). a–d, Same plot as previously shown but for the CAdEx network (spike sorting, a; inhibitory in-degree sorting, b) and HH network (spike sorting, c; inhibitory in-degree sorting, d). Cascade phenomena are still observable in a, b, and d, hence showing its robustness, but not in c, where propagation takes a slightly different form, highlighting the contrast induced by different perspectives on a single-complex dynamics.</figcaption></figure>

order parameter *R*, which gives a naturally normalized measure of instantaneous alignment (or similarity) of the membrane potentials. Both are defined in time, over a class of neurons. As randomness plays a crucial role in our simulations, through network connectivity as well as noise realization, it is important to control how much it affects the results we obtain. To achieve so, we start by fixing the network connectivity while averaging over noise realizations, and then average over connectivities while looking for noise realizations that lead to propagative and non-propagative scenarios for each structure.

### Mean membrane potential in time

In Figure 9, a and b, the mean membrane potential $\mu_V$ defined for each group of excitatory (RS) and inhibitory (FS) neurons, in time. The top and bottom rows, respectively, refer to the averages and SDs over noise realizations (input), as the network connectivity is held fixed. For propagative (Fig. 9a) and non-propagative (Fig. 9b) scenarios, all the data presented from this point were obtained by regrouping neurons having the exact same inhibitory in-degree, thus corresponding to a discrete one-to-one sampling of the input distribution. Note that, given the network architecture under consideration, the number of afferent inhibitory synapses defined over both populations of neurons follows a binomial distribution

with a mean of ~100 connections. From that, we arrange neurons in groups of identical number of inhibitory connections, which gives us ~60 groups (varying with population and connectivities) containing at least one neuron.

To confirm that our results were not depending on the specific connectivity we had, we simulated 50 different networks with different connectivities (which otherwise were identical) and found a couple of noise realizations for each corresponding to propagative and non-propagative scenarios. Those various networks still have the same metastructure and follow similar statistics. They only show that within those specific choices, the variations that exist do not impact the results we show. We applied the same method to create the different groups, but the number of said groups could differ because of the random variability in the connections. Therefore, many “extreme” groups are poorly represented among the various connectivities, which would make them hard to average. We thus discarded them. The average and SD of $\mu_V$ over the 50 different connectivities is shown in Figure 10, a and b. The white lines could be a weak manifestation of the previous effect, which made the SD very high, coupled with the fixed range of color scales, imposed for the need of clarity. We observe that this figure looks very similar to Figure 9, which suggest that the results are not limited to a specific connectivity.
<figure><img src="image_16.png" /><figcaption>**Figure 9.** Mean membrane potential over subgroups of neurons (same network connectivity, different noise realizations) for each group defined as a function of their incoming inhibitory connections, averaged over 50 noise realizations (17 non-propagative and 33 propagative). a, b, Color maps correspond for each group to the average membrane potential (top) and SD (bottom) across noise realizations in the propagative situations (a) and non-propagative situations (b) for both excitatory (RS) and inhibitory (FS) populations. The blue rectangle highlights the (time) region where the system either switches to a propagative regime or remains stable.</figcaption></figure>

We see from Figures 9 and 10 that the inhibitory in-degree profile seems to play a major role in the overall dynamics. Indeed, as the perturbation is growing (starting 250 ms before the maximum at 2 s), we can first observe a strong increase of the mean membrane potential of all excitatory neurons, starting from low indices, then followed by a low-potential cascade,

also starting from low indices and then contaminating higher ones.

This latter effect is much clearer in the case of inhibitory neurons, where the cascade follows very well the input profile, in both propagative and non-propagative scenarios. Note that the low-potential area can be easily understood as a high-firing regime: neurons fire as soon as they
<figure><img src="image_17.png" /><figcaption>**Figure 10.** Mean membrane potential over subgroups of neurons (different network connectivities) for each group defined as a function of their incoming inhibitory connections. Here, we averaged over 50 network connectivities, for which we found a couple of noise realizations corresponding to propagative and non-propagative scenarios. a, b, Color maps correspond for each group to the average membrane potential (top) and SD (bottom) across different connectivities in the propagative (a) and non-propagative scenarios (b) for both excitatory (RS) and inhibitory (FS) populations. The blue rectangle highlight the (time) region where the system either switches to a propagative regime or remains stable.</figcaption></figure>

leave their resting potential, thus displaying very low values of membrane potential when calculated (and sampled) over time.

Interestingly, these pictures show that, up to the decisive point of the seizures, the continuous measures look very similar, thus hinting at an instantaneous finite-size fluctuation causing the whole network to explode. Also, it is noteworthy that the new “hierarchy” set by the cascade is conserved in the non-propagative regime, while the propagative regime seems to have an overall reset effect.

Also, we see from these graphs that there is a particular time window where the variance of the mean membrane potential is larger for the most inhibitory-connected neurons, in both RS and FS populations (although it appears

clearer for RS ones here, because of the need to rescale the FS colorbar to have comparable results). This increase of variance, while still present, is weaker and on a smaller time window in the average over connectivities compared with the average over noise realizations. This suggests that connectivity plays a role in the intensity of the effect, although it remains qualitatively similar. We found that this time window defines the period when the network can actually switch to propagation: the high variance corresponds to different times when various realizations “ explode,” and thus defines a region of instability.

A central point to raise here is that what makes the difference between propagative and non-propagative scenarios is most likely not an infinitesimal instability defined
<figure><img src="image_18.png" /><figcaption>**Figure 11.** Kuramoto *R* of membrane potentials over subgroups of neurons (same network connectivity, different noise realizations) for each group defined as a function of their incoming inhibitory connections, averaged over 50 noise realizations (17 non-propagative and 33 propagative). *a*, *b*, Color maps correspond for each group to the average Kuramoto parameter (top) and SD (bottom) across noise realizations in the propagative (*a*) and non-propagative (*b*) scenarios for both excitatory (RS) and inhibitory (FS) populations. The blue rectangle highlights the (time) region where the system either switches to a propagative regime or remains stable.</figcaption></figure>

from a macroscopic perspective [i.e., that is because of a positive eigenvalue of a Jacobian defined from a large-scale representation (e.g., mean-field)], otherwise the non-propagative behavior would simply not be observable (as, except for chaotic dynamics, we do not observe unstable trajectories in phase space). Indeed, what differs between the various simulations is either the noise, or the connectivity realization, which may, or may not, bring the system to a point of instability. The external Poissonian drive, with finite-size fluctuations is thus constitutive of the scenarios we observe.

To gain more insight into the diversity of dynamics across neuron groups, we turn our attention to a measure of alignment, or synchronization, namely the Kuramoto order parameter *R*.

### Kuramoto order parameter

The Kuramoto parameter represents a degree of alignment, with a value of 0 meaning there is no alignment and a value of 1 meaning everyone is perfectly aligned. We show in Figure 11, *a* and *b*, that the Kuramoto order parameter *R* defined for each group of excitatory (RS) and inhibitory (FS) neurons in time, averaged over noise realizations (Fig. 11a, *b*, top row), and SD over realizations (Fig. 11a, *b*, bottom row), in propagative (Fig. 11a) and non-propagative (Fig. 11b) scenarios (network connectivity held fixed).

Again, we reproduce the results with 50 network connectivities, for both propagative and non-propagative scenarios (Fig. 12a, *b*). We clearly see that the results are qualitatively similar, although with seemingly higher contrast than Figure 11.
<figure><img src="image_19.png" /><figcaption>b Kuramoto R Non-propagation</figcaption></figure>

<figure><img src="image_20.png" /></figure>

**Figure 12.** Kuramoto *R* of membrane potentials over subgroups of neurons (different connectivities) for each group defined as a function of their incoming inhibitory connections. Here, we averaged over 50 network connectivities for which we found a couple of noise realizations corresponding to propagative and non-propagative scenarios. *a*, *b*, Color maps correspond for each group to the average membrane potential (top) and SD (bottom) across different connectivities in the propagative (a) and non-propagative scenarios (b) for both excitatory (RS) and inhibitory (FS) populations. The blue rectangle highlights the (time) region where the system either switches to a propagative regime or remains stable.

The cascade previously observed is clearly visible for the average *R*, in the form of a “desynchronization cascade.” For the propagative scenario, we note here a recruitment process between two radically different regimes having nonetheless alignment features: a fluctuation-driven AI dynamics, where membrane potentials are mostly conditioned by the balance of inhibitory versus excitatory inputs, and a seizure characterized by high spiking and membrane potentials clamped by refractoriness. Interestingly, in the non-propagative scenario, it appears that the misalignment of the inhibitory neuron groups finally attained is fueled by the joint activity (of the network and the input), thus hinting at an

out-of-equilibrium steady state (that continues until the end of the plateau of the perturbation, 1 s later). From the SD perspective, two main features are worth pointing out. First, we again observe the instability window, characterized by high SD between realizations in propagative scenarios. Second, we see that the two types of averaging leads to strikingly similar results, although slightly different quantitatively speaking, the average over connectivities leading to a higher contrast during the cascade. Therefore, our results are independent of both the noise realization and the specific connectivity, although an average over one or the other is useful to observe a typical case.
<figure><img src="image_21.png" /><figcaption>**Figure 13.** a, b, Steady-state and dynamic profiles of RS neurons for $\mu_V$ (a) and Kuramoto R (b) over subgroups of neurons [same network connectivity (unless specified), different noise realizations], for fixed external input. The steady states (Stat.) represent the stable activity without perturbation. They are drawn together with various profiles for different amplitudes of perturbation (Dyn) captured right before seizure onset, at, respectively, 1950 ms (60 Hz), 1950 ms (80 Hz), and 1930 ms (100 Hz, as the seizures develop before 1950 ms). Networks are the same as previously analyzed, except, when stated, Net 2, which represents another network connectivity, for robustness. SEs estimated over noise realizations are shown in shaded areas.</figcaption></figure>

## Dynamic versus static approach

We have seen that changing the slope and the amplitude of the signal alters the chances of triggering a seizure, thus hinting that the time evolution of the perturbation is central. Then we observed a hierarchical structure setting in from the point of view of continuous measures, following the perturbation. However, fundamental questions remain. How much of this latter phenomenon is actually dynamic? Would we find the same structures if we bombarded the network with a fixed input at, say, 80 Hz? Can we observe the same dynamic structures for scenarios that are always, or never, propagative (no matter the noise or connectivity realization)? This would indicate that the structures observed thus far might have little to do with the seizure phenomenology itself but would either be the mere results of strong conditioning of the network by the level of input (if static structures are similar), or simply not yield any explanation for the instability we observe (if always/ never propagative scenarios show similar features).

We now turn our attention to Figure 13a, which displays the static $\mu_V$ profiles in RS population obtained for fixed external inputs ("Stat" curves), together with the profiles captured at the typical onset of the seizure, for various amplitudes: 60 Hz (never propagative), 80 Hz (sometimes propagative), and 100 Hz (always propagative). The network realization is the same as previously analyzed, except when explicitly stated (Net 2), where we refer to another connectivity realization. For the 80 Hz scenarios with the first network (the one we have been investigating so far), we kept the splitting of the realizations between propagative and non-propagative to highlight the potential differences of structures.

First, as previously observed, the profiles obtained for propagative versus non-propagative regimes are very similar for lower values of inhibitory connectivity. Then, we clearly see that the $\mu_V$ profiles extracted from the dynamic situations (hereafter called the "dynamic profiles") are very different from the static ones.

In addition, it is worth pointing out that the profile obtained for an 80 Hz amplitude with a different random realization of the network (where all 50 noise realizations are put together, based on the previous observation that propagative and non-propagative scenarios show very similar structures) is very similar to those already shown, with small SE, which, together with the previous observation that noise and network realizations seem to play similar roles, underlie a robust network phenomenology. Furthermore, we see that the profiles obtained for 60, 80, and 100 Hz amplitudes are different. The nature of their differences is of great interest for low indices, where we observe that 60 and 100 Hz profiles are located on opposite sides of the central 80 Hz profile: their ordering in this region is consistent with that of their response to the perturbation we have observed so far (Fig. 4). This said, the dynamic profiles yet show similar qualitative features: they all are nonmonotonous and display two well separated parts. Indeed, for low indices (until 30) $\mu_V$ is increasing with values starting around the lowest of the static profiles (10 Hz), while their high-indices part is more aligned with high static profiles. Interestingly, we see that for 60 and 80 Hz the right part is well aligned with the static profile obtained for similar inputs. This does not seem to be the case for 100 Hz, although the static input simulation displays some instability, which makes their comparison less relevant. Although it is not straightforward to link $\mu_V$ with the instantaneous regime, we have seen that low values can be associated with high firings (the neurons spending most of their time clamped at -65 mV). This helps in understanding what is happening here: for higher values of amplitude, the less inhibitory-connected neurons are firing more, and can thus lead the rest of the network to higher activities.

Figure 13b shows the Kuramoto order parameter aspect of the latter figure. Here the R profiles display structures quite different from those observed for $\mu_V$. Indeed,
<table><caption>Table 1: Triggered and prevented events</caption><thead><tr><th></th><th>Time of peak</th><th>+1 Hz</th><th>+5 Hz</th><th>-1 Hz</th><th>-5 Hz</th></tr></thead><tbody><tr><td rowspan="6">Percentage of prevented events</td><td>t = 1500 ms</td><td>0.1806</td><td>0.1944</td><td>0.1528</td><td>0.1389</td></tr><tr><td>t = 1850 ms</td><td>0.1389</td><td>0.1944</td><td>0.0972</td><td>0.1528</td></tr><tr><td>t = 1950 ms</td><td>0.1528</td><td>0.2361</td><td>0.125</td><td>0.0694</td></tr><tr><td>t = 1975 ms</td><td>0.0972</td><td>0.0</td><td>0.3472*</td><td>0.3889*</td></tr><tr><td>t = 2000 ms</td><td>0.0139</td><td>0.0</td><td>0.25*</td><td>0.5556**</td></tr><tr><td>t = 2500 ms</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0972</td></tr><tr><td rowspan="6">Percentage of triggered propagation events</td><td>t = 1500 ms</td><td>0.25*</td><td>0.1786</td><td>0.2857*</td><td>0.25*</td></tr><tr><td>t = 1850 ms</td><td>0.178</td><td>0.1786</td><td>0.2143</td><td>0.2143</td></tr><tr><td>t = 1950 ms</td><td>0.0357</td><td>0.6071**</td><td>0.2143</td><td>0.5**</td></tr><tr><td>t = 1975 ms</td><td>0.7143**</td><td>1.0**</td><td>0.25*</td><td>0.28572*</td></tr><tr><td>t = 200 ms</td><td>0.6071**</td><td>1.0**</td><td>0.0</td><td>0.0714</td></tr><tr><td>t = 2500 ms</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td></tr></tbody></table>

The percentage of prevented events refers to 72 initially propagative behaviors. *, ≥25%; **, ≥50%. The time of peak corresponds to the moment where the maximum of the stimulus is reached, and the amplitude corresponds to a variation of the external input (see the main text). The percentage of triggered propagation events refers to an initial number of 38 non-propagative cases. *, ≥25%; **, ≥50%.

the various static profiles do not display such clear variability as for $\mu_v$, although little differences can still be observed: high inputs seem to show more variability in low indices, while ending at higher values for higher indices. More importantly, the dynamic profiles are here very different from the static ones and among themselves. In addition, the simulated propagative and non-propagative scenarios show little differences here as well, and the profiles corresponding to same amplitude (80 Hz) and different network architecture (Net 2) also overlap here. Interestingly, we can also observe that the 60 and 100 Hz profiles are different and located apart from the 80 Hz, although they also show different magnitudes of their inverted peaks. Given that the ordering of these magnitudes is not consistent with the various degrees of instability, we suggest that the position of the peak might be the most discriminating factor to establish whether the scenario is propagative. This would be consistent with the observations we made thus far and would confirm our previously suggested scenario: the more we approach the centre group, the more neurons are considered (binomial distribution). In Fig. 13, the green peak (100Hz) tells us that more neurons in that region have undergone the desynchronization cascade we mentioned earlier; that is, more neurons have already “switched sides” and entered a high-firing regime, thus giving more inertia to the cascade phenomenon. The middle scenario (80 Hz) would then sit on a tipping point, that is a point separating two radically different dynamic regimes of the system.

These latter observations show that, from the perspective of both mean membrane potential and Kuramoto order parameter calculated inside the groups formed from inhibitory in-degree, we are in the presence of a structured behavior that emerges from an intricate interaction between dynamics and architecture, and which cannot be recovered from static approaches.

### Can seizure propagation be controlled by external inputs?

After having established that the structure of the dynamics allows or does not allow the propagation of the paroxysmal perturbation, we now investigate whether we

could use the previous finding of a strong instability window for the 80 Hz dynamic scenario to alter the fate of the AdEx network dynamics. This approach is based on the following reasoning: we have observed, with a detailed analysis, that switching to one scenario or another is determined in a short time window (just before the eventual seizure). Thus, we want to design a stimulation protocol to reduce the chance of seizure propagation, based on this observation, but one that does not require the same level of analysis, hence making it applicable inline and without the need of extensive computational power. To do so, we will study the region around the seizure to determine this relevant time window.

To achieve this, we apply a Gaussian stimulation, with a 10 ms time constant, two different amplitudes (1 and 5 Hz), positive or negative, through a variation of the external excitatory input (which, depending on when the simulation is applied, can be the drive of 6 Hz or the drive plus somewhere on the perturbation of 80 Hz with a time constant of 100 ms). For simulations performed under the same conditions, the stimulations were applied at different times, as detailed in Table 1. Table 1 shows, for a total of 100 simulations (with same network connectivity but different noise realizations), among which 72 were propagative, that the relative percentage of simulations has undergone a triggering and a cancelation of the seizure, respectively.

We see that it is possible to “reverse” the scenario from propagative to non-propagative in the time windows between 1975 and 2000 ms (and vice versa, albeit for a larger time window) thanks to (or because of) the stimulation (as can be seen in Table 1, the Percentage of prevented events section; for the opposite, see Table 1, the Percentage of triggered propagation events section). A notably interesting case is that >50% of the seizures are prevented if a stimulation of -5 Hz is applied in the same time window. This could open interesting leads in furthering qualitative comparisons between computational simulations and real-life situations, and could eventually guide future interventions.

### Discussion

In this computational work, we studied the response of various spiking neural networks to paroxysmal inputs. We
observed that the same networks can display various types of responses, depending on their nature (the neuron model used at its nodes), the shape of the perturbation (here we analyzed particularly a plateau-like input with various slopes and amplitudes), and the realization of the random number generator. In the case of AdEx and CAdEx networks, two radically different responses to a qualitatively similar in- coming excitatory perturbation are observed. Indeed, the lat- ter could either recruit the excitatory population and thus allow the seizure to propagate to efferent areas or be “con- trolled” by the activity of the inhibitory population, keeping the excitatory population at a low activity level, thus prevent- ing further propagation. The response of the network de- pends not only on the amplitude of the perturbation but also on its rising speed. This is consistent with experimental ob- servations (Saggio et al., 2020). Interestingly, in the case of an HH network, our investigations show very different network responses, where mostly the amplitude of the perturbation plays a role and where no variability on noise realizations was observed.

A rich literature shows that seizures can be classified according to their onset/offset features described by bi- furcation types (Jirsa et al., 2014; Saggio et al., 2017, 2020). The most observed bifurcation at the onset of a seizure is a saddle-node bifurcation (Saggio et al., 2017), which is characterized by an abrupt change in the base- line of the electrophysiological signal (Jirsa et al., 2014). We observed in the current work that seizures are propa- gative in AdEx and CAdEx networks when they rise abruptly enough in the network. There is here an interest- ing correspondence revealing the importance of the onset of seizure dynamics, as it has been shown from a clinical point of view (Lagarde et al., 2019). It is worth noting that the absence of such phenomenology in HH networks (for the scenarios we considered) raises interesting questions in the modeling of seizure dynamics, but also more gener- ally in neuronal networks. How do the quantitative differ- ences (number of variables) and qualitative differences (types of processes taken into account) in the single-neu- ron models affect the global dynamics? Are more precise models always the best in all respects? This underlines the importance of the choice of model and of parameters: by modeling a neuronal network and observing a phenomenon that resembles reality, we are not testing whether the specif- ic ingredients we chose are constitutive of this phenomen- on, but how they would be if they were chosen a priori. It is only the systematic cross-model observations and compari- sons that can yield an answer as to which are the necessary and sufficient ingredients to observe a given phenomenon.

Note that, in clinical observations, the most accessible measurements are made on a macroscopic scale. In the study proposed here, we observe the activities at a smaller descriptive scale by building a network of neuron models. We thus have a complex system of very high dimension, rendering impossible, a priori, the obtaining of a simple de- scription of the dynamics, which motivates the statistical ap- proach proposed here. With this type of analysis, we were able to track in time key features of the underlying dynamics, especially those supported by the structure of the network: inhibitory in-degree can be mobilized to explain global

differences in network response. Indeed, we proposed a coarse-grained description of the network dynamics based on inhibitory in-degree, allowing us to capture internal proc- esses that were not visible at first, and that play a significant role in the global out-of-equilibrium dynamics. We chose in- hibitory in-degree as it was found to be the most influential aspect determining the firing rate (Fig. 5a). It is interesting to note that inhibitory neurons were also the ones that had the highest firing rate (~15 Hz), while the excitatory neurons were way lower (~2 Hz) and the Poisson noise was also lower (~6 Hz by construction). That difference could be the reason for the disparity in influence more than the nature of the neurons, and while it could be interesting to investigate, it is not in the scope of this study and does not change the main results as the categorization was only used as a tool to visualize the data. This opens the way to a flexible modeling framework of internal subpopulations, whose precision can be adapted to the most significant level of description, de- pending on the context and the questions asked. This is a first bottom-up step toward a coarser description of the sys- tem and, hence, may guide reliable modeling attempts at larger scales.

We have also established that not only does this struc- ture matter, but also its interaction with instantaneous fi- nite-size fluctuations of the noise and the time evolution of the global dynamics. These are all constitutive of the observed behaviors, and none can be neglected to under- stand them.

Also, our results showed that, for the AdEx network, there exists a time window, characterized by a high var- iance across noise realizations, during which it is possible to reverse the behavior by applying an appropriate stimu- lation. The use of a stimulus to interrupt a seizure has been applied in the past in the case of absence seizure (Rajna and Lona, 1989). These results have been used as bases of computational studies at the scale of the EEG (Taylor et al., 2014). Computational work on the response of a network model to stimuli to disrupt seizure-like activ- ities has shown the importance of the precise timing of the stimulation (Anderson et al., 2007). Then, the use of electrode stimulation has been developed in rodents (Pais-Vieira et al., 2016). These different approaches have been implemented, including deep brain stimulation, vagus nerve stimulation (Boon et al., 2018), and magnetic stimulation (Ye and Kaszuba, 2019). However, experimen- tal recordings of the response to stimuli do not allow us to understand the mechanisms of large populations of neu- rons. Indeed, even if progress in calcium imaging or in multielectrode arrays has made it possible since this last decade to record a large number of neurons simultane- ously, we do not yet have access to the exact structure of the network they constitute. The study presented here is thus a proof of concept, based on a specific network model.

Finally, we also found that it is possible to “control” the propagation of the seizure by appropriate stimulation in a given time window. We think that this constitutes not only an important prediction of the model, but also a poten- tially important possibility for the treatment of some types of intractable focal epilepsies. This prediction could be
tested in future modeling work at the mesoscopic scale, with realistic connectivity between the focus and neighboring areas. Such a model could be used to test the hypothesis that appropriate stimulation in areas adjacent to the focus may prevent the propagation of the seizure.

Perhaps the most exciting perspective is that the same paradigm could be used experimentally to control seizures. This would require a system to detect the onset of the seizure in the focus, and another system to deliver appropriate stimuli in adjacent areas. Such a system could be applied to experimental models of focal seizures, to evaluate whether such a paradigm could reverse the propagation, and thus generalization, of the seizure. This could be another way of controlling seizures, not by suppressing the focus, but by making sure that the paroxysmal activity does not propagate.

## References

Anderson W, Kudela P, Cho J, Bergey G, Franaszczuk P (2007) Studies of stimulus parameters for seizure disruption using neural network simulations. *Biol Cybern* 97:173-194.

Beghi E (2020) The epidemiology of epilepsy. *Neuroepidemiology* 54:185-191.

Boon P, Cock ED, Mertens A, Trinka E (2018) Neurostimulation for drug-resistant epilepsy. *Curr Opin Neurol* 31:198-210.

Brette R, Gerstner W (2005) Adaptive exponential integrate-and-fire model as an effective description of neuronal activity. *J Neurophysiol* 94:3637-3642.

Brunel N (2000) Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. *J Comput Neurosci* 8:183-208.

Carlu M, Chehab O, Dalla Porta L, Depannemaecker D, Héricé C, Jedynak M, Köksal Ersöz E, Muratore P, Souihel S, Capone C, Zerlaut Y, Destexhe A, di Volo M (2020) A mean-field approach to the dynamics of networks of complex neurons, from nonlinear integrate-and-fire to Hodgkin-Huxley models. *J Neurophysiol* 123:1042-1051.

Dehghani N, Peyrache A, Telenczuk B, Le Van Quyen M, Halgren E, Cash SS, Hatsopoulos NG, Destexhe A (2016) Dynamic balance of excitation and inhibition in human and monkey neocortex. *Sci Rep* 6:23176.

Depannemaecker D, Destexhe A, Jirsa V, Bernard C (2021) Modeling seizures: from single neurons to networks. *Seizure* 90:4-8.

Depannemaecker D, Ivanov A, Lillo D, Spek L, Bernard C, Jirsa V (2022) A unified physiological framework of transitions between seizures, sustained ictal activity and depolarization block at the single neuron level. *J Comput Neurosci* 50:33-49.

Destexhe A (2009) Self-sustained asynchronous irregular states and up-down states in thalamic, cortical and thalamocortical networks of nonlinear integrate-and-fire neurons. *J Comput Neurosci* 27:493-506.

El Boustani S, Pospischil M, Rudolph-Lilith M, Destexhe A (2007) Activated cortical states: experiments, analyses and models. *J Physiol Paris* 101:99-109.

Górski T, Depannemaecker D, Destexhe A (2021) Conductance-based adaptive exponential integrate-and-fire model. *Neural Comput* 33:41-66.

Hille B (1992) Ionic channels of excitable membranes. *Sunderland, MA: Sinauer Associates*.

Hodgkin AL, Huxley AF (1952) A quantitative description of membrane current and its application to conduction and excitation in nerve. *J Physiol* 117:500-544.

Jirsa VK, Stacey WC, Quilichini PP, Ivanov AI, Bernard C (2014) On the nature of seizure dynamics. *Brain* 137:2210-2230.

Jiruska P, de Curtis M, Jefferys JGR, Schevon CA, Schiff SJ, Schindler K (2013) Synchronization and desynchronization in epilepsy: controversies and hypotheses. *J Physiol* 591:787-797.

Lagarde S, Buzori S, Trebuchon A, Carron R, Scavarda D, Milh M, McGonigal A, Bartolomei F (2019) The repertoire of seizure onset patterns in human focal epilepsies: determinants and prognostic values. *Epilepsia* 60:85-95.

Nadler JV, Spencer DD (2014) What is a seizure focus? *Adv Exp Med Biol* 813:55-62.

Naud R, Marcille N, Clopath C, Gerstner W (2008) Firing patterns in the adaptive exponential integrate-and-fire model. *Biol Cybern* 99:335.

Neumann AR, Raedt R, Steenland HW, Sprengers M, Bzymek K, Navratilova Z, Mesina L, Xie J, Lapointe V, Kloosterman F, Vonck K, Boon PAJM, Soltesz I, McNaughton BL, Luczak A (2017) Involvement of fast-spiking cells in ictal sequences during spontaneous seizures in rats with chronic temporal lobe epilepsy. *Brain* 140:2355-2369.

Pais-Vieira M, Yadav AP, Moreira D, Guggenmos D, Santos A, Lebedev M, Nicolelis MAL (2016) A closed loop brain-machine interface for epilepsy control using dorsal column electrical stimulation. *Sci Rep* 6:32814.

Paulk AC, Kfir Y, Khanna AR, Mustroph ML, Trautmann EM, Soper DJ, Stavisky SD, Welkenhuysen M, Dutta B, Shenoy KV, Hochberg LR, Richardson RM, Williams ZM, Cash SS (2022) Large-scale neural recordings with single neuron resolution using neuropixels probes in human cortex. *Nat Neurosci* 25:252-263.

Peyrache A, Dehghani N, Eskandar EN, Madsen JR, Anderson WS, Donoghue JA, Hochberg LR, Halgren E, Cash SS, Destexhe A (2012) Spatiotemporal dynamics of neocortical excitation and inhibition during human sleep. *Proc Natl Acad Sci U.S.A* 109:1731-1736.

Rajna P, Lona C (1989) Sensory stimulation for inhibition of epileptic seizures. *Epilepsia* 30:168-174.

Saggio ML, Spiegler A, Bernard C, Jirsa VK (2017) Fast-slow bursters in the unfolding of a high codimension singularity and the ultra-slow transitions of classes. *J Math Neurosc* 7:7.

Saggio ML, Crisp D, Scott JM, Karoly P, Kuhlmann L, Nakatani M, Murai T, Dümpelmann M, Schulze-Bonhage A, Ikeda A, Cook M, Gliske SV, Lin J, Bernard C, Jirsa V, Stacey WC (2020) A taxonomy of seizure dynamotypes. *Elife* 9:e55632.

Soltesz I, Staley K (Eds) (2008) Computational neuroscience in epilepsy. Available at https://www.sciencedirect.com/book/9780123736499/computational-neuroscience-in-epilepsy.

Taylor PN, Wang Y, Goodfellow M, Dauwels J, Moeller F, Stephani U, Baier G (2014) A computational study of stimulus driven epileptic seizure abatement. *PLoS One* 9:e114316.

Ye H, Kaszuba S (2019) Neuromodulation with electromagnetic stimulation for seizure suppression: from electrode to magnetic coil. *IBRO Rep* 7:26-33.

Zerlaut Y, Chemla S, Chavane F, Destexhe A (2018) Modeling mesoscopic cortical dynamics using a mean-field model of conductance-based networks of adaptive exponential integrate-and-fire neurons. *J Comput Neurosci* 44:45-61.