# Expansion of epileptogenic networks via neuroplasticity in neural mass models

Elif Köksal Ersöz, Pascal Benquet, Fabrice Wendling

► To cite this version:

Elif Köksal Ersöz, Pascal Benquet, Fabrice Wendling. Expansion of epileptogenic networks via neuroplasticity in neural mass models. PLoS Computational Biology, 2024, 20 (12), pp.e1012666. 10.1371/journal.pcbi.1012666. hal-04833654

HAL Id: hal-04833654

https://hal.science/hal-04833654v1

Submitted on 12 Jun 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# Expansion of epileptogenic networks via neuroplasticity in neural mass models

Elif Köksal-Ersöz*, Pascal Benquet, Fabrice Wendling

Univ Rennes, INSERM, LTSI UMR 1099, Rennes, France

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

## Abstract

Neuroplasticity refers to functional and structural changes in brain regions in response to healthy and pathological activity. Activity dependent plasticity induced by epileptic activity can involve healthy brain regions into the epileptogenic network by perturbing their excitation/inhibition balance. In this article, we present a new neural mass model, which accounts for neuroplasticity, for investigating the possible mechanisms underlying the epileptogenic network expansion. Our multiple-timescale model is inspired by physiological calcium-mediated synaptic plasticity and pathological extrasynaptic N-methyl-D-aspartate (NMDA) dependent plasticity dynamics. The model highlights that synaptic plasticity at excitatory connections and structural changes in the inhibitory system can transform a healthy region into a secondary epileptic focus under recurrent seizures and interictal activity occurring in the primary focus. Our results suggest that the latent period of this transformation can provide a window of opportunity to prevent the expansion of epileptogenic networks, formation of an epileptic focus, or other comorbidities associated with epileptic activity.

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

OPEN ACCESS

**Citation:** Köksal-Ersöz E, Benquet P, Wendling F (2024) Expansion of epileptogenic networks via neuroplasticity in neural mass models. PLoS Comput Biol 20(12): e1012666. https://doi.org/10.1371/journal.pcbi.1012666

**Editor:** Lyle J. Graham, Centre National de la Recherche Scientifique, FRANCE

**Received:** February 14, 2024

**Accepted:** November 21, 2024

**Published:** December 3, 2024

**Copyright:** © 2024 Köksal-Ersöz et al. This is an open access article distributed under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** Simulation files are available at https://github.com/elifkoksal/plasticNMM.

**Funding:** This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 855109 to FW). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared that no competing interests exist.

## Author summary

Repetitive epileptic activity leads to a large variety of plastic changes in brain. Understanding the mechanism and its time course can indicate windows of opportunity for preventing the devastating impact of focal epileptic activity on brain networks. Animal models of focal epilepsy have underlined the role of NMDA receptors and deregulations in the GABAergic system in the expansion of epileptic activity. Here we develop a plastic neural mass model integrating physiological and pathologic plasticity. Our results show that neuroplasticity mediated by epileptic activity increases synaptic efficacy between epileptic and healthy regions, and causes irreversible changes in the GABAergic system in the healthy region. The model gives mechanistic insights into the physio-pathological impact of focal epileptic activity on large-scale networks and suggests that seizures may beget seizures.

## Introduction

The most general definition of neuroplasticity is the ability of neural networks to reorganize over time, which can take place under physiological and pathological conditions. Activity

<footnote>□ Current address: Inria, Lyon Neuroscience Research Center, Bron, France</footnote>

<footnote>* elif.koksal-ersoz@inserm.fr</footnote>
dependent plasticity is the ability of synapses to strengthen or weaken, in response to increases or decreases in their activity. Under physiological conditions, it is accompanied by many regulatory mechanisms through which the excitation/inhibition balance is preserved (for examples [1–6]). However, epileptiform discharges drive not only physiologic synaptic plasticity but also pathological plasticity, which disturbs the excitation/inhibition balance, and can cause expansion of epileptogenic networks during epileptogenesis [7–9]. The aim of this article is to investigate the expansion of epileptogenic networks with a mathematical model of physiopathological plasticity under recurrent epileptic activity.

Synaptic plasticity occurs on multiple timescales at both presynaptic and postsynaptic sites [10]. Short-term plasticity (STP, in milliseconds) refers to transient changes in the probability of neurotransmitter release at the presynaptic site [11]. Long-term synaptic plasticity (in hours or days), although generally considered as postsynaptic alterations, is related to both (i) the insertion/removal of alpha-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA) receptors (AMPARs) at the postsynaptic site and (ii) the change in the probability of neurotransmitter release at the presynaptic site [12].

Postsynaptic long-term potentiation (LTP) and long-term depression (LTD) depend on the calcium concentration at the postsynaptic site [13]. In excitatory (glutamatergic) synapses, an increase in postsynaptic calcium concentration through activation of synaptic N-methyl-D-aspartate (NMDA) receptors (NMDARs), metabotropic glutamate receptors, voltage-gated calcium channels, etc. is necessary and sufficient for synaptic plasticity (see [12,14] and references therein). Calcium-dependent activation of secondary messengers, such as nitric oxide, endocannabinoids, etc., triggers presynaptic long-term plasticity [15,16]. The late phase of long-term plasticity, referred to as maintenance or consolidation, occurs via new protein synthesis and gene transcription [12]. We refer to these processes of postsynaptic long-term plasticity, presynaptic long-term plasticity and consolidation as *physiological plasticity*.

Epileptic discharges that cause glutamatergic spillover have a twofold impact. First, epileptiform activity strengthens excitatory synapses by increasing the number of AMPARs [17,18], as in the case of LTP. Second, a high amount of glutamate released during epileptiform activity can activate extrasynaptic NMDARs that cause pathological changes in GABAergic signaling, through inactivation of functional potassium chloride co-transporters 2 (KCC2) and internalization of GABAergic receptors [19,20]. KCC2 controls the outflow of the chloride ions from the intracellular space to the extracellular space and maintains the GABAergic reversal potential. Calcium ions flowing through extrasynaptic NMDARs bind on calpain that inactivates KCC2 [20–22]. Reducing the number of functional KCC2 results in the increase of the intracellular chloride concentration and changes GABAergic responses from inhibitory to excitatory. KCC2 dysfunction is related to epileptic activity in the naive site [7], and it is one of the factors contributing to epileptic activity in the mature brain [23].

The extrasynaptic NMDAR activity also affects tonic inhibition that is mainly mediated by the $\alpha_5$ GABA$_A$ receptors in the hippocampus. Overexpression of GluN2B NMDA receptors, which are located mostly in extrasynaptic regions, downregulates the expression of $\alpha_5$ GABA$_A$ receptors expression, promotes their internalization, and therefore, reduces tonic GABAergic [24] (but see [25] for limbic epileptogenesis). We refer to these two mechanisms, i.e. reduced KCC2 and internalization of extrasynaptic GABAR, as *pathological plasticity*.

After several seizures physiological plasticity (LTP, LTS) seems to be occluded. However, pathological plasticity can increase the susceptibility of a brain region to epileptic discharges, expand epileptic network, and generate a secondary (mirror) focus. Secondary focus refers to epileptogenesis in a naive area induced by repetitively uncontrolled epileptic seizures from a primary seizure focus [26,27]. Independent epileptic activity in the contralateral area in humans suggests that mirror focus develops via a kindling process mediated by the primary
focus. In animal models of epilepsy, for instance, under kindling [28] or kainite application [7] to one of the hippocampi (primary focus), secondary focus appears in the contralateral hippocampus due to recurrent epileptic activity in the primary focus. Khalilov's elegant work in the immature brain [7] showed that activation of NMDARs and long-term alterations in GABAergic synapses are responsible for chronic epileptogenesis in the secondary focus. Despite the key role of plasticity in expansion of epileptogenic networks, there is a lack of computational modeling studies of neuroplasticity under epileptic activity.

Previous modeling studies of calcium-mediated postsynaptic LTP/LTD have mainly focused on cellular activity [29–34] (see [35] for a mean-field approximation), as well as presynaptic long-term plasticity [30,33,36,37]. Recently, Chindemi et al. [30] developed a detailed cellular model of the neocortex that integrates the calcium-mediated presynaptic and postsynaptic plasticity mechanisms. However, the literature related to modeling calcium-mediated plasticity with neuro-physiologically plausible neural mass model (NMMs) remains limited, and up to date up to date no computational model implements the activity dependent pathological plasticity. We identified two studies on this phenomenon in the context of magnetic stimulation [38,39]. The originality of this study is to differentiate the impact of physiological plasticity and the pathological plasticity induced by epileptic discharges on epileptogenicity of NMMs. Our computational model concentrates on the physio-pathological plasticity under epileptic activity in a network of two unidirectionally coupled NMMs [40], each representing an epileptic and a non-epileptic brain region, where the epileptic one perturbs the non-epileptic one. First, we study how synaptic connections evolve under epileptic activity. We then investigate the pathological plasticity triggered by the physiological plasticity. The results suggest that the interaction between physiological and pathological plasticity can explain the expansion of epileptogenic networks and the formation of secondary epileptic foci.

# Methods

## Mesoscopic model of physiological short- and long-term plasticity

The seminal study by Shouval et al. [29] has proposed a calcium control hypothesis for postsynaptic plasticity, in which intra-cellular calcium concentration controls the rate of change of synaptic efficacy. The proposed model in [29] was then developed to account for the synaptic consolidation in [32] by introducing a variable $\rho$ that describes the state of the synaptic efficacy. We combine both formulations as:

$$ \frac{d\rho}{dt} = \frac{1}{\tau_{\rho}} \left( -\rho(1-\rho)(\rho_{*} - \rho) + (1-\rho)\Omega_{p}([Ca]) - \rho\Omega_{d}([Ca]) \right), \quad (1) $$

where $\tau_{\rho}$ is the time constant of the synaptic variations, $\rho_{*} = 0.5$ is the boundary between the depressed ($\rho = 0$) and potentiated ($\rho = 1$) states. The functions $\Omega_{p,d}([Ca])$ describe how the synaptic efficacy is controlled by the calcium concentration:

$$ \Omega_{p,d}([Ca]) = \frac{\gamma_{p,d}}{1 + \exp(-\beta_{p,d}([Ca] - \theta_{p,d}))}, \quad (2) $$

where $\gamma_{p}>\gamma_{d}$ are potentiation and depression rates, $\beta_{p,d}$ are the slopes at potentiation and depression boundaries $\theta_{p}>\theta_{d}$ (Fig 1A). The depression rate is smaller than the potentiation rate for provoking a transition from depressed to potentiated state for a high probability [29,30,32]. We followed [29] for the value of slopes. We studied the impact of the boundaries on the long-term plasticity in the manuscript.
<figure><img src="image_4.png" /><figcaption>**Fig 1.** Synaptic efficacy and activation functions. (A) Calcium dependent synaptic efficacy function, Ω([Ca]) = −Ω_d([Ca]) + Ω_p([Ca]) (2), and (B) activation function, H(V_post) (5) for the parameter set given in Table 1.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012666.g001

The calcium concentration is controlled by the NMDAR-mediated current I<sub>NMDA</sub><sup>[Ca]</sup>:

$$ \frac{d[Ca]}{dt} = I_{NMDA}^{[Ca]} - \frac{[Ca]}{\tau_{ca}}, \qquad (3) $$

where τ<sub>ca</sub> is the decay time constant of calcium influx. The NMDAR-mediated calcium current I<sub>NMDA</sub><sup>[Ca]</sup> depend on both the presynaptic activity via the kinetics of glutamate binding and the postsynaptic membrane potential [29] expressed as,

$$ I_{NMDA}^{[Ca]} = G(F_{pre})h_{ca}H(V_{post}), \qquad (4) $$

where G(F<sub>pre</sub>) formulates the kinetics of the glutamate binding depending on the firing rate of the presynaptic population, F<sub>pre</sub>. The function H(V<sub>post</sub>) (Fig 1B) formulates the relation between the NMDA currents and the postsynaptic (population level) depolarization:

$$ H(V_{post}) = \frac{1}{1 + \exp(-\mu(V_{post} - V_{th}))}, \qquad (5) $$

where μ is the slope at V<sub>th</sub> = V<sub>post</sub>. The function H(V<sub>post</sub>) can be interpreted as the magnesium blockage at the population level. Parameter h<sub>ca</sub> is a scaling factor, which distinguishes the NMDAR-mediated postsynaptic depolarization from calcium entry.

Plasticity of excitatory synapses is modeled at both presynaptic and postsynaptic levels. Presynaptic plasticity corresponds to variations in glutamate release probability, running on short- and long-timescales. The presynaptic STP rule follows [41]:

$$ \frac{dr}{dt} = \frac{1-r}{\tau_r} - u r F_{pre} \qquad (6A) $$

$$ \frac{du}{dt} = \frac{U_s - u}{\tau_f} + U_s(1 - u)F_{pre} \qquad (6B) $$

where r is the amount of available neurotransmitters, u is the utilization factor, τ<sub>r</sub> is the
recovery timescale, $\tau_f$ is the facilitation timescale and $U_s$ is the neurotransmitter release probability. Since the interaction between glutamatergic synapses is depression dominated [42,43], we consider $\tau_r > \tau_f$.

The presynaptic long-term plasticity is the long-term variation in the neurotransmitter release probability $U_s$ controlled by the calcium-dependent retrograde signaling [30],

$$ \frac{dU_s}{dt} = \frac{1}{\tau_U} (U_s^d - U_s + \rho(U_s^p - U_s^d)). \qquad (7) $$

Parameters $U_s^d$ and $U_s^p$ represent the depressed and potentiated states of the release probability $U_s$, and $\tau_U$ is the time constant of the presynaptic long-term plasticity.

Postsynaptic long-term plasticity corresponds to an insertion/removal of glutamatergic AMPAR,

$$ \frac{d\tilde{C}_{\text{AMPA}}}{dt} = \frac{1}{\tau_{C_{\text{AMPA}}}} (C_{\text{AMPA}}^{d} - \tilde{C}_{\text{AMPA}} + \rho(C_{\text{AMPA}}^{p} - C_{\text{AMPA}}^{d})). \qquad (8) $$

Parameters $C_{\text{AMPA}}^d$ and $C_{\text{AMPA}}^p$ represent the depressed and potentiated states of the AMPAergic coupling strength $\tilde{C}_{\text{AMPA}}(t)$, and $\tau_{C_{\text{AMPA}}}$ is the time-constant of the postsynaptic long-term plasticity. For simplicity, we assume that presynaptic and postsynaptic long-term changes have the same time constants $\tau_U = \tau_{C_{\text{AMPA}}}$ . Since the expression of synaptic plasticity is slower than its induction [44], we consider $\tau_p < \tau_U = \tau_{C_{\text{AMPA}}}$ as in [30].

The whole model of pre- and postsynaptic plasticity reads as follows:

$$ \frac{dr}{dt} = \frac{1-r}{\tau_r} - u r F_{\text{pre}}, \qquad (9A) $$

$$ \frac{du}{dt} = \frac{U_s - u}{\tau_f} + U_s(1 - u)F_{\text{pre}}, \qquad (9B) $$

$$ \frac{d[Ca]}{dt} = I_{\text{NMDA}}^{[Ca]} - \frac{[Ca]}{\tau_{ca}}, \qquad (9C) $$

$$ \frac{d\rho}{dt} = \frac{1}{\tau_{\rho}} \left( -\rho(1-\rho)(\rho_{*} - \rho) + (1-\rho)\Omega_{p}([Ca]) - \rho\Omega_{d}([Ca]) \right), \qquad (9D) $$

$$ \frac{dU_s}{dt} = \frac{1}{\tau_U} (U_s^d - U_s + \rho(U_s^p - U_s^d)), \qquad (9E) $$

$$ \frac{d\tilde{C}_{\text{AMPA}}}{dt} = \frac{1}{\tau_{C_{\text{AMPA}}}} (C_{\text{AMPA}}^{d} - \tilde{C}_{\text{AMPA}} + \rho(C_{\text{AMPA}}^{p} - C_{\text{AMPA}}^{d})). \qquad (9F) $$

The parameter values of (9) are given in Table 1. Except the values of time constants as mentioned above, the parameter values are chosen in accordance with the dynamics of the interacting NMMs detailed in the following sections.

## The revisited neural mass model

We revisited a minimal NMM of the CA1 region of the hippocampus [40] for obtaining an autonomous model that can undergoes seizures for a fixed parameter set. The original model
<table><caption>Table 1. Variables and parameter values of (9) unless otherwise stated.</caption><thead><tr><th></th><th>Variables</th><th>Equation number</th></tr></thead><tbody><tr><td>r(t)</td><td>amount of available neurotransmitters</td><td>Eq 9A</td></tr><tr><td>u(t)</td><td>utilization factor of neurotransmitters</td><td>Eq 9B</td></tr><tr><td>[Ca](t)</td><td>calcium concentration</td><td>Eq 9C</td></tr><tr><td>ρ(t)</td><td>state of the synaptic efficacy</td><td>Eq 9D</td></tr><tr><td>C&#773;<sub>AMPA</sub>(t)</td><td>AMPAergic coupling strength</td><td>Eq 9E</td></tr><tr><td>U<sub>s</sub>(t)</td><td>neurotransmitter release probability</td><td>Eq 9F</td></tr><tr><td>θ<sub>d</sub></td><td>depression threshold</td><td>0.1</td></tr><tr><td>θ<sub>p</sub></td><td>potentiation threshold</td><td>0.4</td></tr><tr><td>γ<sub>d</sub></td><td>depression amplitude</td><td>1</td></tr><tr><td>γ<sub>p</sub></td><td>potentiation amplitude</td><td>5</td></tr><tr><td>β<sub>d</sub></td><td>depression slope</td><td>80</td></tr><tr><td>β<sub>p</sub></td><td>potentiation slope</td><td>80</td></tr><tr><td>μ</td><td>slope of NMDA-dependent postsynaptic depolarization</td><td>1</td></tr><tr><td>V<sub>th</sub></td><td>threshold of NMDA-dependent postsynaptic depolarization</td><td>5 V</td></tr><tr><td>h<sub>ca</sub></td><td>calcium current scaling factor</td><td>10</td></tr><tr><td>τ<sub>d</sub></td><td>presynaptic short-term depression time constant</td><td>0.200 sec</td></tr><tr><td>τ<sub>f</sub></td><td>presynaptic short-term facilitation time constant</td><td>0.05 sec</td></tr><tr><td>τ<sub>ca</sub></td><td>the calcium influx decay time constant</td><td>0.05 sec</td></tr><tr><td>τ<sub>ρ</sub></td><td>synaptic efficacy time constant</td><td>50 sec</td></tr><tr><td>τ<sub>U</sub></td><td>presynaptic long-term plasticity time constant</td><td>100 sec</td></tr><tr><td>τ<sub>CA</sub></td><td>postsynaptic long-term plasticity time constant</td><td>100 sec</td></tr><tr><td>U<sub>s</sub><sup>d</sup></td><td>Value of depressed release probability</td><td>0.4</td></tr><tr><td>U<sub>s</sub><sup>p</sup></td><td>Value of potentiated release probability</td><td>0.8</td></tr><tr><td>C<sub>AMPA</sub><sup>d</sup></td><td>Value of depressed AMPAergic coupling strength</td><td>50</td></tr><tr><td>C<sub>AMPA</sub><sup>p</sup></td><td>Value of potentiated AMPAergic coupling strength</td><td>100</td></tr></tbody></table>

https://doi.org/10.1371/journal.pcbi.1012666.t001

includes four interacting neuronal subpopulations: two interconnected subpopulations of glutamatergic pyramidal neurons (P, P'), and two subpopulations of GABAergic inhibitory interneurons (somatostatin positive (SOM), and parvalbumin positive (PV), also called dendrite-projecting slow and soma-projecting fast interneurons, respectively). The activity of each subpopulation is described by a “wave to pulse” function, $S(v) = 5/(1+\exp(0.56(6-v)))$, transforming the net membrane polarization in response to synaptic inputs into a firing rate. The synaptic interactions between the subpopulations is described by a “pulse to wave” that converts the input average firing rate into a postsynaptic potential (PSP), which can be either excitatory-EPSP—or inhibitory—IPSP) at the input of each subpopulation, that is $h(t) = W/\tau_w t \exp(-t/\tau_w)$, where W represents the average PSP amplitude and $\tau_w$ is the time constant. This linear filter, also known as the alpha-function, introduces a second order ordinary differential equation $\frac{d^2v}{dt^2} = \frac{W}{\tau_w}S(v) - \frac{2}{\tau_w}\frac{dy}{dt} - \frac{1}{\tau_w^2}y$ where $y(t)$ is the PSP in response to the input $S(v)$. The equations of an unconnected population i NMMi reads:

$$ \frac{d^2 y_P^{(i)}}{dt^2} = \frac{A^{(i)}}{\tau_a} S(V_P^{(i)}) - \frac{2}{\tau_a} \frac{dy_P^{(i)}}{dt} - \frac{1}{\tau_a^2} y_P^{(i)}, \quad (10A) $$
$$ \frac{d^2 y_{p'}^{(i)}}{dt^2} = \frac{A^{(i)}}{\tau_a} (p^{(i)}(t) + C_{p',p}^{(i)} S(C_{p,p'}^{(i)} y_{p'}^{(i)})) - \frac{2}{\tau_a} \frac{dy_{p'}^{(i)}}{dt} - \frac{1}{\tau_a^2} y_{p'}^{(i)}, \quad (10B) $$

$$ \frac{d^2 y_{SOM}^{(i)}}{dt^2} = \frac{B^{(i)}}{\tau_b} S(C_{P,SOM}^{(i)} y_{P}^{(i)}) - \frac{2}{\tau_b} \frac{dy_{SOM}^{(i)}}{dt} - \frac{1}{\tau_b^2} y_{SOM}^{(i)}, \quad (10C) $$

$$ \frac{d^2 y_{PV}^{(i)}}{dt^2} = \frac{G^{(i)}}{\tau_g} S(C_{P,PV}^{(i)} y_{P}^{(i)} - C_{SOM,PV}^{(i)} y_{SOM}^{(i)}) - \frac{2}{\tau_g} \frac{dy_{PV}^{(i)}}{dt} - \frac{1}{\tau_g^2} y_{PV}^{(i)}, \quad (10D) $$

with variables $y_k^{(i)}$ representing the PSP generated by a subpopulation $k$. Parameters $C_{k,p}^{(i)}$ scales the impact from subpopulation $k$ to subpopulation $p$ within the population NMM$_i$. The net polarization of P subpopulation of NMM$_i$, $V_P^{(i)}$, reads,

$$ V_P^{(1)} = \Sigma(C_{k,P}^{(1)} y_k^{(1)}) $$

for NMM$_1$, and,

$$ V_P^{(2)} = \Sigma(C_{k,P}^{(2)} y_k^{(2)}) + \Sigma \text{Input}_{NMM_1 \to NMM_2} $$

for NMM$_2$, which is subject to the excitatory AMRAergic and NMDAergic signals from NMM$_1$ on its P subpopulation, $\Sigma \text{Input}_{NMM_1 \to NMM_2}$. We consider $V_P^{(i)}$ as the model output, which is a simplified proxy for the local field potential recorded by an extracelluar electrode positioned in the neuronal population. The model diagram is given in Fig 2A.

Intracranial recordings from epileptic patients with focal epilepsy show four distinct phases of activity: interictal phase with sporadic spikes, preictal phase with periodic preictal spikes, tonic low-voltage fast onset marked by a gamma-band activity, and clonic discharges [45]. Recordings show sequential transitions between these phases, that is, interictal, preictal, then, ictal phase (so-called seizure, which starts with a tonic low-voltage fast onset phase followed by clonic discharges). For an appropriate choice of model parameters, the transitions between these phases in (10) can be obtained by varying the IPSP amplitude of the SOM interneurons, $B^{(i)}$, (e.g. [40,46]). Fig 2B shows the bifurcation diagram of (10) for NMM$_1$ as a function of $B^{(1)}$ with the parameter set in Table 2. The branch of equilibrium points undergoes four Hopf bifurcations at $B^{(1)} \approx \{0.47,2.66,9.98,32.14\}$ and saddle-node bifurcations at $B^{(1)} \approx \{2.60,2.92,32.01,50.38\}$. The gamma-band oscillations lie in $B^{(1)} \approx (2.66,0.47)$ and the tonic lower frequency oscillations in $B^{(1)} \approx (9.98,32.14)$. Under a stochastic input, such as $p^{(i)}(t) = N(p_m^{(i)}, p_s)$ in (10B), the gamma-band oscillations can be observed for $0 < B^{(1)} < 4$, the tonic phase for $4 < B^{(1)} < 32$. The low-frequency periodic preictal spikes and aperiodic interictal spikes can be located at the boundary between tonic oscillations and stable equilibrium points for $32 < B^{(1)} < 50$.

One can ensure an autonomous transition between the phases of epileptic activity by modeling $B^{(i)}$ as a variable that varies in the range where these phases are observed. Here, we introduce a slow-subsystem where the IPSP amplitude of the SOM interneurons $B^{(i)}$ is a
<figure><img src="image_5.png" /><figcaption>**Fig 2.** Unidirectionally interacting NMM₁ and NMM₂ and the dynamics of NMM₁ in the interictal and ictal regimes. (A) Model diagram of unidirectionally coupled NMM₁ and NMM₂. Each NMM includes four interacting subpopulations of pyramidal neurons (P and P') and inhibitory interneurons (PV and SOM) with excitatory (red connections) and inhibitory (blue connections) interactions between parametrized by the coupling coefficients $C_{k,p}^{(i)}$ for $i \in \{1,2\}$ with $k$ and $p$ representing the pre- and postsynaptic subpopulations of NMMᵢ, respectively. The PSP generated by a subpopulation $k$ are denoted by $y_k^{(i)}$. NMM₂ is subject to an excitatory input from NMM₁, denoted by $Input_{NMM_1 \to NMM_2}$. (B-E) Dynamics of $Input_{NMM_1 \to NMM_2}$ for $i \in \{1,2\}$ in the interictal regime for $b_{thr}^{(1)} = 34$. (B) Bifurcation</figcaption></figure>
diagram of (10) where the amplitude of $y_p^{(1)}$ is presented as a function of $B^{(1)}$. The blue curve shows the branch of equilibrium points (bold for stable and dashed for unstable equilibrium points). The red curves show the amplitude of $y_p^{(1)}$ in the oscillatory regime. The Hopf bifurcations along the branch of equilibrium points are denoted by red dots and saddle-node bifurcation by blue dots. The range of $B^{(1)}$ values that corresponds to the fast onset, ictal and interictal periods are marked by purple, yellow and cyan patches, respectively. The time solution (black curve) is superimposed on the bifurcation diagram. (C) Time trace for $V_p^{(1)}$ showing interictal spikes. (D) Phase plane of (11) with the $B^{(1)}$-nullcline (blue curve, bold for stable and dashed for branches) and the $n^{(1)}$-nullcline (orange curve). The range of $B^{(1)}$ values that corresponds to the fast onset, ictal and interictal periods in (10) are marked by green, yellow and cyan patches, respectively. The time solution (black curve) is superimposed on the bifurcation diagram. (E) Time trace for $B^{(1)}$. (F-I) Dynamics of (10)-(11) for NMM$_i$ in the ictal regime for $b_{thr}^{(1)} = 32$. Same color codes and markers in (B-E) is used. Arrows in panel (G) indicates the direction of the trajectory (black curve) on the $(B^{(1)},n^{(1)})$-plane with double arrows indicating the fast transitions. The fast onset, ictal and interictal periods are marked on the time traces of $V_p^{(1)}$ and $B^{(1)}$ on panels H and I, respectively.

https://doi.org/10.1371/journal.pcbi.1012666.g002

variable:

$$ \frac{dB^{(i)}}{dt} = \delta \left( n^{(i)} - \left( - \frac{m_1 (B^{(i)} - p_1)^2}{1 + \exp(-p_1 + B^{(i)})} + \frac{1}{1 + \exp(-p_2 + B^{(i)})} + \frac{m_3 (B^{(i)} - p_3)^2}{1 + \exp(-p_3 + B^{(i)})} \right) \right), \quad (11A) $$

$$ \frac{dn^{(i)}}{dt} = \varepsilon \left( -n^{(i)} + n_k + \frac{n_p}{1 + \exp(-n_r(b_{thr}^{(i)} - B^{(i)}))} \right). \quad (11B) $$

Parameters of the $B^{(i)}$- and $n^{(i)}$-nullclines can be adjusted with respect to the critical (bifurcation) points of (10), desired type and range of behavior. Among these parameters, $p_1$ and $p_3$ correspond to the jump points in the $B^{(i)}$-space from ictal to interictal and from interictal to ictal, respectively. Parameters $m_1$ and $m_3$ control the stiffness of the left and right branches of the $B^{(i)}$-nullcline. Parameter $m_1$ is crucial for covering the different phases of the ictal regime starting with gamma-band oscillations and continuing by tonic spiking. Parameter $m_3$ contributes to the duration of the silent phase following the jump from ictal to interictal phase. Parameters $m_1, p_1, \delta$ and $\epsilon$ together determine the duration of the ictal phase. The parameter $b_{thr}^{(i)}$ is the main parameter controlling the excitability, hence, the susceptibility of seizing of NMM$_i$. Parameter values of (11) are given in Table 2.

Note that the system (11) is a slow-fast system itself, where $B^{(i)}$ is the fast and $n^{(i)}$ is the slow variable. The $B^{(i)}$-nullcline of (11) has an N-shape with stable outer and unstable middle branches (Fig 2D). System (10) governs the interictal and preictal phases for the $B^{(i)}$ values on the right branch $B^{(i)}$-nullcline (Fig 2B-2E). The $B^{(i)}$ values on the left and middle branches $B^{(i)}$-nullcline correspond to the ictal phase (fast-onset and tonic discharges) of (10). Depending on the position of the equilibrium point of (11), defined as the intersection between the N-shaped $B^{(i)}$-nullcline and sigmoidal $n^{(i)}$-nullcline, the system (11) follows either a steady state (Fig 2D) or an oscillatory solution (Fig 2G). If the stable equilibrium point is close to the right fold point of the $B^{(i)}$-nullcline, then (11) is in an excitable state, i.e any perturbation on $B^{(i)}$ can trigger a jump to the left branch of the fast nullcline, which corresponds to the ictal regime starting with a fast-onset activity. As the trajectory of (11) moves along the left branch of the $B^{(i)}$-nullcline, the system (10) exhibits the fast-onset, then tonic spiking before jumping back to the right branch, that is the interictal regime (Fig 2F-2I).

Table 2 shows three main differences between the parameters of uncoupled NMM$_1$ and NMM$_2$. The first difference concerns the mean values of the external inputs ($p_m^{(i)}$), which determines the range of the IPSP amplitude of the SOM interneurons ($B^{(i)}(t)$) for which sporadic and tonic spiking are observed in the model. The second difference is the IPSP amplitude of the PV interneurons ($G^{(i)}$), where higher values give gamma-band activity in the model. The
<table><caption>Table 2. Variables and parameters values for (10)-(15) unless otherwise stated.</caption><thead><tr><th colspan="2">Variables of (10)</th><th>Equation number</th></tr></thead><tbody><tr><td>y<sup>(i)</sup><sub>P</sub>(t)</td><td>EPSP emitted by the subpopulation P</td><td>Eq 10A</td></tr><tr><td>y<sup>(i)</sup><sub>P'</sub>(t)</td><td>EPSP emitted by the subpopulation P'</td><td>Eq 10B</td></tr><tr><td>y<sup>(i)</sup><sub>SOM</sub>(t)</td><td>IPSP emitted by the subpopulation SOM</td><td>Eq 10C</td></tr><tr><td>y<sup>(i)</sup><sub>PV</sub>(t)</td><td>IPSP emitted by the subpopulation PV</td><td>Eq 10D</td></tr><tr><td></td><td>Parameters of (10)</td><td>NMM<sub>1</sub> NMM<sub>2</sub></td></tr><tr><td>A<sup>(i)</sup></td><td>EPSP amplitude [mV]</td><td>5 5</td></tr><tr><td>G<sup>(i)</sup></td><td>IPSP<sub>PV</sub> amplitude [mV]</td><td>20 2</td></tr><tr><td>τ<sub>a</sub></td><td>EPSP time constant [s]</td><td>0.01 0.01</td></tr><tr><td>τ<sub>b</sub></td><td>IPSP<sub>SOM</sub> time constant [s]</td><td>0.03 0.03</td></tr><tr><td>τ<sub>g</sub></td><td>IPSP<sub>PV</sub> time constant [s]</td><td>0.003 0.003</td></tr><tr><td>C<sup>(i)</sup><sub>P,P'</sub></td><td>P → P' coupling coefficient</td><td>135 135</td></tr><tr><td>C<sup>(i)</sup><sub>P',P</sub></td><td>P' → P coupling coefficient</td><td>108 108</td></tr><tr><td>C<sup>(i)</sup><sub>SOM,P</sub></td><td>SOM → P coupling coefficient</td><td>35 35</td></tr><tr><td>C<sup>(i)</sup><sub>PV,P</sub></td><td>PV → P coupling coefficient</td><td>200 200</td></tr><tr><td>C<sup>(i)</sup><sub>P,SOM</sub></td><td>P → SOM coupling coefficient</td><td>25 25</td></tr><tr><td>C<sup>(i)</sup><sub>P,PV</sub></td><td>P → PV coupling coefficient</td><td>200 200</td></tr><tr><td>C<sup>(i)</sup><sub>SOM,PV</sub></td><td>SOM → PV coupling coefficient</td><td>120 120</td></tr><tr><td>p<sup>(i)</sup><sub>m</sub></td><td>mean of the Gaussian input p<sup>(i)</sup>(t)</td><td>90 70</td></tr><tr><td>p<sup>(i)</sup><sub>s</sub></td><td>standard deviation of the Gaussian input p<sup>(i)</sup>(t)</td><td>2 2</td></tr><tr><td></td><td>Variables of (11)</td><td>Equation number</td></tr><tr><td>B<sup>(i)</sup>(t)</td><td>IPSP<sub>SOM</sub> amplitude [mV]</td><td>Eq 11A</td></tr><tr><td>n<sup>(i)</sup>(t)</td><td>Auxiliary variable</td><td>Eq 11B</td></tr><tr><td></td><td>Parameters of (11)</td><td>NMM<sub>1</sub> NMM<sub>2</sub></td></tr><tr><td>p<sub>1</sub></td><td>local maximum of the B<sup>(i)</sup>-nullcline</td><td>25 25</td></tr><tr><td>p<sub>2</sub></td><td>inflection point of the middle branch of the B<sup>(i)</sup>-nullcline</td><td>30 30</td></tr><tr><td>p<sub>3</sub></td><td>local minimum of the B<sup>(i)</sup>-nullcline</td><td>33 33</td></tr><tr><td>m<sub>1</sub></td><td>slope of the of the left branch of B<sup>(i)</sup>-nullcline</td><td>0.0015 0.0015</td></tr><tr><td>m<sub>3</sub></td><td>slope of the of the right branch of B<sup>(i)</sup>-nullcline</td><td>0.003 0.003</td></tr><tr><td>n<sub>k</sub></td><td>minimum of the n-nullcline</td><td>-0.2 -0.2</td></tr><tr><td>n<sub>p</sub></td><td>maximum of the n-nullcline sigmoid</td><td>1.4 1.4</td></tr><tr><td>n<sub>r</sub></td><td>slope of the n-nullcline at the inflection point</td><td>2 2</td></tr><tr><td>b<sup>(i)</sup><sub>thr</sub></td><td>inflection point of the n-nullcline</td><td>34 44</td></tr><tr><td>δ</td><td>timescale parameter</td><td>50 50</td></tr><tr><td>ε</td><td>timescale parameter</td><td>0.05 0.05</td></tr><tr><td></td><td>Variables of (12)-(15)</td><td>Equation number</td></tr><tr><td>y<sub>AMPA</sub>(t)</td><td>AMPAergic-EPSP emitted by the input from NMM<sub>1</sub> onto NMM<sub>2</sub></td><td>Eq 12</td></tr><tr><td>y<sub>NMDA</sub>(t)</td><td>NMDAergic-EPSP emitted by the input from NMM<sub>1</sub> onto NMM<sub>2</sub></td><td>Eq 13</td></tr><tr><td>y<sub>NMDA,ext</sub>(t)</td><td>NMDAergic-EPSP emitted by the activation of extrasynaptic NMDAR</td><td>Eq 14</td></tr><tr><td>K(t)</td><td>Auxiliary variable representing the changes in GABAergic activity</td><td>Eq 15</td></tr><tr><td></td><td>Parameter values of (12)-(15)</td><td>NMM<sub>1</sub> NMM<sub>2</sub></td></tr><tr><td>A<sup>(2)</sup><sub>AMPA</sub></td><td>AMPAergic-EPSP amplitude [mV]</td><td>- 10</td></tr><tr><td>A<sup>(2)</sup><sub>NMDA</sub></td><td>NMDAergic-EPSP amplitude [mV]</td><td>- 2</td></tr><tr><td>A<sup>(2)</sup><sub>NMDA,ext</sub></td><td>Extrasynaptic NMDAergic-EPSP amplitude [mV]</td><td>- 1</td></tr><tr><td>τ<sub>AMPA</sub></td><td>AMPAergic-EPSP time constant [s]</td><td>- 0.005</td></tr><tr><td>τ<sub>NMDA</sub></td><td>NMDAergic-EPSP time constant [s]</td><td>- 0.02</td></tr></tbody></table>

(Continued)
<table><caption>**Table 2.** (Continued)</caption><thead><tr><th></th><th>Variables of (10)</th><th colspan="2">Equation number</th></tr></thead><tbody><tr><td>τ<sub>NMDA,ext</sub></td><td>Extrasynaptic NMDAergic-EPSP time constant [s]</td><td>-</td><td>0.04</td></tr><tr><td>C<sub>NMDA</sub></td><td>NMDAergic coupling coefficient</td><td>-</td><td>C<sub>AMPA,min</sub></td></tr><tr><td>τ<sub>K</sub></td><td>Extrasynaptic NMDA integration time constant</td><td>-</td><td>10</td></tr></tbody></table>

https://doi.org/10.1371/journal.pcbi.1012666.t002

third difference is the excitability ($b_{thr}^{(i)}$), which controls both $B^{(i)}(t)$ and the susceptibility of an NMM to jump from the interictal phase to seizure. The lower $b_{thr}^{(i)}$, the more susceptible the NMM is to exhibiting a seizure. For simplicity, we considered the same values for the rest of the parameters. The bifurcation diagram of uncoupled NMM$_2$ and its response to the epileptic activity when coupled to NMM$_1$ (Fig 2) are given in S1 and S2 Figs.

## Modeling plasticity with coupled neural mass models

We consider two populations NMM$_1$ and NMM$_2$ modelled by the NMM formulation given in (10)-(11). NMM$_2$ receives excitatory inputs from NMM$_1$ through AMPAergic and NMDAergic signaling. Let $F_p^{(1)}$ be the firing rate of the P subpopulation of the NMM$_1$, $r^{(i)}(t)$ and $u^{(i)}(t)$ are the variables describing the presynaptic neurotransmitter release (6). Then the AMPAergic EPSP is given by

$$ \frac{d^2 y_{\text{AMPA}}}{d^2 t} = \frac{r^{(1)} u^{(1)} A_{\text{AMPA}}^{(2)}}{\tau_{\text{AMPA}}} F_p^{(1)} - \frac{2}{\tau_{\text{AMPA}}} \frac{dy_{\text{AMPA}}}{dt} - \frac{1}{\tau_{\text{AMPA}}^2} y_{\text{AMPA}}, \quad (12) $$

where $τ_{AMPA}$ is the synaptic time constant and $A_{AMPA}^{(2)}$ is the amplitude that is modulated by STP variables $r^{(1)}(t)$ and $u^{(1)}(t)$. The NMDAergic-EPSP is given by

$$ \frac{d^2 y_{\text{NMDA}}}{d^2 t} = \frac{r^{(1)} u^{(1)} A_{\text{NMDA}}^{(2)}}{\tau_{\text{NMDA}}} F_p^{(1)} - \frac{2}{\tau_{\text{NMDA}}} \frac{dy_{\text{NMDA}}}{dt} - \frac{1}{\tau_{\text{NMDA}}^2} y_{\text{NMDA}}, \quad (13) $$

where $τ_{NMDA}$ is the synaptic time constant and $A_{NMDA}^{(2)}$ is the amplitude. The NMDAergic EPSP is slower than the AMPAergic EPSP [47], thus $τ_{AMPA} < τ_{NMDA}$. The net depolarization of the postsynaptic subpopulation due to the activation of NMDAR depends on the postsynaptic potential, $V_p^{(2)}$ as well as the presynaptic entry $C_{NMDA}y_{NMDA}$, i.e. $C_{NMDA}y_{NMDA}H(V_p^{(2)})$, with $C_{NMDA}$ is the coupling strength.

The pathological plasticity is activated when both the probability of vesicle opening $U_s$ and the AMPAR insertion $\tilde{C}_{AMPA}$ are potentiated (9). Under this condition, an afferent epileptic spike from NMM$_1$, which releases high levels of glutamate, activates the extrasynaptic NMDAR. The EPSP mediated by extrasynaptic NMDR, $y_{NMDA,ext}$, has slower kinetics than the EPSP mediated by synaptic NMDAR [48] (i.e. $τ_{NMDA} < τ_{NMDA,ext}$) and $y_{NMDA,ext}$ is given by

$$ \frac{d^2 y_{\text{NMDA,ext}}}{d^2 t} = \frac{r^{(1)} u^{(1)} A_{\text{NMDA,ext}}^{(2)}}{\tau_{\text{NMDA,ext}}} F_p^{(1)} - \frac{2}{\tau_{\text{NMDA,ext}}} \frac{dy_{\text{NMDA,ext}}}{dt} - \frac{1}{\tau_{\text{NMDA,ext}}^2} y_{\text{NMDA,ext}}, \quad (14) $$

with the postsynaptic depolarization $C_{NMDA}y_{NMDA,ext}H(V_p^{(2)})$. The net postsynaptic polarization of the P subpopulation of the NMM$_2$ then reads,

$$ V_p^{(2)} = \Sigma(C_{k,p}^{(2)}y_k^{(2)}) + \tilde{C}_{AMPA}y_{AMPA} + C_{NMDA}(y_{NMDA}H(V_p^{(2)}) + y_{NMDA,ext}H(V_p^{(2)})) $$

with the plastic $\tilde{C}_{AMPA}(t)$ given in (9) and constant $C_{NMDA}$.
The changes in GABAergic activity caused by the activation of extrasynaptic NMDAR is modeled via an auxiliary variable $K(t)$,

$$ \frac{dK}{dt} = \frac{1}{\tau_K} (-K(0.5 - K)(1 - K) - C_{NMDA}y_{NMDA,ext}H(V_P^{(2)})). \quad (15) $$

When extrasynaptic NMDAR is inactive (i.e. $y_{NMDA,ext} = 0$), the Eq 15 has three equilibrium points at $K^* = \{0,0.5,1\}$. The equilibria at $K^* = \{0,1\}$ are stable, but the equilibrium at $K^* = 0.5$ is unstable. The extrasynaptic NMDAR activation (i.e. $y_{NMDA,ext} \neq 0$) can lead to a saddle-node bifurcation of the equilibria $K^* = 0.5$ and $K^* = 1$ leaving (15) with a single equilibrium point at $K^* \approx 0$, which becomes the global attractor of (15). Inactivation of the extrasynaptic NMDAR reintroduces the three equilibria but the transition from $K \approx 0$ to $K \approx 1$ state would not be possible in the model unless $y_{NMDA,ext}$ is negative. The latter is not biologically plausible or computationally possible because NMDA-mediated currents are excitatory. Therefore, the transition from $K \approx 0$ to $K \approx 1$ is irreversible. This phenomenological formulation is based on experimental studies in animal models of epilepsy, which suggest that when the secondary focus “matures”, abolition of the primary focus or disruption of cortical connections does not abolish the secondary epileptic zone [27,49].

The variable $K(t)$ modulates the excitability of NMM$_2$ as,

$$ b_{thr}^{(2)}(t) = b_{thr}^{(2)} - k_B^{(2)}(1 - K(t)), $$

where the parameter $k_B^{(2)}$ scales the impact of $K(t)$ on the excitability. Since the gamma-band activity observed at the seizure onset can be obtained by increasing the IPSP amplitude of PV interneurons [40,46], $K(t)$ modulates $G^{(2)}$ as

$$ G^{(2)}(t) = G^{(2)} + k_G^{(2)}(1 - K(t)), $$

where the parameter $k_G^{(2)}$ scales the impact of $K(t)$ on the IPSP amplitude. We assumed that, in the absence of extrasynaptic NMDAR activity, $K(t) = 1$, for which IPSP amplitude is too small for yielding a gamma-band activity in (10) and the equilibrium point of (11) is far from the critical point of the seizure transition (i.e. the $n^{(2)}$-nullcline intersects the $B^{(2)}$-nullcline on its right branch). As $K(t) \to 0$ with the activation of the extrasynaptic NMDAR, $G^{(2)}(t)$ increases to a suitable level for obtaining gamma-band activity, whereas $b_{thr}^{(2)}(t)$ decreases and drives the equilibrium point of (10) towards the critical point of seizure transition in (10)-(11). In other words, how far (11) is from a Hopf bifurcation, hence how far (10) is from seizing, is controlled by the extrasynaptic NMDAR activity.

The parameter values of Eqs 12–15 are given in Table 2. We followed the order of magnitude between different receptor types [47,48]. The amplitude of PSPs and coupling strengths are chosen to achieve self-consistency. Stochastic differential equations were iterated using Euler-Maruyama method with a step size $dt = 10e-5$ second. Simulation files are available at https://github.com/elifkoksal/plasticNMM. Bifurcation analysis was done with AUTO-07p [50].

## Reduced system of physiological plasticity

The system of equations describing the physiological plasticity is a multiple timescale system with ($\tau_d, \tau_f) < \tau_{ca} < \tau_p < \tau_U = \tau_{CAMPA}$ (Table 1). We benefit from the multiple timescale structure of (9) to reduce the complexity of the full model to a lower dimensional system of equations for analyzing the synaptic plasticity in the long term. We use the elements from singular perturbation theory to decompose the system into fast and slow subsystems and focus on the
dynamic in the slow regime [51]. We assume that the slow regime is governed by the variables concerning LTD/LTP and consolidation, e.g. ($\rho$, $U_s$, $\tilde{C}_{AMPA}$). In the slow regime, the fast variables describing STP, $r^{(1)}$ and $u^{(1)}$, are controlled by the slow variables, as:

$$r^{(1)s} = \frac{1}{1 + \tau_d u^{(1)s} F_p^{(1)}},$$

$$u^{(1)s} = \frac{U_s(1 + \tau_f F_p^{(1)})}{1 + \tau_f U_s F_p^{(1)}}$$

where superscript (.)s indicates the slow regime. The NMDA-mediated calcium dynamics in the slow regime is given by

$$[Ca]^s = \tau_{ca} G(F_p^{(1)}) h_{ca} H(V_p^{(2)}) = \tau_{ca} C_{NMDA} y_{NMDA}(t) h_{ca} H(V_p^{(2)}).$$

If the NMDAergic-EPSP $y_{NMDA}(t)$ is approximated by $y_{NMDA} \approx \tau_{NMDA} A_{NMDA}^{(2)} r^{(1)s} u^{(1)s} F_p^{(1)}$ and the AMPAergic-EPSP $y_{AMPA}(t)$ by, $y_{AMPA} \approx \tau_{AMPA} A_{AMPA}^{(2)} r^{(1)s} u^{(1)s} F_p^{(1)}$, the postsynaptic potential $V_p^{(2)}$ reads

$$V_p^{(2)} = \Sigma(C_{k,P}^{(2)} y_k^{(2)}) + \tilde{C}_{AMPA} \tau_{AMPA} A_{AMPA}^{(2)} r^{(1)s} u^{(1)s} F_p^{(1)} + C_{NMDA} \tau_{NMDA} A_{NMDA}^{(2)} r^{(1)s} u^{(1)s} F_p^{(1)}.$$

For the NMDA-mediated calcium entry and postsynaptic depolarization to occur, the postsynaptic membrane should be sufficiently depolarized, i.e. $V_p^{(2)} > V_{th}$. This depolarization depends on the AMPAergic input. Assuming that the postsynaptic system is balanced for $F_p^{(1)}$ and intra-population interactions are independent of the external input, i.e. $\Sigma(C_{k,P}^{(2)} y_k^{(2)}) = 0$, the postsynaptic polarization in response to the AMPAergic input from $NMM_1$ can be reduced to

$$V_p^{(2)} \approx \tilde{C}_{AMPA} \tau_{AMPA} A_{AMPA}^{(2)} r^{(1)s} u^{(1)s} F_p^{(1)},$$

and the calcium concentration to

$$[Ca]^s = C_{NMDA} \tau_{NMDA} A_{NMDA}^{(2)} r^{(1)s} u^{(1)s} F_p^{(1)} h_{ca}, \text{ if } V_p^{(2)} > V_{th}, \text{ otherwise } [Ca]^s \approx 0.$$

As the fast variables ($r^{(1)}, u^{(1)}$), $[Ca]$ are expressed as a function of the slow variables, ($\rho$, $U_s$, $\tilde{C}_{AMPA}$), the reduced (slow) system reads:

$$\frac{d\rho}{dt} = -\rho(1-\rho)(\rho_* - \rho) + (1-\rho)\Omega_p([Ca]) - \rho\Omega_d([Ca]), \quad (16A)$$

$$\frac{dU_s}{dt} = \frac{\tau_p}{\tau_U} (U_s^d - U_s + \rho(U_s^p - U_s^d)). \quad (16B)$$

$$\frac{d\tilde{C}_{AMPA}}{dt} = \frac{\tau_p}{\tau_{C_{AMPA}}} (C_{AMPA}^d - \tilde{C}_{AMPA} + \rho(C_{AMPA}^p - C_{AMPA}^d)), \quad (16C)$$

where $\tilde{t} = \frac{t}{\tau_p}$. If $\frac{\tau_p}{\tau_U} = \frac{\tau_p}{\tau_{C_{AMPA}}} \ll 1$, then we can take them to zero limit and treat ($U_s$, $\tilde{C}_{AMPA}$) as parameters to investigate the variation of $\rho$.
# Results

## Physiological plasticity

The model of physiological plasticity includes presynaptic STP, calcium-driven pre- and post-synaptic LTP/LTD and consolidation ((9), Fig 3A). These mechanisms run in different timescales, with STP being the fastest, and consolidation being the slowest. In Section Reduced system of physiological plasticity, we applied singular perturbation theory to obtain a reduced system of physiological plasticity. Here we assume that the timescales of calcium-driven long-term plasticity and consolidation are different enough ($\tau_{\rho} \ll \tau_{U} = \tau_{C_{AMPA}}$), to investigate the variation of $\rho$, which, in turn, controls the long-term variations in the neurotransmitter release probability, $U_s$, and the AMPARergic coupling strength, $\tilde{C}_{AMPA}$, by treating these two slow variables as parameters. For a simple investigation of $\rho$ (16a) for $\frac{\tau_{\rho}}{\tau_{U}} = \frac{\tau_{\rho}}{\tau_{C_{AMPA}}} = 0$, we assume that $\tilde{C}_{AMPA} \approx C_{AMPA}^d + \kappa U_s$, ($\kappa > 0$), which is a reasonable assumption since both $\tilde{C}_{AMPA}$ and $U_s$ depend on $\rho$ in the same manner. Fig 4A and 4B show $V_p^{(2)}$ and $[Ca]^s$ as a function of $F_p^{(1)}$ and $U_s$ with the parameter set given in Table 1. The postsynaptic subpopulation depolarizes monotonically for increasing $F_p^{(1)}$ and $U_s$ (Fig 4A). The change in the calcium concentration follows a sigmoidal shape (Fig 4B). It remains close to zero for low values of $U_s$ and $F_p^{(1)}$. It increases sharply for medium values of $U_s$ and $F_p^{(1)}$, then slowly for higher values of $U_s$ and $F_p^{(1)}$ as it approaches a plateau.

We then compute the equilibrium points of (16) for $\frac{dU_s}{dt} = \frac{d\tilde{C}_{AMPA}}{dt} = 0$. In the absence of any pre- or postsynaptic activity, (16a) has three equilibrium points at $\rho^* = \{0,0.5,1\}$ with $\rho^* = \{0,1\}$ being stable and $\rho^* = 0.5$ unstable. If the variations in $U_s$ and $F_p^{(1)}$ are not sufficient to change the number of equilibrium points of (16A), the depressed synapses will remain depressed, and the potentiated ones will remain potentiated. Otherwise, (16A) has a single stable equilibrium point, such that, the synapses will either depress or potentiate. Fig 4C shows $\rho^*$ as a function of

<figure><img src="image_6.png" /><figcaption>**Fig 3. Physiological plasticity at presynaptic and postsynaptic sites.** (A) Presynaptic glutamatergic release and subsequent postsynaptic depolarization results in calcium influx via NMDARs. Calcium signaling activates independent biochemical pathways, leading to postsynaptic and presynaptic long-term changes due to the insertion of new AMPARs and an increase in vesicle release probability. Under physiological conditions, excitation is balanced by sufficient chloride trafficking through GABAR and KCC2. (B) Pathological plasticity at postsynaptic site. When the released glutamate is too high to be used by synaptic AMPAR and NMDAR, it diffuses to the extrasynaptic site and activates the extrasynaptic NMDAR. Activation of NMDAR causes internalization of extrasynaptic GABAR. Calcium influx via the extrasynaptic NMDAR activates calpain that perturbs the KCC2 functioning.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012666.g003
<figure><img src="image_7.png" /><figcaption>**Fig 4.** Postsynaptic depolarization, calcium concentrating and long-term plasticity in the slow regime. (A) Postsynaptic polarization $V_p^{(2)}$ and (B) calcium concentration $[Ca]^s$ as a function of $(F_p^{(1)}, U_s)$ for the parameter set given in Table 1 and $\gamma_p = 2$. (C) Equilibrium points of (16a) as function of $(F_p^{(1)}, U_s)$ for different values of $(\theta_p, \theta_d)$. In the blue regions (16a) has three equilibrium points at $\rho_{1,2,3}^* \in [0, 1]$, a single equilibrium point at $\rho^* = 0$ in the salmon regions, and a single equilibrium point at $0.5 < \rho^* \le 1$ in the red regions.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012666.g004
$(U_s, F_p^{(1)})$ as the LTD and LTP thresholds, $\theta_d$ and $\theta_p$, vary on the calcium concentration surface (Fig 4B). The dark blue regions indicate that (16a) has three equilibrium points at $\rho^* = \{0,0.5,1\}$, therefore, the synapses do not change. When the potentiation threshold $\theta_p$ is too high ($\theta_p=1.5$, upper panels), the depressed state (salmon regions) is the only attractor for moderate to high input firing rates $F_p^{(1)}$ and release probability $U_s$ depending on $\theta_d$. Lowering the potentiation threshold $\theta_p$ introduces a range of $(U_s, F_p^{(1)})$ for which the potentiated state is an attractor (red regions in Fig 4C for $\theta_p = \{0.4,0.7\}$). Consequently, both depressed and potentiated synapses are preserved for a low presynaptic activity. Lowering the potentiation threshold expands the range of LTP and, increasing $\theta_d$ shrinks the LTD range.

Physiological plasticity under epileptic activity

In our model, long-term plasticity depends on the time lag between the interictal spikes of NMM₁ and NMM₂, as well. Here, the term “spike” refers to shape of the model output assumed as the sum of the input signals to the subpopulation P of NMMᵢ (see Section The revisited neural mass model). In this section we assume that only the physiological synaptic plasticity mechanism (Fig 3A) is active. Fig 5 exemplifies how $\rho$ changes in response to single spikes of NMM₁ and NMM₂ elicited by a pulse stimulation when $\rho$ is initiated at $\rho(0) = 0.5$. The synaptic variable $\rho$ potentiates regardless of the pulse delay for $(\theta_d, \theta_p) = (0.6,0.7)$, depresses $(\theta_d, \theta_p) = (0.1,0.7)$ (Fig 5A). We obtain a depression-potentiation-depression curve of $\Delta\rho$ for $(\theta_d, \theta_p) = (0.3,0.7)$, for which the system response is shown in Fig 5B–5D. When NMM₂ spikes $\Delta t = 0.02$ sec before NMM₁, the calcium concentration stays between $(\theta_d, \theta_p)$ and $\rho$ decreases (Fig 5B). When NMM₂ and NMM₁ spike simultaneously, the calcium concentration exceeds $\theta_p$ and $\rho$ increases (Fig 5C). Finally, when NMM₂ spikes $\Delta t = 0.02$ sec after NMM₁, $\rho$ decreases since the calcium concentration stays between $(\theta_d, \theta_p)$ (Fig 5D).

Depending on the potentiation and depression thresholds, and intrinsic dynamics of the interacting populations, aperiodic but continuous interictal epileptic discharges can cause long-term changes in the presynaptic release probability and coupling strength between the interacting populations (Fig 6). Here, we assume that NMM₁ is in interictal spiking mode as in Fig 2B, and the less excitable NMM₂ responds to NMM₁ by generating spikes. If the potentiation threshold is high, then initially depressed synapses between NMM₁ and NMM₂ remain depressed. For a lower potentiation threshold, the depolarization in NMM₂ is high enough to increase the calcium concentration above $\theta_p$, hence causing a positive change in the synaptic variable $\rho(\Delta\rho>0)$. If this regime is maintained long enough and $(\theta_d, \theta_p)$ are in a suitable range for LTP, $\rho$ can cross the threshold separating depressed to potentiated states, and a pre- and postsynaptic LTP occurs (Fig 6B). Because of the LTP, the time lag between the NMM₁ and NMM₂ spikes decreases. In the same configuration but with potentiated synapses, increasing the difference between $\theta_d$ and $\theta_p$ causes pre- and postsynaptic LTD since the calcium concentration remains in $(\theta_d, \theta_p)$ (Fig 6C). Consequently, the time lag between the NMM₁ and NMM₂ spikes increases. If both depression and potentiation thresholds are high, despite the interaction between NMM₁ and NMM₂, synaptic variables do not change.

Experimental studies suggest that a high rate of synchronized activity, during an epileptic seizure for instance, increases the glutamate release. Indeed, the model also suggests that an increased firing rate during a seizure in NMM₁ can drive the calcium concentration above $\theta_p$, and cause LTP, while the interictal activity with the same values of $(\theta_d, \theta_p)$ does not change synaptic efficacy (Fig 6A vs Fig 7A). Alternatively, if calcium concentration is already above $\theta_p$, then the increased firing rate can accelerate LTP in NMM₂ (Fig 6B vs Fig 7B). In Fig 7 the repeated seizures in NMM₁ are triggered by a stochastic input to the $B^{(1)}$ variable.
<figure><img src="image_8.png" /><figcaption>**Fig 5.** Change in synaptic efficacy under epileptic spikes. (A) Rate of change in the synaptic efficacy variable $\rho$ under single epileptic spike in NMM₁ and NMM₂ when the system (9)-(11) is initialized from $\rho(0) = 0.5$ with $(U_s, \tilde{C}_{\text{AMPA}}) = (U_s^d, C_{\text{AMPA}}^d)$. Epileptic spikes are triggered by a pulse stimulation applied with a time delay of $\Delta t$ (sec) for $(\theta_d, \theta_p) = (0.6,0.7)$ (green line), $(\theta_d, \theta_p) = (0.3,0.7)$ (orange line) and $(\theta_d, \theta_p) = (0.1,0.7)$ (blue line). (B-D) Pre- and post-spiking is triggered by a single pulse with $\Delta t = -0.02$ sec in (B), $\Delta t = 0$ in (C) and $\Delta t = 0.02$ in (D) for $(\theta_d, \theta_p) = (0.3,0.7)$. The corresponding behaviour of [Ca] and $\rho$ are shown in the middle and bottom panels, respectively. The orange dotted and green dotted lines in the middle panels mark $\theta_d$ and $\theta_p$, respectively.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012666.g005
<figure><img src="image_9.png" /><figcaption>**Fig 6. Physiological plasticity under interictal epileptic discharges.** (A) Potentiation of initially depressed synapses for ($\theta_{d}, \theta_{p}$) = (0.3,0.4). (B) Depression of initially potentiated synapses for ($\theta_{i,d}, \theta_{p}$) = (0.1,1.5).</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012666.g006

Expansion of the epileptogenic network

Physiological plasticity alone increases coupling strength by increasing the probability of neuro-
transmitter release at the presynaptic site and inserting AMPAR at the postsynaptic site. Unless
the extrasynaptic NMDAR are activated at the postsynaptic site, the postsynaptic population
can preserve its excitation/inhibition balance, hence, does not undergo seizures. However, the
extrasynaptic NMDAR activation triggers different processes that disturb the excitation/inhibi-
tion balance by reducing the effect of the GABARergic activity [7,52,53] (Fig 3B).

We assumed that the extrasynaptic NMDAR (14) subtype is activated after the synapses are
potentiated, i.e. for $U_s>0.7$. The extrasynaptic NMDAR activity changes the amplitude of the
PV interneurons and the excitability of NMM₂ via the auxiliary variable $K(t)$ (15). We investi-
gate the impact of these changes in the whole model given by (9)-(15). We simulate the whole
system (9)-(15) starting from the last point of the solution in Fig 7B for ($\theta_{d}, \theta_{p}$) = (0.3,0.4).
When NMM₁ undergoes seizures, and NMM₂ responds by epileptic spikes since NMM₂ has a
low excitability. Activation of the extrasynaptic NMDAR drives the variable $K(t)$ to 0 as it
causes “structural” changes within NMM₂. In particular, the amount of the GABAergic loss in
NMM₂, which is represented by the parameter $k_B^{(2)}$, can decrease the excitability threshold of
NMM₂ and cause seizures in NMM₂. The relation between the GABAergic loss and seizure in
NMM₂ is represented in the phase space of ($B^{(2)}$, $n^{(2)}$) in Fig 8A, and the time traces are shown
in Fig 8B. For $0 \le k_B^{(2)} < 9.5$, the equilibrium point of the ($B^{(2)}$, $n^{(2)}$)-subsystem remains on the
<figure><img src="image_10.png" /><figcaption>**Fig 7.** Impact of seizure to physiological plasticity. Potentiation of depressed synapses for $(\theta_d, \theta_p) = (0.3,0.7)$ (A) and for $(\theta_d, \theta_p) = (0.3,0.4)$ (B).

https://doi.org/10.1371/journal.pcbi.1012666.g007</figcaption></figure>

right of the critical point of the seizure transitions (i.e. right fold of the $B^{(2)}$-nullcline). For $0 \le k_B^{(2)} < 7.4$, NMM$_2$ responds to the seizures of NMM$_1$ by interictal spikes only. For $7.4 \le k_B^{(2)} < 9.5$, the fluctuations in $K(t)$ due to the activity of NMM$_1$ cause $B^{(2)}$ to cross the critical point and to jump to the left branch of the $B^{(2)}$-nullcline, which correspond to an ictal regime in NMM$_2$. After following the left branch of the $B^{(2)}$-nullcline, $B^{(2)}$ jumps back to the right branch of the $B^{(2)}$-nullcline, hence NMM$_2$ to the interictal regime. For $9.5 \le k_B^{(2)} < 10.5$, NMM$_2$ generates spontaneous and periodic seizures as the equilibrium point of the $(B^{(2)}, n^{(2)})$-subsystem passes to the left of the critical point of the seizure transition.

One of the markers of the seizure onset zone in the case of focal epilepsy is the gamma-band activity observed at the seizure onset [45]. Parameter $k_G^{(2)}$ scales the impact of the variable $K(t)$, hence of the activation of extrasynaptic NMDAR, on the fast GABAergic interneurons. NMM$_2$ generates a gamma-band activity at the seizure onset for $k_G^{(2)} > 8$ as demonstrated in the power spectral density (PSD) of $V_p^{(2)}$ during a seizure triggered by NMM$_1$ for $k_B^{(2)} = 9$ (Fig 8C). Fig 8E and 8F show the response of NMM$_2$ for $(k_B^{(2)}, k_G^{(2)}) = (0, 0)$ and $(k_B^{(2)}, k_G^{(2)}) = (9, 20)$ for the epileptic spikes and seizure of NMM$_1$, respectively. NMM$_2$ exhibits a seizure in $t = [180,220]$ after spiking as a response to the first seizure in NMM$_1$ in $t = [160,197]$, and a second one in $t = [412,451]$ as a response to the second seizure in NMM$_1$ in $t = [408,447]$. The time delay between the first seizure in NMM$_1$ and NMM$_2$'s response to it $(\Delta t_{seizure1} = t_{seizure1}^{(1)} - t_{seizure1}^{(2)})$ varies between -25 sec and -20 sec, whereas the delay between the second seizure in NMM$_1$ and NMM$_2$'s response $(\Delta t_{seizure2} = t_{seizure2}^{(1)} - t_{seizure2}^{(2)})$ varies between -16 sec and -1 sec (Fig 8G). This variation in the response times is due to the difference between the pathological plasticity stages: the first seizure occurs in the early stage (before the transition from $K \approx 0$ to $K \approx 1$ in (15)), whereas the second one is in the late stage (after the transition from $K \approx 0$ to $K \approx 1$ in (15)), highlighting the progressive nature of the condition and its impact on neural dynamics. We also notice that $\Delta t_{seizure}$ decreases with increasing $k_B^{(2)}$ (Fig 8G), hence as NMM$_2$ loses its GABAergic integrity.
<figure><img src="image_11.png" /><figcaption>**Fig 8.** Effect of the pathological plasticity on NMM₂. (A) Phase portrait of (B⁽²⁾,n⁽²⁾)-subsystem showing the effect of the parameter *k*<sub>B</sub>⁽²⁾ on the excitability of NMM₂. NMM₂ crosses the critical point as *k*<sub>B</sub>⁽²⁾ increases. (B) Time traces of B⁽²⁾ shown in (A). (C) Normalized PSD of V<sub>P</sub>⁽²⁾ during seizures shows the effect of *k*<sub>G</sub>⁽²⁾ for *k*<sub>B</sub>⁽²⁾ = 9. Gamma-band activity appears for *k*<sub>G</sub>⁽²⁾ > 8. (D-E) Example time trace of (9)-(15). (D) Dynamics of NMM₁ that exhibits two seizures in *t* = [160,197] and in *t* = [408,447]. (E) Dynamics of NMM₂ for (*k*<sub>B</sub>⁽²⁾, *k*<sub>G</sub>⁽²⁾) = (0, 0). (F) Dynamics of NMM₂ for (*k*<sub>B</sub>⁽²⁾, *k*<sub>G</sub>⁽²⁾) = (9, 20), where NMM₂ exhibits two seizures in *t* = [180,220] and in *t* = [412,415]. (G) Time differences between the seizures exhibited by NMM₂.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012666.g008

Activity in the secondary focus after silencing the epileptic zone

The activation of extrasynaptic NMDAR induces irreversible structural changes in NMM₂, in particular, in its GABAergic structure. How significant these changes are can be observed from the difference between the activities of NMM₁ and NMM₂. For this we compute the spike frequencies of NMMs before the first seizure in NMM₁ in *t* = [0,160], which corresponds
<figure><img src="image_12.png" /><figcaption>**Fig 9. Resection of the primary focus.** (A) Spike frequencies of NMM₁ and NMM₂ for different values of $k_B^{(2)}$ before NMM₁'s first seizure between $t = [0,160]$ and between NMM₁'s first and second seizures in $t = [210,400]$ in Fig 8. The simulations are continued starting from $t = 500\text{sec}$ of the solutions in Fig 8. After the cessation of any form of activity in NMM₁ at $t = 500\text{sec}$ (B), NMM₂ stops exhibiting seizures for $k_B^{(2)} = 9$ (C) but continues spontaneous seizures for $k_B^{(2)} = 10$ (D).</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012666.g009

to the beginning of the pathological plasticity process, and between the first and second seizures in NMM₁ in $t = [200,400]$, which corresponds to the late phase of the pathological plasticity process. The spike frequency of NMM₂ is computed for the intervals $t = [0,160]$ and $t = [210,400]$, so after its response to the NMM₁'s seizures (Fig 9A). The spike frequencies of NMM₁ in these time windows remain intact. The spike frequencies in NMM₂ are indistinguishable from the ones of NMM₁ for $0 \le k_B^{(2)} \le 7$, but they are much higher for $k_B^{(2)} = \{8, 9, 10\}$. This is because as $k_B^{(2)}$ increases, the excitability threshold of NMM₂ decreases and the amplitude of fluctuations in $K(t)$ increases. These two actors trigger additional spikes (Fig 8E) and seizures in NMM₂ (Fig 8B for $k_B^{(2)} = 10$). Then the question is what happens to NMM₂ if NMM₁ is silenced. To answer this question, we reinitiate the system (9)-(15) from the last points of the solutions in Fig 8 for $A^{(1)} = 0$ (Fig 9B–9D). Removing the input from NMM₁ to NMM₂ abolishes interictal spikes in NMM₂ both for $k_B^{(2)} = 9$ and $k_B^{(2)} = 10$ and seizures for $k_B^{(2)} = 9$. If the GABAergic loss is important, as it is for $k_B^{(2)} = 10$, stochastic inputs to $B^{(2)}(t)$ can trigger spontaneous seizures NMM₂, despite a completely silent NMM₁.
Discussion

In this study, we presented a NMM of physiological and pathological plasticity under epileptic activity. Physiological plasticity takes into account short- and long-term variations in the glutamatergic interactions at both presynaptic and postsynaptic sites. Pathological plasticity considers the mechanisms that can contribute to abnormal changes resulting in exaggerated hyperexcitability associated with epileptic activity, such as GABAergic modulations at the postsynaptic site. Under epileptic activity, increased connectivity strength resulting from the physiological plasticity triggers the pathological plasticity, which can cause a secondary focus. These mechanisms have been implemented in an autonomous (self-driven) NMM, which generates the typical focal seizure pattern from interictal to seizure termination.

In neurobiological studies, spillover has been observed during synchronous activation of a large number of glutamatergic fibers [54]. NMMs, by definition, represent synchronous activity of a neural population and therefore synchronicity is implicitly coded in the model [55]. Physiological synaptic plasticity in the model integrates three main mechanisms. The first one is the presynaptic STP. The second one is the calcium-mediated long-term synaptic plasticity that changes the probability of glutamate release at the presynaptic site and insertion of AMPAR at the postsynaptic site. The third mechanism is the consolidation of the long-term plasticity. Although the description of consolidation remained phenomenological as in [30,32], its integration into the model allowed for a sustained change in synaptic connectivity.

We assumed that presynaptic and postsynaptic long-term plasticity are proportional. While this simplified assumption is true for retrograde nitric oxide signaling, retrograde endocannabinoid release works in the opposite direction. Costa et al. [37] suggested a phenomenological model integrating the effects of opposing proteins on long-term presynaptic plasticity. Such effects can be integrated in our model, for instance, by introducing a second variable to track the calcium variations, which modifies the neurotransmitter release probability.

Our approach of modeling NMDA currents is quasi-physiological, similar to [38], but without explicitly modeling the magnesium gate. Yet, it does include the essential ingredients mediating the NMDAR activity, which are the presynaptic release and postsynaptic potentiation. Conductance-based modeling of synaptic receptors in the neural mass formulation has been proposed in [56], and more recently, in [57]. While the conductance-based neural mass formulation may provide a more physiological way of modeling neural dynamics, we argue that our approach is consistent with the neural mass formulation we considered and captures the key ingredients necessary for the problem studied. Interictal epileptic discharges in the presynaptic population and the postsynaptic response to them can induce LTP or LTD depending on the plasticity thresholds. The synchronization between pre- and postsynaptic interictal epileptic discharges increases (decreases) with LTP (LTD). Ictal activity in the presynaptic population speeds up the LTP process, and if frequent, it can cause a transition from LTD-trend to LTP-trend. LTP then triggers pathological mechanisms as the postsynaptic region receives higher glutamatergic input from the presynaptic region. The model suggests that preventing LTP by targeting both the presynaptic region (decreasing the epileptic activity) and postsynaptic region (blocking NMRA receptors, increasing the potentiation threshold, or applying GABAergic agonists, etc.) can be crucial to avoid pathological consequences.

The pathological plasticity under epileptic activity includes disruption of the GABAergic pathway caused by the activation of the extrasynaptic NMDAR. The model captures possible changes without detailing the process, such as how calcium entry through the extrasynaptic NMDAR acts on calpain, how a decreasing number of KCC2 changes the GABA reversal potential, or how GABAR expression is downregulated, etc. Rather, the impact of the extrasynaptic NMDAR activation is reflected by the auxiliary variable *K(t)* that acts on the
perisomatic and dendritic GABAergic interneurons of the NMM. The former mimics the effect of KCC2 downregulation at the level of interactions between PV interneurons and from PV interneurons to PYR cells that induces a gamma-band activity, which is considered a signature of an epileptogenic zone in focal seizures [58–60]. The latter mimics the effect of GABAR downregulation that reduces the excitability threshold of the population targeted by the epileptic population. Indeed, the dendritic SOM interneurons control the calcium dynamics, for instance by directly inhibiting calcium influx [61–63]. Low dendritic SOM projection may lead to accumulation of intracellular calcium and aggravate KCC2 dysfunctions. The impact of the extrasynaptic NMDAR activation on the GABAergic system is scaled in the model. This scaling factor can be interpreted as an intrinsic property of the region subject to epileptic activity that can maintain the GABAergic integrity. In our simulations, we observe a delay between the seizures when a seizure in the postsynaptic region is triggered by a presynaptic seizure. Such propagation delays have been reported in animal models of focal epilepsy [64], and interestingly, it can decrease during epileptogenesis [65].

In this study we extended a widely used model of the hippocampus [40] for ensuring autonomous epileptic activity from interictal to ictal and seizure termination. Essentially, this was achieved by introducing a slow subsystem that modulates the IPSP amplitude of the slow GABAergic interneurons. In other words, we introduced a slow-fast process for obtaining an autonomous interictal-to-ictal transition. Except for very few studies [59,66], the parameter variations leading to an interictal-to-ictal transition has been performed manually, in particular by varying the IPSP amplitude of the slow GABAergic interneurons. In [58,59] a dynamic chloride accumulation in pyramidal neurons caused by GABAergic activity modulates the IPSP amplitude. However, the authors of [59] do not propose a mechanism for seizure termination. In [66] this modulation was obtained through feedback from the activity of the pyramidal cells. Our approach is similar to the phenomenological approach of [66], except the feedback from the pyramidal neurons and the dimension of the slow subsystem. In the transition to seizures, not only pyramidal neurons but also interneurons can play a critical role [67]. Our model can be improved further to reflect the involvement of specific cellular and network mechanisms by considering appropriate feedback mechanisms. Furthermore, the fast subsystem representing neuronal subpopulations and the slow subsystem controlling the excitability and the IPSP amplitude have the essential ingredients for studying different bifurcation types, which can be inherent to specific seizure types (absence seizures, focal seizures) with different seizure onset patterns [45]. Furthermore, the slow-fast system can be used to study the effect of electrical stimulation that is a general procedure in clinics to understand the excitability of the brain region to identify epileptogenic networks.

Our study of secondary epileptogenesis confirms the so-called *seizures can beget seizures* phenomena: epileptic activity in one brain region can recruit healthy regions and cause the formation of a complex epileptic network [68]. The model suggests that epileptic seizures boosts physio-pathological plasticity, as suggested in the literature [8,69]. Once a healthy zone is driven to the critical point, then it generates epileptic discharges and/or seizures. The phenomenon of secondary epileptogenesis was first demonstrated in animal models [49], including: (1) the electrical stimulation kindling-like model in frogs, rabbits, rats, guinea pigs, and cats; (2) the kainic acid model in rats and mice; and (3) several others [27]. In humans, this phenomenon can be observed in clinical practice [27], and it can be related to the organization of epileptic networks. The network concept in epilepsy can be a key factor in identifying the anatomic distribution of the epileptogenic process and in the clinical expression of dynamic course of seizures [70,71]. For example, development of the epileptic networks in temporal lobe epilepsy is correlated with time (epilepsy duration) and has been argued to be due to secondary epileptogenesis processes [72,73]. Therefore, the timing of interference with epileptic
activity may be critical [74]. Indeed, in clinical practice, acute symptomatic seizures are treated
with anti-seizure medication [75]. Furthermore, the secondary epileptogenesis is not only lim-
ited to the “creation” of a secondary focus or the expansion of an epileptic network but also to
the expansion of the dysfunctional tissues. For instance, the epileptic zone can be larger that
MRI-visible lesion, which can require larger resections than the lesion [76,77], and related the
abnormalities in surrounding cortex [78]. In addition to the role of epileptic activity in epilep-
togenesis, epileptic activity can interfere with physiological oscillatory activity in distributed
neural networks and cause cognitive comorbidity [79–83]. Generalizing our model to person-
alized dynamic brain network modeling for planning the treatment of epilepsy is a future
research direction [84].

On the other hand, epilepsy is a complex and patient-specific disease, and it is challenging to propose a “one-fits-all” theory. The “seizures beget seizures” phenomenon has been questioned recently [85,86]. The correlation between epilepsy duration and number of high epileptogenic structures was not reported for every case. This suggests that mechanisms of epileptic network’s extension (and most likely of secondary epileptogenesis) rely on many factors, such as the intrinsic properties of underlying lesion and the anatomical substrate of the location of the pri- mary epileptogenic zone. Our model has many limitations to address this challenge. First, it only considers glutamatergic plasticity mediated by epileptic activity between unidirectionally coupled neuronal populations. A growing body of literature suggests that inhibitory GABAergic synapses exhibit long-term plasticity, which can have a protective role against the expansion of epileptogenic networks [87,88]. Local GABAergic mechanisms, such as blanket inhibition [89,90] and local circuity [91], can serve to balance excitation and prevent epilepsy. The model does not address mGLUR mediated plasticity [92,93], homeostatic effects [94,95] or other mechanisms, such as epigenetic changes [96]. Another limiting factor can be the structural and pathological differences between the primarily and secondary regions. Indeed, the impact of the pathological changes are scaled in the model, which can be considered to study the propagation of the epileptic activity in a network of connected nodes with different levels of GABAergic integrity. Yet here, we restricted ourselves to a minimal network of glutamatergic feedforward connectivity, where the projection is unidirectional on excitatory neuronal subpopulation. The glutamatergic interactions between brain regions are often bidirectional targeting both excit- atory neurons and inhibitory interneurons. Inhibitory populations can also have long-range inhibitory projections [97,98]. These inhibitory/excitatory bidirectional interactions would eventually impact the network dynamics under epileptic activity. In addition, the considered NMM is the simplest formulation for the CA1 region of the hippocampus. A laminar NMM can be considered for the distant synaptic rules related to dendritic location, which can be useful for studying the plasticity induced in neocortical regions by non-invasive brain stimulation. Modeling experimental recording could be an asset. Finally, the model simulates epileptiform discharges, whereas the physiological plasticity is also mediated by theta rhythm and theta- gamma oscillations in the hippocampus [99]. These points will be addressed in future works.

Conclusion

We have proposed a conceptual framework based on neural mass modeling for studying physio-pathological plasticity under epileptic activity. We have shown how strengthened connectivity between an epileptic and non-epileptic brain regions can lead to lead to GABAergic dysfunctions in the latter, and it can form an epileptogenic network. Our study can benefit the development of *plastic* large-scale brain models for studying the expansion of epileptogenic networks, formation of secondary foci, and the design of invasive or non-invasive brain stimulation protocols to control the epileptogenic networks.
**Supporting information**

S1 Fig. Dynamics of NMM₁ and NMM₂ during the interictal phase under a unidirectional interaction from NMM₁ to NMM₂ for the parameter set given in Table 2 of the main text.

(A) Bifurcation diagram of NMM₁ where the amplitude of $y_p^{(1)}$ is presented as a function of $B^{(1)}$. The blue curve shows the branch of equilibrium points (bold for stable and dashed for unstable equilibrium points). The red curves show the amplitude of $y_p^{(1)}$ in the oscillatory regime. The Hopf bifurcations along the branch of equilibrium points are denoted by red dots and saddle-node bifurcation by blue dots. The dynamical regimes that correspond to the fast onset, ictal and interictal periods are marked by purple, yellow and cyan patches, respectively. The time solution (black curve) is superimposed on the bifurcation diagram. (B) Time trace for $V_p^{(1)}$ showing interictal spikes. (C) Phase plane of the $(B^{(1)},n^{(1)})$-subsystem with the $B^{(1)}$-nullcline (blue curve, bold for stable and dashed for branches) and the $n^{(1)}$-nullcline (orange curve). The $B^{(1)}$ values that correspond to the dynamical regimes in (A) are marked by the same color code. The time solution (black curve) is superimposed on the phase plane. (E) Bifurcation diagram of the uncoupled NMM₂ where the amplitude of $y_p^{(2)}$ is presented as a function of $B^{(2)}$. The blue curve shows the branch of equilibrium points (bold for stable and dashed for unstable equilibrium points). The red curves show the amplitude of $y_p^{(2)}$ in the oscillatory regime. The Hopf bifurcations along the branch of equilibrium points are denoted by red dots and saddle-node bifurcation by blue dots. The dynamical regimes that correspond to the oscillatory states is marked in yellow, the steady states in orange for low values of $B^{(2)}$ and in cyan for high values of $B^{(2)}$. The time solution (black curve) is superimposed on the bifurcation diagram. (F) Time trace for $V_p^{(2)}$ showing interictal spikes. (G) Phase plane of the $(B^{(2)}, n^{(2)})$-subsystem with the $B^{(2)}$-nullcline (blue curve, bold for stable and dashed for branches) and the $n^{(1)}$-nullcline (orange curve). The $B^{(2)}$ values that correspond to the dynamical regimes in (E) are marked by the same color code. The time solution (black curve) is superimposed on the phase plane.

(TIF)

S2 Fig. Dynamics of NMM₁ and NMM₂ during the ictal phase under a unidirectional interaction from NMM₁ to NMM₂ for the parameter set given in Table 2 of the main text. (A)

Bifurcation diagram of NMM₁ where the amplitude of $y_p^{(1)}$ is presented as a function of $B^{(1)}$. The blue curve shows the branch of equilibrium points (bold for stable and dashed for unstable equilibrium points). The red curves show the amplitude of $y_p^{(1)}$ in the oscillatory regime. The Hopf bifurcations along the branch of equilibrium points are denoted by red dots and saddle-node bifurcation by blue dots. The dynamical regimes that correspond to the fast onset, ictal and interictal periods are marked by purple, yellow and cyan patches, respectively. The time solution (black curve) is superimposed on the bifurcation diagram. (B) Time trace for $V_p^{(1)}$ showing interictal spikes. (C) Phase plane of the $(B^{(1)},n^{(1)})$-subsystem with the $B^{(1)}$-nullcline (blue curve, bold for stable and dashed for branches) and the $n^{(1)}$-nullcline (orange curve). The $B^{(1)}$ values that correspond to the dynamical regimes in (A) are marked by the same color code. The time solution (black curve) is superimposed on the bifurcation diagram. (E) Bifurcation diagram of the uncoupled NMM₂ where the amplitude of $y_p^{(2)}$ is presented as a function of $B^{(2)}$. The blue curve shows the branch of equilibrium points (bold for stable and dashed for unstable equilibrium points). The red curves show the amplitude of $y_p^{(2)}$ in the oscillatory regime. The Hopf bifurcations along the branch of equilibrium points are denoted by red dots and saddle-node bifurcation by blue dots. The dynamical regimes that correspond to the oscillatory states is marked in yellow, the steady states in orange for low values of $B^{(2)}$ and in cyan for high values of $B^{(2)}$.
for high values of $B^{(2)}$. The time solution (black curve) is superimposed on the phase plane. (F)
Time trace for $V_p^{(2)}$ showing interictal spikes. (G) Phase plane of the $(B^{(2)}, n^{(2)})$-subsystem with
the $B^{(2)}$-nullcline (blue curve, bold for stable and dashed for branches) and the $n^{(1)}$-nullcline
(orange curve). The $B^{(2)}$ values that correspond to the dynamical regimes in (A) are marked by
the same color code. The time solution (black curve) is superimposed on the phase plane.
(TIF)

## Acknowledgments

Authors thank Fabrice Bartolomei (APHM, Timone Hospital, Epileptology and Cerebral Rhythmology Department, Marseille, France) for valuable discussions.

## Author Contributions

**Conceptualization:** Elif Köksal-Ersöz, Pascal Benquet.

**Data curation:** Elif Köksal-Ersöz.

**Formal analysis:** Elif Köksal-Ersöz.

**Funding acquisition:** Fabrice Wendling.

**Investigation:** Elif Köksal-Ersöz.

**Methodology:** Elif Köksal-Ersöz.

**Project administration:** Elif Köksal-Ersöz, Pascal Benquet, Fabrice Wendling.

**Resources:** Fabrice Wendling.

**Software:** Elif Köksal-Ersöz.

**Supervision:** Elif Köksal-Ersöz, Pascal Benquet.

**Validation:** Elif Köksal-Ersöz.

**Visualization:** Elif Köksal-Ersöz.

**Writing – original draft:** Elif Köksal-Ersöz.

**Writing – review & editing:** Elif Köksal-Ersöz, Pascal Benquet, Fabrice Wendling.

## References

1. Turrigiano GG, Nelson SB. Homeostatic plasticity in the developing nervous system. Nat Rev Neurosci. 2004 Feb; 5(2):97–107. https://doi.org/10.1038/nrn1327 PMID: 14735113

2. Zhang W, Linden DJ. The other side of the engram: experience-driven changes in neuronal intrinsic excitability. Nat Rev Neurosci. 2003 Nov; 4(11):885–900. https://doi.org/10.1038/nrn1248 PMID: 14595400

3. Abraham WC, Bear MF. Metaplasticity: the plasticity of synaptic plasticity. Trends Neurosci. 1996 Apr 1; 19(4):126–30. https://doi.org/10.1016/s0166-2236(96)80018-x PMID: 8658594

4. Scimemi A. Structure, function, and plasticity of GABA transporters. Front Cell Neurosci [Internet]. 2014 [cited 2023 Dec 12];8. Available from: https://www.frontiersin.org/articles/10.3389/fncel.2014.00161 PMID: 24987330

5. Hennequin G, Agnes EJ, Vogels TP. Inhibitory Plasticity: Balance, Control, and Codependence. Annu Rev Neurosci. 2017; 40(1):557–79. https://doi.org/10.1146/annurev-neuro-072116-031005 PMID: 28598717

6. Chiu CQ, Barberis A, Higley MJ. Preserving the balance: diverse forms of long-term GABAergic synaptic plasticity. Nat Rev Neurosci. 2019 May; 20(5):272–81. https://doi.org/10.1038/s41583-019-0141-5 PMID: 30837689
7. Khalilov I, Holmes GL, Ben-Ari Y. In vitro formation of a secondary epileptogenic mirror focus by inter-hippocampal propagation of seizures. Nat Neurosci. 2003 Oct; 6(10):1079–85. https://doi.org/10.1038/nn1125 PMID: 14502289

8. Ben-Ari Y. Epilepsies and neuronal plasticity: for better or for worse? Dialogues Clin Neurosci. 2008 Mar 31; 10(1):17–27. https://doi.org/10.31887/DCNS.2008.10.1/ybenari PMID: 18472481

9. Ben-Ari Y, Holmes GL. Effects of seizures on developmental processes in the immature brain. Lancet Neurol. 2006 Dec 1; 5(12):1055–63. https://doi.org/10.1016/S1474-4422(06)70626-3 PMID: 17110286

10. Citri A, Malenka RC. Synaptic Plasticity: Multiple Forms, Functions, and Mechanisms. Neuropsychopharmacology. 2008 Jan; 33(1):18–41. https://doi.org/10.1038/sj.npp.1301559 PMID: 17728696

11. Zucker RS, Regehr WG. Short-Term Synaptic Plasticity. Annu Rev Physiol. 2002; 64(1):355–405. https://doi.org/10.1146/annurev.physiol.64.092501.114547 PMID: 11826273

12. Malenka RC, Bear MF. LTP and LTD: An Embarrassment of Riches. Neuron. 2004 Sep 30; 44(1):5–21. https://doi.org/10.1016/j.neuron.2004.09.012 PMID: 15450156

13. Nevian T, Sakmann B. Spine Ca2+ Signaling in Spike-Timing-Dependent Plasticity. J Neurosci. 2006 Oct 25; 26(43):11001–13. https://doi.org/10.1523/JNEUROSCI.1749-06.2006 PMID: 17065442

14. Graupner M, Brunel N. Mechanisms of induction and maintenance of spike-timing dependent plasticity in biophysical synapse models. Front Comput Neurosci [Internet]. 2010 [cited 2022 May 18]; 4. Available from: https://www.frontiersin.org/article/10.3389/fncom.2010.00136

15. Schulz PE. Long-term potentiation involves increases in the probability of neurotransmitter release. Proc Natl Acad Sci. 1997 May 27; 94(11):5888–93. https://doi.org/10.1073/pnas.94.11.5888 PMID: 9159170

16. Yang Y, Calakos N. Presynaptic long-term plasticity. Front Synaptic Neurosci [Internet]. 2013 [cited 2023 Jan 24];5. Available from: https://www.frontiersin.org/articles/10.3389/fnsyn.2013.00008 PMID: 24146648

17. Abegg MH, Savic N, Ehrengruber MU, McKinney RA, Gåhwiler BH. Epileptiform activity in rat hippocampus strengthens excitatory synapses. J Physiol. 2004; 554(2):439–48. https://doi.org/10.1113/jphysiol.2003.052662 PMID: 14594985

18. Goddard GV, Douglas RM. Does the Engram of Kindling Model the Engram of Normal Long Term Memory? Can J Neurol Sci J Can Sci Neurol. 1975 Nov; 2(4):385–94. https://doi.org/10.1017/s0317167100020539 PMID: 811349

19. Wu K, Castellano D, Tian Q, Lu W. Distinct regulation of tonic GABAergic inhibition by NMDA receptor subtypes. Cell Rep. 2021 Nov 9; 37(6):109960. https://doi.org/10.1016/j.celrep.2021.109960 PMID: 34758303

20. Hardingham GE, Bading H. Synaptic versus extrasynaptic NMDA receptor signalling: implications for neurodegenerative disorders. Nat Rev Neurosci. 2010 Oct; 11(10):682–96. https://doi.org/10.1038/nrn2911 PMID: 20842175

21. Lee HHC, Deeb TZ, Walker JA, Davies PA, Moss SJ. NMDA receptor activity downregulates KCC2 resulting in depolarizing GABAA receptor-mediated currents. Nat Neurosci. 2011 Jun; 14(6):736–43. https://doi.org/10.1038/nn.2806 PMID: 21532577

22. Kahle KT, Deeb TZ, Puskarjov M, Silayeva L, Liang B, Kaila K, et al. Modulation of neuronal activity by phosphorylation of the K–Cl cotransporter KCC2. Trends Neurosci. 2013 Dec 1; 36(12):726–37. https://doi.org/10.1016/j.tins.2013.08.006 PMID: 24139641

23. Kahle KT, Khanna AR, Duan J, Staley KJ, Delpire E, Poduri A. The KCC2 Cotransporter and Human Epilepsy: Getting Excited About Inhibition. The Neuroscientist. 2016 Dec 1; 22(6):555–62. https://doi.org/10.1177/1073858416645087 PMID: 27130838

24. Wu C, Sun D. GABA receptors in brain development, function, and injury. Metab Brain Dis. 2015 Apr; 30(2):367–79. https://doi.org/10.1007/s11011-014-9560-1 PMID: 24820774

25. Chen Q, He S, Hu XL, Yu J, Zhou Y, Zheng J, et al. Differential Roles of NR2A- and NR2B-Containing NMDA Receptors in Activity-Dependent Brain-Derived Neurotrophic Factor Gene Regulation and Limbic Epileptogenesis. J Neurosci. 2007 Jan 17; 27(3):542–52. https://doi.org/10.1523/JNEUROSCI.3607-06.2007 PMID: 17234586

26. Goldensohn ES. The Relevance of Secondary Epileptogenesis to the Treatment of Epilepsy: Kindling and the Mirror Focus. Epilepsia. 1984; 25(s2):S156–68.

27. Shen Y, Gong Y, Ruan Y, Chen Z, Xu C. Secondary Epileptogenesis: Common to See, but Possible to Treat? Front Neurol [Internet]. 2021 [cited 2024 Jan 25]; 12. Available from: https://www.frontiersin.org/articles/10.3389/fneur.2021.747372 PMID: 34938259

28. Kuang Y, Xu C, Zhang Y, Wang Y, Wu X, Wang Y, et al. Low-frequency stimulation of the primary focus retards positive transfer of secondary focus. Sci Rep. 2017 Mar 23; 7(1):345. https://doi.org/10.1038/s41598-017-00479-z PMID: 28336934
29. Shouval HZ, Bear MF, Cooper LN. A unified model of NMDA receptor-dependent bidirectional synaptic plasticity. Proc Natl Acad Sci. 2002 Aug 6; 99(16):10831-6. https://doi.org/10.1073/pnas.152343099 PMID: 12136127

30. Chindemi G, Abdellah M, Amsalem O, Benavides-Piccione R, Delattre V, Doron M, et al. A calcium-based plasticity model for predicting long-term potentiation and depression in the neocortex. Nat Commun. 2022 Jun 1; 13(1):3038. https://doi.org/10.1038/s41467-022-30214-w PMID: 35650191

31. Kumar A, Mehta M. Frequency-Dependent Changes in NMDAR-Dependent Synaptic Plasticity. Front Comput Neurosci [Internet]. 2011 [cited 2022 May 18]; 5. Available from: https://www.frontiersin.org/article/10.3389/fncom.2011.00038 PMID: 21994493

32. Graupner M, Brunel N. Calcium-based plasticity model explains sensitivity of synaptic changes to spike pattern, rate, and dendritic location. Proc Natl Acad Sci. 2012 Mar 6; 109(10):3991-6. https://doi.org/10.1073/pnas.1109359109 PMID: 22357758

33. Cai Y, Gavornik JP, Cooper LN, Yeung LC, Shouval HZ. Effect of Stochastic Synaptic and Dendritic Dynamics on Synaptic Plasticity in Visual Cortex and Hippocampus. J Neurophysiol. 2007 Jan; 97 (1):375-86. https://doi.org/10.1152/jn.00895.2006 PMID: 17035360

34. Inglebert Y, Aljadeff J, Brunel N, Debanne D. Synaptic plasticity rules with physiological calcium levels. Proc Natl Acad Sci. 2020 Dec 29; 117(52):33639-48. https://doi.org/10.1073/pnas.2013663117 PMID: 33328274

35. Shouval HZ, Castellani GC, Blais BS, Yeung LC, Cooper LN. Converging evidence for a simplified bio-physical model of synaptic plasticity. Biol Cybern. 2002 Dec 1; 87(5):383-91. https://doi.org/10.1007/s00422-002-0362-x PMID: 12461628

36. Senn W, Markram H, Tsodyks M. An Algorithm for Modifying Neurotransmitter Release Probability Based on Pre- and Postsynaptic Spike Timing. Neural Comput. 2001 Jan 1; 13(1):35-67. https://doi.org/10.1162/089976601300014628 PMID: 11177427

37. Costa RP, Froemke RC, Sjöström PJ, van Rossum MC. Unified pre- and postsynaptic long-term plasticity enables reliable and flexible learning. Nelson SB, editor. eLife. 2015 Aug 26; 4:e09457. https://doi.org/10.7554/eLife.09457 PMID: 26308579

38. Fung PK, Robinson PA. Neural field theory of calcium dependent plasticity with applications to transepial magnetic stimulation. J Theor Biol. 2013 May 7; 324:72-83. https://doi.org/10.1016/j.jtbi.2013.01.013 PMID: 23376643

39. Modolo J, Thomas A, Legros A. Neural mass modeling of power-line magnetic fields effects on brain activity. Front Comput Neurosci [Internet]. 2013 [cited 2023 Feb 24]; 7. Available from: https://www.frontiersin.org/articles/10.3389/fncom.2013.00034 PMID: 23596412

40. Wendling F, Bartolomei F, Bellanger JJ, Chauvel P. Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition. Eur J Neurosci. 2002; 15(9):1499-508. https://doi.org/10.1046/j.1460-9568.2002.01985.x PMID: 12028360

41. Tsodyks M, Pawelzik K, Markram H. Neural Networks with Dynamic Synapses. Neural Comput. 1998 May; 10(4):821-35. https://doi.org/10.1162/089976698300017502 PMID: 9573407

42. Markram H, Wang Y, Tsodyks M. Differential signaling via the same axon of neocortical pyramidal neurons. Proc Natl Acad Sci. 1998 Apr 28; 95(9):5323-8. https://doi.org/10.1073/pnas.95.9.5323 PMID: 9560274

43. Wang Y, Markram H, Goodman PH, Berger TK, Ma J, Goldman-Rakic PS. Heterogeneity in the pyramidal network of the medial prefrontal cortex. Nat Neurosci. 2006 Apr; 9(4):534-42. https://doi.org/10.1038/nn1670 PMID: 16547512

44. Sjöström PJ, Turrigiano GG, Nelson SB. Rate, Timing, and Cooperativity Jointly Determine Cortical Synaptic Plasticity. Neuron. 2001 Dec; 32(6):1149-64. https://doi.org/10.1016/s0896-6273(01)00542-6 PMID: 11754844

45. Lagarde S, Buzori S, Trebuchon A, Carron R, Scavarda D, Milh M, et al. The repertoire of seizure onset patterns in human focal epilepsies: Determinants and prognostic values. Epilepsia. 2019 Jan; 60(1):85-95. https://doi.org/10.1111/epi.14604 PMID: 30426477

46. Köksal Ersöz E, Modolo J, Bartolomei F, Wendling F. Neural mass modeling of slow-fast dynamics of seizure initiation and abortion. PLOS Comput Biol. 2020 Nov 9; 16(11):e1008430. https://doi.org/10.1371/journal.pcbi.1008430 PMID: 33166277

47. Zito K, Scheuss V. NMDA Receptor Function and Physiological Modulation. In: Squire LR, editor. Encyclopaedia of Neuroscience [Internet]. Oxford: Academic Press; 2009 [cited 2023 Feb 3]. p. 1157-64. Available from: https://www.sciencedirect.com/science/article/pii/B9780080450469012250

48. Müller L, Tokay T, Porath K, Köhling R, Kirschstein T. Enhanced NMDA receptor-dependent LTP in the epileptic CA1 area via upregulation of NR2B. Neurobiol Dis. 2013 Jun 1; 54:183-93.
49. Wilder BJ. The mirror focus and secondary epileptogenesis. In: International Review of Neurobiology [Internet]. Academic Press; 2001 [cited 2024 Jul 2]. p. 435–46. (Brain Plasticity and Epilepsy; vol. 45). Available from: https://www.sciencedirect.com/science/article/pii/S0074774201450227

50. Doedel EJ, Champneys A, Fairgrieve TF, Yu AB, Kuznetsov A.P., Oldeman BE, et al. Auto-07p: Continuation and bifurcation software for ordinary differential equations [Internet]. 2007. Available from: http://cmvl.cs.concordia.ca/auto/

51. Kuehn C. Multiple Time Scale Dynamics [Internet]. Springer International Publishing; 2015 [cited 2020 Apr 9]. (Applied Mathematical Sciences). Available from: https://www.springer.com/fr/book/9783319123158

52. Szente MB, Boda B. Cellular mechanisms of neocortical secondary epileptogenesis. Brain Res. 1994 Jun 20; 648(2):203–14. https://doi.org/10.1016/0006-8993(94)91119-3 PMID: 7922535

53. Gilmore R, Morris H, Van Ness PC, Gilmore-Pollak W, Estes M. Mirror focus: function of seizure frequency and influence on outcome after surgery. Epilepsia. 1994; 35(2):258–63. https://doi.org/10.1111/j.1528-1157.1994.tb02429.x PMID: 8156943

54. Arnth-Jensen N, Jabaudon D, Scanziani M. Cooperation between independent hippocampal synapses is controlled by glutamate uptake. Nat Neurosci. 2002 Apr; 5(4):325–31. https://doi.org/10.1038/nn825 PMID: 11896395

55. Wendling F, Benquet P, Bartolomei F, Jirsa V. Computational models of epileptiform activity. J Neurosci Methods. 2016 Feb; 260:233–51. https://doi.org/10.1016/j.jneumeth.2015.03.027 PMID: 25843066

56. Moran RJ, Stephan KE, Dolan RJ, Friston KJ. Consistent spectral predictors for dynamic causal models of steady-state responses. NeuroImage. 2011 Apr 15; 55(4):1694–708. https://doi.org/10.1016/j.neuroimage.2011.01.012 PMID: 21238593

57. Sheheitli H, Jirsa V. Incorporating slow NMDA-type receptors with nonlinear voltage-dependent magnesium block in a next generation neural mass model: derivation and dynamics [Internet]. bioRxiv; 2023 [cited 2023 Nov 23]. p. 2023.07.03.547465. Available from: https://www.biorxiv.org/content/10.1101/2023.07.03.547465v1

58. Kurbatova P, Wendling F, Kaminska A, Rosati A, Nabbout R, Guerrini R, et al. Dynamic changes of depolarizing GABA in a computational model of epileptogenic brain: Insight for Dravet syndrome. Exp Neurol. 2016 Sep; 283:57–72. https://doi.org/10.1016/j.expneurol.2016.05.037 PMID: 27246997

59. Lopez-Sola E, Sanchez-Todo R, Lleal È, Köksal-Ersöz E, Yochum M, Makhalova J, et al. A personalizable autonomous neural mass model of epileptic seizures. J Neural Eng. 2022 Sep; 19(5):055002. https://doi.org/10.1088/1741-2552/ac8ba8 PMID: 35995031

60. Magloire V, Cornford J, Lieb A, Kullmann DM, Pavlov I. KCC2 overexpression prevents the paradoxical seizure-promoting action of somatic inhibition. Nat Commun. 2019 Mar 15; 10(1):1225. https://doi.org/10.1038/s41467-019-08933-4 PMID: 30874549

61. Marlin JJ, Carter AG. GABA-A Receptor Inhibition of Local Calcium Signaling in Spines and Dendrites. J Neurosci. 2014 Nov 26; 34(48):15898–911. https://doi.org/10.1523/JNEUROSCI.0869-13.2014 PMID: 25429132

62. Urban-Ciecko J, Barth AL. Somatostatin-expressing neurons in cortical networks. Nat Rev Neurosci. 2016 Jul; 17(7):401–9. https://doi.org/10.1038/nrn.2016.53 PMID: 27225074

63. Goldberg JH, Lacefield CO, Yuste R. Global dendritic calcium spikes in mouse layer 5 low threshold spiking interneurones: implications for control of pyramidal cell bursting. J Physiol. 2004 Jul; 558 (2):465–78. https://doi.org/10.1113/jphysiol.2004.064519 PMID: 15146046

64. Mihaly A, Szente M, Dubravcsik Z, Boda B, Kiraly E, Nagy T, et al. Parvalbumin- and calbindin-containing neurons express c-fos protein in primary and secondary (mirror) epileptic foci of the rat neocortex. Brain Res. 1997 Jun 27; 761(1):135–45.

65. Dhaher R, Gruenbaum SE, Sandhu MRS, Ottestad-Hansen S, Tu N, Wang Y, et al. Network-Related Changes in Neurotransmitters and Seizure Propagation During Rodent Epileptogenesis. Neurology [Internet]. 2021 May 4 [cited 2024 Jul 16];96(18). Available from: https://www.neurology.org/doi/ https://doi.org/10.1212/WNL.0000000000011846 PMID: 33722994

66. Desroches M, Faugeras O, Krupa M. Slow-Fast Transitions to Seizure States in the Wendling-Chauvel Neural Mass Model. Opera Medica Physiol. 2015 Dec 1; 2(3–4):228–34.

67. Gnatkovsky V, Librizzi L, Trombin F, Curtis M de. Fast activity at seizure onset is mediated by inhibitory circuits in the entorhinal cortex in vitro. Ann Neurol. 2008; 64(6):674–86.

68. Ben-Ari Y, Crepel V, Represa A. Seizures Beget Seizures in Temporal Lobe Epilepsies: The Boomerang Effects of Newly Formed Aberrant Kainatergic Synapses. Epilepsy Curr. 2008 May 1; 8(3):68–72. https://doi.org/10.1111/j.1535-7511.2008.00241.x PMID: 18488058
69. Ben-Ari Y, Dudek FE. Primary and Secondary Mechanisms of Epileptogenesis in the Temporal Lobe: There is a before and an After. Epilepsy Curr. 2010 Sep 1; 10(5):118–25. https://doi.org/10.1111/j. 1535-7511.2010.01376.x PMID: 20944823

70. Bartolomei F, Lagarde S, Wendling F, McGonigal A, Jirsa V, Guye M, et al. Defining epileptogenic networks: Contribution of SEEG and signal analysis. Epilepsia. 2017; 58(7):1131–47. https://doi.org/10. 1111/epi.13791 PMID: 28543030

71. Lehnertz K, Bröhl T, Wrede R von. Epileptic-network-based prediction and control of seizures in humans. Neurobiol Dis. 2023 Jun 1; 181:106098.

72. Bartolomei F, Cosandier-Rimele D, McGonigal A, Aubert S, Régis J, Gavaret M, et al. From mesial temporal lobe to temporoperisylvian seizures: A quantified study of temporal lobe seizure networks. Epilepsia. 2010; 51(10):2147–58. https://doi.org/10.1111/j.1528-1167.2010.02690.x PMID: 20738379

73. Bartolomei F, Chauvel P, Wendling F. Epileptogenicity of brain structures in human temporal lobe epilepsy: a quantified study from intracerebral EEG. Brain. 2008 Jul 1; 131(7):1818–30.

74. Bjellvi J, Olsson I, Malmgren K, Wilbe Ramsay K. Epilepsy duration and seizure outcome in epilepsy surgery. Neurology. 2019 Jul 9; 93(2):e159–66.

75. Mauritz M, Hirsch LJ, Camfield P, Chin R, Nardone R, Lattanzi S, et al. Acute symptomatic seizures: an educational, evidence-based review. Epileptic Disord. 2022; 24(1):26–49. https://doi.org/10.1684/epd. 2021.1376 PMID: 34789447

76. Alashjaie R, Kerr EN, AlShoumer A, Hawkins C, Yau I, Weiss S, et al. Surgical outcomes in children with drug-resistant epilepsy and hippocampal sclerosis. Epilepsy Res. 2024 Jul 1; 203:107367. https:// doi.org/10.1016/j.eplepsyres.2024.107367 PMID: 38703703

77. Ying Z, Wang I, Blümcke I, Bulacio J, Alexopoulos A, Jehi L, et al. A comprehensive clinico-pathological and genetic evaluation of bottom-of-sulcus focal cortical dysplasia in patients with difficult-to-localize focal epilepsy. Epileptic Disord. 2019; 21(1):65–77. https://doi.org/10.1684/epd.2019.1028 PMID: 30782578

78. Macdonald-Laurs E, Warren AEL, Lee WS, Yang JYM, MacGregor D, Lockhart PJ, et al. Intrinsic and secondary epileptogenicity in focal cortical dysplasia type II. Epilepsia. 2023; 64(2):348–63. https://doi. org/10.1111/epi.17495 PMID: 36527426

79. Halász P, Bódizs R, Ujma PP, Fabó D, Szűcs A. Strong relationship between NREM sleep, epilepsy and plastic functions—A conceptual review on the neurophysiology background. Epilepsy Res. 2019 Feb 1; 150:95–105. https://doi.org/10.1016/j.eplepsyres.2018.11.008 PMID: 30712997

80. Gelinas JN, Khodagholy D, Thesen T, Devinsky O, Buzsáki G. Interictal epileptiform discharges induce hippocampal-cortical coupling in temporal lobe epilepsy. Nat Med. 2016 Jun; 22(6):641–8. https://doi. org/10.1038/nm.4084 PMID: 27111281

81. Dahal P, Ghani N, Flinker A, Dugan P, Friedman D, Doyle W, et al. Interictal epileptiform discharges shape large-scale intercortical communication. Brain. 2019 Nov 1; 142(11):3502–13. https://doi.org/10. 1093/brain/awz269 PMID: 31501850

82. Yu H, Kim W, Park DK, Phi JH, Lim BC, Chae JH, et al. Interaction of interictal epileptiform activity with sleep spindles is associated with cognitive deficits and adverse surgical outcome in pediatric focal epilepsy. Epilepsia [Internet]. 2023 [cited 2023 Nov 24]; n/a(n/a). Available from: https://onlinelibrary.wiley. com/doi/abs/10.1111/epi.17810 PMID: 37983643

83. Khalife MR, Scott RC, Hernan AE. Mechanisms for Cognitive Impairment in Epilepsy: Moving Beyond Seizures. Front Neurol [Internet]. 2022 [cited 2022 Aug 22]; 13. Available from: https://www.frontiersin. org/articles/10.3389/fneur.2022.878991 PMID: 35645970

84. Dallmer-Zerbe I, Jiruska P, Hlinka J. Personalized dynamic network models of the human brain as a future tool for planning and optimizing epilepsy therapy. Epilepsia. 2023; 64(9):2221–38. https://doi.org/ 10.1111/epi.17690 PMID: 37340565

85. Jiruska P, Freestone D, Gnatkovsky V, Wang Y. An update on the seizures beget seizures theory. Epilepsia. 2023;00(n/a):1–12. https://doi.org/10.1111/epi.17721 PMID: 37466948

86. Berg AT, Shinnar S. Do Seizures Beget Seizures? An Assessment of the Clinical Evidence in Humans. J Clin Neurophysiol. 1997 Mar; 14(2):102. https://doi.org/10.1097/00004691-199703000-00003 PMID: 9165405

87. Kullmann DM, Lamsa KP. Long-term synaptic plasticity in hippocampal interneurons. Nat Rev Neurosci. 2007 Sep; 8(9):687–99.

88. Honoré E, Khlaifia A, Bosson A, Lacaille JC. Hippocampal Somatostatin Interneurons, Long-Term Synaptic Plasticity and Memory. Front Neural Circuits [Internet]. 2021 [cited 2023 Jul 28]; 15. Available from: https://www.frontiersin.org/articles/10.3389/fncir.2021.687558 PMID: 34149368
89. Karnani MM, Agetsuma M, Yuste R. A blanket of inhibition: functional inferences from dense inhibitory connectivity. Curr Opin Neurobiol. 2014 Jun 1; 26:96–102. https://doi.org/10.1016/j.conb.2013.12.015 PMID: 24440415

90. Wenzel M, Hamm JP, Peterka DS, Yuste R. Acute focal seizures start as local synchronizations of neuronal ensembles. J Neurosci. 2019 Aug 19;3176–18. https://doi.org/10.1523/JNEUROSCI.3176-18. 2019 PMID: 31427393

91. Magloire V, Mercier MS, Kullmann DM, Pavlov I. GABAergic Interneurons in Seizures: Investigating Causality With Optogenetics. The Neuroscientist. 2019 Aug 1; 25(4):344–58. https://doi.org/10.1177/1073858418805002 PMID: 30317911

92. Zarnadze S, Bäuerle P, Santos-Torres J, Böhm C, Schmitz D, Geiger JR, et al. Cell-specific synaptic plasticity induced by network oscillations. Bartos M, editor. eLife. 2016 May 24; 5:e14912. https://doi.org/10.7554/eLife.14912 PMID: 27218453

93. Hadler MD, Tzilivaki A, Schmitz D, Alle H, Geiger JRP. Gamma-Oscillation Plasticity Is Mediated by Parvalbumin Interneurons [Internet]. bioRxiv; 2023 [cited 2023 Aug 25], p. 2023.06.21.545901. Available from: https://www.biorxiv.org/content/10.1101/2023.06.21.545901v1

94. Litwin-Kumar A, Doiron B. Formation and maintenance of neuronal assemblies through synaptic plasticity. Nat Commun. 2014 Nov 14; 5(1):5319. https://doi.org/10.1038/ncomms6319 PMID: 25395015

95. Issa NP, Nunn KC, Wu S, Haider HA, Tao JX. Putative roles for homeostatic plasticity in epileptogenesis. Epilepsia. 2023; 64(3):539–52. https://doi.org/10.1111/epi.17500 PMID: 36617338

96. Abraham WC, Jones OD, Glanzman DL. Is plasticity of synapses the mechanism of long-term memory storage? Npj Sci Learn. 2019 Jul 2; 4(1):1–10. https://doi.org/10.1038/s41539-019-0048-y PMID: 31285847

97. Caputi A, Melzer S, Michael M, Monyer H. The long and short of GABAergic neurons. Curr Opin Neurobiol. 2013 Apr 1; 23(2):179–86. https://doi.org/10.1016/j.conb.2013.01.021 PMID: 23394773

98. Cho KKA, Shi J, Phensy AJ, Turner ML, Sohal VS. Long-range inhibition synchronizes and updates prefrontal task activity. Nature. 2023 May; 617(7961):548–54. https://doi.org/10.1038/s41586-023-06012-9 PMID: 37100905

99. Bikbaev A, Manahan-Vaughan D. Relationship of hippocampal theta and gamma oscillations to potentiation of synaptic transmission. Front Neurosci [Internet]. 2008 [cited 2023 Mar 2]; 2. Available from: https://www.frontiersin.org/articles/10.3389/neuro.01.010.2008 PMID: 18982107