Neural mass modeling of slow-fast dynamics of seizure initiation and abortion

Abstract

Epilepsy is a dynamic and complex neurological disease affecting about 1% of the worldwide population, among which 30% of the patients are drug-resistant. Epilepsy is characterized by recurrent episodes of paroxysmal neural discharges (the so-called seizures), which manifest themselves through a large-amplitude rhythmic activity observed in depth-EEG recordings, in particular in local field potentials (LFPs). The signature characterizing the transition to seizures involves complex oscillatory patterns, which could serve as a marker to prevent seizure initiation by triggering appropriate therapeutic neurostimulation methods. To investigate such protocols, neurophysiological lumped-parameter models at the mesoscopic scale, namely neural mass models, are powerful tools that not only mimic the LFP signals but also give insights on the neural mechanisms related to different stages of seizures. Here, we analyze the multiple time-scale dynamics of a neural mass model and explain the underlying structure of the complex oscillations observed before seizure initiation. We investigate population-specific effects of the stimulation and the dependence of stimulation parameters on synaptic timescales. In particular, we show that intermediate stimulation frequencies (>20 Hz) can abort seizures if the timescale difference is pronounced. Those results have the potential in the design of therapeutic brain stimulation protocols based on the neurophysiological properties of tissue.

Fig 1. SEEG signals recorded in two patients with epilepsy during the interictal to ictal transition and simulated signals.

(a) Epileptic seizure recorded in a patient showing the typical pre-ictal spiking pattern with three phases: sporadic spikes, pre-ictal bursts, and fast onset. (b) Zoom into each phase of the actual SEEG signal. The pre-ictal burst type-1 is followed by the pre-ictal burst type-2. (c) Simulated signals corresponding to each phase. (d) Epileptic seizure recorded in a different patient showing the typical pre-ictal spiking pattern with three phases. The region of pre-ictal bursting is zoomed in panel (e).

Fig 2. Bifurcation diagrams of the system (1) showing the different pre-ictal stages numbered 1–3 in Fig 1.

(a) Amplitude of the PSP of the pyramidal cell subpopulation is plotted as a function of the synaptic gain of the SOM+ interneuron subpopulation B. The other system parameters are given in Table 1 with G = 35. Bold and dashed lines correspond to stable and unstable solutions, respectively. Region 1 (blue) corresponds to sporadic bursts, region 2 (orange) to sustained bursts, and region 3 (purple) to low voltage fast onset activity. The system yields large amplitude ≈30 Hz oscillations in the unnumbered green shaded region. The unnumbered yellow shaded area corresponds to high amplitude stable equilibrium points, and white corresponds to physiological background activity. The arrow shows the route from background to low voltage fast onset activity in the parameter space. In order to preserve the readability of the diagram, the bifurcations along the branches of periodic solutions are not shown. (b) Co-dimension 2 diagram of the Hopf (H) and limit point (LP) bifurcations marked on panel (a) in the parameter space of B and G (synaptic gain of the PV+ interneuron subpopulation). The LP₁ and LP₂ points merge on a cusp (CP) bifurcation, and the H₁ and H₂ merge on a zero-Hopf (ZH) bifurcation.

Fig 3. Bursting orbit of system (3).

(a) Solution of (3) (blue orbit) and L⁰ (red curve) on the critical surface S⁰(green surface) projected on the (v₀, v₂, v₃)-space. Single-headed, double-headed and triple-headed arrows indicate the flow direction during superslow, slow and fast time-scales, respectively. LP denotes limit point bifurcation. The L⁰ curve changes stability at two limit points, LP₁ and LP₂ (red dots). The middle branch of the L⁰ curve between these limit points is unstable (dashed). (b) Time course of the variables (v₃, v₀, v₂) of the orbit plotted in panel (a). (c) Solution of (3) projected on the bifurcation diagram (black curve) of (4) for ε = 0 where v₂ is threaded as a parameter. Arrows show the direction of the flow with respective time-scales. Bold and dashed lines correspond to stable and unstable solutions, respectively. H donates a Hopf bifurcation, LP a limit point bifurcation. The equilibrium points along the black Z-shaped curve are unstable on the middle branch of the curve, between LP₁ at $v_2^{L{P}_1}=4.778$ and LP₂ at $v_2^{L{P}_2}=20.66$ (black dots), and on the upper branch between H₁ at $v_2^{H_1}=0.27$ and H₂ at $v_2^{H_2}=14.27$ (green dots). The amplitude of the stable limit cycles is bounded by the green continuous curves connecting the H₁ and H₂ in the ε = 0 limit.

Fig 4. Spike number as a function of C₅.

(a) Solution of (4) with 3 spikes for C₅ = 500 projected on the bifurcation diagram of the fast system (5) as a function of v₂. Arrows indicate the direction of the flow. (b) Solution of (4) with 1 spike for C₅ = 300 projected on the bifurcation diagram of the fast system (5) as a function of v₂. Arrows indicate the direction of the flow. (c) Co-dimension 2 diagrams of the Hopf (H) points (green) and the limit points (LP, black) in the (v₂, C₅) parameter space marked on the left and middle panels. As C₅ decreases, H₂ moves leftwards and eventually the spike number decreases. For C₅ = 139, H₂ and LP₁ are aligned at v₂ = 4.778. A further decrease in C₅ places H₂ on the left of LP₁ and leaves no chances for a bursting solution.

Fig 5. Geometrical analyses of a constant input.

Constant input is applied to SOM+ interneurons (a), to pyramidal cells (b) and to both SOM+ interneurons and pyramidal cells (c). Left panels show the projection of the nullsurfaces, critical slow manifold and the orbit of the reduced model (7). Right panels show the LFP signal of the full system (3) subject to the constant inputs analyzed on the left. All parameters are as given in Table 1, except B = 15. (a) The y₅-nullsurface Θ (blue surface for k_P = 0), and y₇-nullsurface Σ (red surface for k_SOM = −1, black surface for k_SOM = 0, green surface for k_SOM = 1) are projected on the (v₂, y₅, v₀)-space. The blue curve L⁰ (stable on the bold, unstable on the dashed) is on the intersection between Θ and the {y₅ = 0}-hyperplane. The black and red orbits are the solutions of the system for k_SOM = 0 and k_SOM = −1, respectively. For k_SOM = 1, the solution approaches to the green stable equilibrium point on the intersection between Σ(k_SOM = 1) and L⁰. Panel (a1) shows the time series for k_SOM = {0, 1}, and panel (a2) for k_SOM = −1. (b) The y₅-nullsurface Θ (red surface for k_P = 1, green surface for k_P = −1), and y₇-nullsurface Σ (black surface for k_SOM = 0) are projected on the (v₂, y₅, v₀)-space. The red curve L⁰ (stable on the bold, unstable on the dashed) is on the intersection between Θ(k_P = 1) and the {y₅ = 0}-hyperplane. The green curve L⁰ (stable on the bold, unstable on the dashed) is on the intersection between Θ(k_P = −1) and the y₅ = 0 hyperplane. The green and red orbits are the solutions of the system for k_P = 1 and k_P = −1, respectively. Panel (b1) shows time series for k_P = 1, and panel (b2) for k_P = −1. (c) The y₅-nullsurface Θ (red surface for k_P = 1) and y₇-nullsurface Σ (green surface for k_SOM = 1, blue surface for k_SOM = 2) are projected on the (v₂, y₅, v₀)-space. The red curve L⁰ (stable on the bold, unstable on the dashed) is on the intersection between Θ(k_P = 1) and the y₅ = 0 hyperplane. The green curve L⁰ (stable on the bold, unstable on the dashed) is on the intersection between Θ(k_P = 1) and the {y₅ = 0}-hyperplane. The green orbit is the solution of the system for (k_P, k_SOM) = (1,1). For (k_P, k_SOM) = (1,2) the solution approaches to the cyan stable equilibrium point on the intersection between Σ(k_SOM = 2) and L⁰. Panel (c1) shows time series for (k_P, k_SOM) = (1,1), and panel (c2) for (k_P, k_SOM) = (1,2).

Fig 6. System (1) under stimulation.

Biphasic stimulation with a 0.5 ms pulse width (total pulse duration is 1 ms) is applied to the system in type-1 bursting. Panels (a) and (b) show the energy map of the simulated LFP signal that is lower in the blue region than the yellow region (see the color bar on the right). The energy of the LFP signal was computed by using $E={\sum }_{n=1}^{n=N}{\left|x\left(n\right)\right|}^2,$ where n stands for the index of and N for the size of the discrete signal x(n). (a) Only the SOM+ interneurons receive the biphasic perturbation. (b) The pyramidal cell, SOM+ interneurons and PV+ interneurons receive the same biphasic perturbation. Panels (c), (d) and (e) show the time course of the marked stimulation on panels (a) and (b). (c) 15 Hz biphasic stimulation with 10 arb. unit amplitude is applied to the SOM+ interneurons (k_SOM = 1, k_P = k_PV = 0). (d) 15 Hz biphasic stimulation with 10 arb. unit amplitude is homogenously applied to all subpopulations (k_SOM = k_P = k_PV = 1). (e) 25 Hz biphasic stimulation with 10 arb. unit amplitude is homogenously applied to all subpopulations (k_SOM = k_P = k_PV = 1).