Adaptation and Inhibition Control Pathological Synchronization in a Model of Focal Epileptic Seizure

Visual Abstract


Introduction
Epilepsy is the fourth most common neurologic disorder, and is responsible for a greater total global burden of disease than any neurologic conditions except for stroke and migraine (Beghi et al., 2005;Rothstein et al., 2005;Chin and Vora, 2014). Epileptic seizures are characterized by the increased excitability/excitation in the brain's recurrently coupled neuronal networks (Lytton, 2008). Typically, experimental seizure models assume that seizures occur due to decreased inhibition (Karnup and Stelzer, 1999;Sivakumaran et al., 2015) or increased excitation in the neural networks (Ursino and la Cara, 2006;Hall and Kuhlmann, 2013).
There is also evidence that interneurons increase their firing at seizure initiation (Lillis et al., 2012) and are active during the time course of the epileptic activity , suggesting that the activity of interneurons contributes importantly to aspects of seizure dynamics. The activity-dependent interplay between the pyramidal cells and interneurons could play an essential role for seizure generation mechanisms (Krishnan and Bazhenov, 2011;Naze et al., 2015;Buchin et al., 2016b). In neural mass models, neuron populations are often treated as rate units lacking intrinsic adaptation (Touboul et al., 2011). The dy-namic behavior of the neural populations is determined by the balance between excitation and inhibition. Despite the simplicity of these models, they can be successfully used to reproduce resting and interictal states as well as ictal discharges by producing time series comparable with macroscopic measurements such as electroencephalogram signals and field potentials (Demont-Guignard et al., 2009).
However, not all types of epileptic seizures can be explained by looking only at the balance between excitation and inhibition (Traub et al., 2005); intrinsic excitability changes on the single-neuron level also play an important role (Krishnan and Bazhenov, 2011). Studies on human subiculum tissue showed that the complete blockade of type A GABAergic neurotransmission (and thus inactivation of the effects of inhibitory population) precludes seizure emergence while, if applied after seizure initiation, it abolishes rather than enhances the seizure activity. These manipulations usually bring back the neural network in the slice toward pre-ictal events, which have substantially different frequency content than seizure activity , and which in this case fail to trigger ictal events. In human epileptic tissues, including peritumoral neocortex (Pallud et al., 2014), interictal discharges are generated spontaneously. These events are triggered by interneurons which depolarize pyramidal cells with impaired chloride regulation, leading to depolarizing effects of GABA. Once activated, pyramidal cells excite other cells via AMPA-mediated glutamatergic transmission. In these tissues, seizures can be produced by increasing local excitability using modified bathing media. The transition to seizures is characterized by the emergence of specific pre-ictal events initiated by pyramidal cells which synchronize local neurons by AMPA synapses. These pre-ictal events cluster before seizure initiation which requires functional AMPA, NMDA as well as GABA A signals. The conventional neural mass models are unable to explain these pre-ictal oscillations because they require the excitatory population to generate periodic oscillations in the absence of inhibition. The second motivation for incorporating intrinsic excitability into neural mass models is that in epileptogenic areas, such as human subiculum, there is a substantial proportion of neurons with non-trivial intrinsic properties such as spike-frequency adaptation (Jensen et al., 1994;Huberfeld et al., 2007). To take these properties into account, neural mass models need to be enriched by the addition of components such as slow potassium currents (Pinsky and Rinzel, 1994).
In addition, seizures are typically accompanied by high potassium concentrations (Dietzel and Heinemann, 1986;Xiong and Stringer, 1999;Fröhlich and Bazhenov, 2006;Florence et al., 2009), which in turn activate calcium currents (Bazhenov et al., 2004;Fröhlich et al., 2008), which in turn affect spike-frequency adaptation and intrinsic bursting. These properties are likely to modulate the single neuron firing and thus further influence the neuronal dynamics. These findings motivate the development of neural mass models that can capture the intrinsic excitability in coupled neural populations.
In this work, we developed a novel neural mass model consisting of an inhibitory neural population and an adaptive excitatory neuronal population (Buchin and Chizhov, 2010b). We calibrated the parameters of the model to local field potential (LFP) data recorded in human subiculum slices during rest, seizure, and full disinhibition in pre-ictal condition. We then analyzed the model as calibrated to each of these three regimes. Our results emphasize the role of intrinsic excitability such as adaptation in the excitatory population, which help explain the transitions between rest, seizure, and full disinhibition states.

Epileptic tissue
Temporal lobe tissue blocks containing the hippocampus, subiculum, and part of the entorhinal cortex were obtained from 45 people of both sexes with pharmacoresistant medial temporal lobe epilepsies associated with hippocampal sclerosis (age, 18 -52 years; seizures for 3-35 years) undergoing resection of the amygdala, the hippocampus, and the anterior parahippocampal gyrus. All of the individuals gave their written informed consent and the study was approved by the Comité Consultatif National d'Ethique.

Tissue preparation
The post-surgical tissue was transported in a cold, oxygenated solution containing 248 mM D-sucrose, 26 mM NaHCO 3 , 1 mM KCl, 1 mM CaCl 2 , 10 mM MgCl 2 , and 10 mM D-glucose, equilibrated with 5% CO 2 in 95% O 2 . Hippocampal-subicular-entorhinal cortical slices or isolated subicular slices (400-m thickness, 3 ϫ 12 mm length and width) were cut with a vibratome (HM650 V, Microm). They were maintained at 37°C and equilibrated with 5% CO 2 in 95% O 2 in an interface chamber perfused with a solution containing 124 mM NaCl, 26 mM NaHCO 3 , 4 mM KCl, 2 mM MgCl 2 , 2 mM CaCl 2 , and 10 mM D-glucose. Bicuculline or picrotoxin was used to block GABA A receptors. Ictal-like activity was induced by increasing the external K ϩ concentration to 8 mM and reducing the Mg 2ϩ concentration to 0.25 mM to increase the cellular excitability (similar to Huberfeld et al., 2011).

Recordings
Up to four tungsten electrodes etched to a tip diameter of ϳ5 m were used for the extracellular recordings. The signals were amplified 1000-fold and filtered to pass frequencies of 0.1 Hz to 10 kHz (AM Systems, 1700). The extracellular signals were digitized at 10 kHz with a 12-bit, 16-channel A-D converter (Digidata 1200A, Molecular Devices) and monitored and saved to a PC with Axoscope (Molecular Devices).

Data analysis
Records were analyzed using pCLAMP 10 software and scripts written in MATLAB 2016a. Power spectrum estimation was performed using fast Fourier transforms. The major frequencies of oscillations were computed via the multitaper method (Thomson, 1982).

Simulations and analysis
Neural population model simulations were performed in XPPAUT 8.0 using the direct Euler method of integration, with a time step of 0.05 ms. Smaller time steps were tested and provided substantially similar results. In all simulations the initial conditions were systematically varied to check stability of numerical results. The data for the model was taken from one representative patient in the brain slice demonstrating resting state, seizure and preictal oscillations.

Neural mass model
In the model we considered interacting excitatory and inhibitory neural populations coupled by AMPA and GABA A synapses. All model parameters and variables are presented in Tables 1, 2. Each population was characterized by the average membrane potential of a population of leaky integrate-and-fire (LIF) neurons (similar to Chizhov and Graham, 2007;Touboul et al., 2011) with approximations for adaptive currents taken from Buchin and Chizhov (2010b): 1.5 mS/cm 2 Excitatory-to-excitatory conductance g EI 1 mS/cm 2 Excitatory-to-inhibitory conductance g IE 2; 0.5; 1 mS/cm 2 Inhibitory-to-excitatory conductance g II 0.2 mS/cm 2 Inhibitory-to-inhibitory conductance U reset E -65 mV Reset membrane potential (Chizhov and Graham, 2007;Buchin and Chizhov, 2010b) V thr E -55 mV Threshold membrane potential (Chizhov and Graham, 2007;Buchin and Chizhov, 2010b)  8.3 ms GABA-A decay time (Chizhov et al., 2002) GABA2 0.2 ms GABA-A rise time (Chizhov, 2002)  In all simulations E/I ͑t͒ has been approximated by the following sigmoid function: The population firing rate determines the adaptive (a), excitatory (e), and inhibitory (i) gating variables. Their dynamics are computed using the second-order approximation Chizhov, 2014): To mimic the afferent excitatory input, the excitatory population also received stochastic excitatory input modeled as an Ornstein-Uhlenbeck process (Buchin and Chizhov, 2010a): To mimic elevated extracellular potassium from epileptogenic slice experiments, in the population model, we increased potassium reversal potential in both populations V K E/I from -90 to -75 mV, i.e., from K o ϭ 4 mM to K 0 ϭ 8 mM. This value of V K E/I was computed based RT / F ϭ 26.64 mV and K i ϭ 138 mM (Krishnan and Bazhenov, 2011). All model parameter values and variable names are present in Table 1, 2. The initial parameter set was chosen manually to reproduce the pre-ictal like oscillations due to balance between g EE and g AHP , seizure and resting state were fit such that g EI parameter variations would make a transition between seizure and resting state.

LFP model
The LFP was calculated based on the activity of the excitatory population. We assumed that pyramidal cells activity dominates the extracellular field (Buzsáki et al., 2012). The dominant theory is that the LFP component is dominated by the single neuron dipole contribution (Buzsáki et al., 2012). Since the neural mass model averages over single neurons, the dipole moment cannot be directly modeled. Thus, to approximate the LFP being recorded near somas of the excitatory populations, we used the assumption that the average membrane potential of the excitatory population is proportional to the LFP, i.e., LFP ϰ U E (Ursino and la Cara, 2006; Demont-Guignard et al., 2009;Ratnadurai-Giridharan et al., 2014).

Construction of the population model
We developed model of interacting excitatory and inhibitory population inspired by Wilson-Cowan approach (Wilson and Cowan, 1972), which consists of excitatory and inhibitory populations coupled by synaptic connections (Fig. 1A). The firing rate in each population depends on the average membrane potential U E/I , which is governed by the subthreshold dynamics of LIF neuron population similar to (Gerstner and Kistler, 2002;Chizhov, 2014; as explained in Materials and Methods). Firing rates of the excitatory and inhibitory populations are determined using the values of U E/I put through function A ͑U E/I ͒ (Johannesma, 1968;Gerstner and Kistler, 2002). To make the model numerically stable and amenable to bifurcation analysis we used a sigmoid function to estimate the population firing rate provided by the A͑U E/I ͒ approximation. To justify the choice of sigmoid parameters, we used least-squares to match it with the analytical solution (Johannesma, 1968;Fig. 1C,D). The sigmoid approximation allows one to efficiently take into account zero and linear parts of the potential-to-rate transfer functions v E/I ͑t͒, and provides saturation due to the single neuron refractory period (Renart et al., 2004). The sigmoid functions of excitatory and inhibitory populations are shown in Figure 1C,D. The difference between the excitatory and inhibitory populations was taken into account by adjusting passive conductances for sodium, potassium, and chloride leak currents estimated in Krishnan and Bazhenov, (2011) based on dynamic ion concentration model.
The subthreshold U E/I dynamics determine the synaptic g EE e͑ E ͒, g EI e͑ E/I ͒, g II i͑ I ͒, g IE i͑ I ͒, and intrinsic g AHP e͑ E/I ͒ conductances (Fig. 1A), computed according to the population firing-rates E/I . Similar to spiking neural network models (Bazhenov et al., 2004;Ratnadurai-Giridharan et al., 2014), adaption in our population model reduces neural firing in the excitatory population after periods of activity. Excitatory population receives external random synaptic input to model excitation from the rest of the brain similar to (Jansen and Rit, 1995;Touboul et al., 2011). To mimic the experimental epileptogenic conditions of human subiculum slice experiments, the potassium reversal potential was elevated from -95 to -75 mV both in the excitatory and inhibitory populations to provide excitatory drive to reproduce the experimental conditions. Elevation of extracellular potassium also leads the increase of intracellular chloride reducing the efficiency of inhibition due to elevated GABA A reversal potential (Huberfeld et al., 2007;Buchin et al., 2016b). To generate the model output comparable with experimental data, we computed the LFP generated by the excitatory population (Buzsáki et al., 2012). This approximation assumes that all pyramidal cells in the excitatory population contribute equally to the recorded LFP signal (Fig. 1B). Thus, the total LFP near somas depends on the average value of the membrane poten-tial in the excitatory population with a certain dimensionality constant, i.e., LFP ϰ kU E/I .

Reproduction of epileptic oscillations
When the excitatory and inhibitory synaptic currents were dynamically balanced, the activity stayed in the low-firing regime, as indicated by LFP power spectrum (Fig. 2). The recorded pyramidal cell during this period demonstrated sparse firing activity, partially time-locked with the discharges on the LFP. We call this activity in the model the balanced or resting state ( Fig. 2A). In this regime, the model does not generate epileptic oscillations. To evaluate the model performance in this resting state, we compared the synthetic LFP with the experimental LFP recorded between seizures ( Fig. 2A). Similar to the experimental data, we found that in the resting state, the model generates broadband oscillations, with the highest power in the 1-to 15-Hz frequency band. In this regime, the average membrane potential of the excitatory population U E/I stays in the range from -60 to -50 mV.
We found that the model was not capable of generating interictal discharges using this parameter set. It has been recently suggested that interneurons play the key role in generating interictal activity (Cohen et al., 2002;Huberfeld et al., 2011). In the presence of GABA A blockade these events were completely blocked, indicating that they depend on combination of GABAergic and glutamatergic signaling. In the recent population model (Chizhov et al., 2017), it was proposed that interictal discharges could be initiated by the inhibitory population, thus explaining in-terneuron firing before pyramidal cell firing . In our model we have not explored this scenario, i.e., when the inhibitory population is also receiving the background synaptic input. These mechanisms would likely play an important role for seizure initiation; however, incorporating all mechanisms at once would make the model impossible to study analytically. Therefore, we have not considered interictal discharges before seizure, while aiming to specifically describe other types of oscillations.
To reproduce the seizure state in the model, we reduced the synaptic inhibition of the excitatory population by decreasing the synaptic conductance parameter g IE (Fig. 2B, black arrow). All other parameters of the model remained the same. In this case the model moved into an oscillatory regime in which the power spectrum of the oscillations changed dramatically to include strong oscillations in the 1-to 4-Hz frequency band, which is typical for ictal discharges . We compared the model power spectrum with the measured LFP recorded during the initial phase of the ictal discharge with the hypersynchronous activity onset. During this activity regime the recorded pyramidal cells generated strong bursts of spikes temporally locked to the LFP (Fig. 2B). The population model displayed discharges with the same frequency band as in the LFP, indicating large amount of synchrony in the excitatory population (Buzsáki et al., 2012). Note that we considered only the initial phase of the seizure (the whole ictal event is shown in Fig. 3E). g EE , g EI : excitatory to excitatory and excitatory to inhibitory maximal conductances; g II , g IE : inhibitory-to-inhibitory and inhibitory-toexcitatory maximal conductance; g AHP : adaptation conductance in the excitatory population;I E ͑t͒: synaptic noise input to the excitatory population; AHP, afterhyperpolarization current (Buchin and Chizhov, 2010b). B, LFP model: i : contribution of a single excitatory cell; N: the number of neurons; U E : the average membrane potential in the excitatory population. C, D, Sigmoid approximation of potential-to-rate function (Johannesma, 1968) of the excitatory (C) and inhibitory population (D).
To further test the validity of our model, we explored its dynamics with inhibitory activity completely blocked (Fig.  2C). In these simulations the initial conditions were set to the resting state and parameter values of the model were set to the seizure state, but with the conductance g IE (from the inhibitory to the excitatory population) set to zero to mimic the experimental conditions. In this case the GABAergic effects of the inhibitory population in the slice has been fully blocked by bicuculine after seizures have been previously established . In response to this change, the activity in the slice became highly synchronized and reduced to regular pre-ictal discharges. During these oscillations the pyramidal cells generated large bursts of activity, temporally coupled with the LFP (Fig. 2C). In the model, similarly to the experimental preparation, the blockade of the GABAergic signaling mimicked by the abolition of the inhibitory population led to the development of a slow oscillatory rhythm with a peak frequency around 1 Hz. These events have been previously reported as pre-ictal discharges . This rhythm has much slower frequency than seizures, and is usually within the 1-to 4-Hz frequency range Buchin et al., 2016b). In addition, these events recur regularly for long periods with very limited modulation.
We call this regime of activity pre-ictal discharges because similar activity takes place before transition toward an ictal state . In this regime, the dynamics of the excitatory population are determined only by the balance between self-excitation, g EE e͑ E/I ͒, afterhyperpolarization current (AHP; Chizhov and Graham, 2008;Buchin and Chizhov, 2010b), g AHP e͑ E/I ͒, and the afferent LFP is present together with intracellular recording from the pyramidal cell. Each plot contains the model scheme, power spectrum, and time traces provided by the excitatory population U E as well as experimental LFP. Red traces correspond to the model, blue traces to the experiment, and green traces to the intracellular recordings from the pyramidal cells. Model parameters for (A): g EE ϭ 1.5 mS/cm 2 ; g EI ϭ 1 mS/cm 2 ; g IE ϭ 2 mS/cm 2 ; g II ϭ 0.2 mS/cm 2 ; g AHP ϭ 1.6 mS/cm 2 ; (B): g EE ϭ 1.5 mS/cm 2 ; g EI ϭ 1 mS/cm 2 ; g IE ϭ 0.5 mS/cm 2 ; g II ϭ 0.2 mS/cm 2 ; g AHP ϭ 1.6 mS/cm 2 ; (C): g EE ϭ 1.5 mS/cm 2 ; g EI ϭ 1 mS/cm 2 ; g IE ϭ 0 mS/cm 2 ; g II ϭ 0.2 mS/cm 2 ; g AHP ϭ 1.6 mS/cm 2 . synaptic current I E ͑t͒. Hence, these pre-ictal oscillations in the model are driven by the synaptic noise and adaptation. The excitatory input to the excitatory population I E ͑t͒ drives the upswings of U E due to recurrent excitatory synapses, with activity then being terminated by AHP currents. These transitions take place randomly due to stochastic nature of the synaptic input.
For quantitative comparisons between the model and experiment we used the linear fit to the power spectrum over frequencies and peak estimation (Table 3). We found that there is substantial intersection between linear fits applied to the power spectrums in resting, seizure, and pre-ictal states (Fig. 2). We found that there is a substantial overlap between these frequencies, providing validation for the model. Note that we compared the overall spectral characteristics between the model and experiment by variation of only one parameter, g IE to reproduce transitions between the pre-ictal, resting and seizure states. If more parameters are varied at the same time, it would be possible to get a better match between the model and experiment.
Overall oscillations in our population model are controlled by the balance between synaptic currents, adaptation and external synaptic input. When synaptic and  Fig. 2A, g IE ϭ 2 mS/cm 2 ), red corresponds to early seizure (Fig. 2B, g IE ϭ 0.5 mS/cm 2 ), yellow corresponds to late seizure (g IE ϭ 0.25 mS/cm 2 ), and purple corresponds to the disinhibition state (Fig. 2C, intrinsic conductances are balanced, the population demonstrates resting state activity, characterized by a flat power spectrum. When there is an imbalance between excitation and inhibition, populations start developing oscillatory rhythms associated with ictal discharges with a frequency of 3-4 Hz. However, complete loss of inhibition leads to the development of another population rhythm, pre-ictal discharges with 1-Hz frequency, controlled by adaptation and recurrent excitation. Thus, the dynamic state of a neural population depends on the interplay between the intrinsic and synaptic excitability within populations as well as external synaptic input.

Analysis of the population model
To delineate the mechanisms giving rise to the different oscillatory modes in the model, we used continuation techniques and bifurcation analysis. Since it is impossible to use the standard techniques to identify bifurcations in the presence of noise, we analyzed the model in the absence of an external input I E ͑t͒. This allowed us to compute the model behavior in the stationary regime and characterize bifurcations happening during transitions between different oscillatory regimes. The initial parameters were chosen to correspond to the resting state. The parameter variations were calculated around this point in the parameter space for g EE , g EI , g IE , and g II bifurcation diagrams, with other parameters held fixed. Analysis of g AHP and V GABA variations was implemented for another parameter set, where g IE ϭ 0.5 S/cm 2 and g IE ϭ 1 S/cm 2 ; other parameters remained the same.
The frequency of seizure oscillations depends on the strength of the synaptic currents in the population model. There is a nonlinear relationship between seizure major frequency and the recurrent excitatory conductance g EE (Fig. 3A). When the g EE is increased up to 2.8 mS/cm 2 , the model responds with an oscillatory frequency near 7.5 Hz. When self-excitation is further increased up to 4 mS/cm 2 , seizure-related oscillations disappear since the system moves to the high activity state due to sigmoidal saturation of the transfer function (Fig. 1C,D). The amount of stimulation of the inhibitory population also influences the oscillatory frequency. When g EI is in the range of 0 to 0.29 mS/cm 2 (Fig. 3B), the population model generates seizure activity with frequencies of 1.2-2.5 Hz. Note that seizure oscillations are possible even when g EI ϭ 0 mS/cm 2 .
Inhibitory synaptic connections also affect the oscillatory frequency of seizure activity. When g IE is as low as ϳ0.6 mS/cm 2 (Fig. 3C), the seizure activity starts around 3 Hz; it decreases to ϳ1 Hz when g EI is close to zero (when g IE ϭ 0 mS/cm 2 , there is no seizure activity in the model). The amount of recurrent inhibition also determines the seizure oscillatory frequency (Fig. 3D). Seizure activity can be initiated by sufficient self-inhibition, i.e., when g II is near 2 mS/cm 2 , seizures of 2.5 Hz are observed. When g II increases, the seizure frequency decreases; for example, at 10 mS/cm 2 , seizure activity is ϳ1.8 Hz.
In the previous sections, the population model was calibrated to data for short periods of seizure activity, where the frequency was not substantially changing (Fig.  2B). Yet, one can see that in the experiment, seizure activity is not stationary and its frequency changes over time. The time course of a typical seizure is shown in Figure 3E. Before the seizure starts there is a resting state, characterized by occasional interictal (Cohen et al., 2002) and pre-ictal discharges . When seizure starts at 22 s, it is characterized by fast oscillations of the extracellular field in the range of 5-6 Hz in the initial phase. During the time course of seizure activity, it gradually decreases to 1-Hz frequency, and from 52 s, it gradually stops.
We aimed to reproduce this aspect of seizure activity using the population model (Fig. 3E). First, the model was initialized in the resting state ( Fig. 2A), green trace. Second, we reduced the amount of inhibitory-to-excitatory coupling (to g IE ϭ 0.5 mS/cm 2 ) to reproduce the seizure state, red trace. Third, we gradually reduced the coupling parameter (to g IE ϭ 0.25 mS/cm 2 ) to reduce the oscillation frequency, yellow trace. Fourth, to model the slow oscillations in the end of seizure, we set the coupling parameter to zero (g IE ϭ 0 mS/cm 2 ), violet trace. Finally, we restored it to the original value to bring the model back to rest (g IE ϭ 2 mS/cm 2 ), green trace. This example illustrates how the transition toward seizure in the population model can be achieved by varying only one parameter, the inhibitory-to-excitatory conductance g IE .
To study the amplitude of pathologic oscillations, we performed a bifurcation analysis and tracked changes of the average membrane potential in the excitatory population, U E (Fig. 4), the self-excitation conductance g EE (Fig.  4A). We found that increasing g EE leads to the development of ictal oscillations when its value increases beyond ϳ2.8 mS/cm 2 . During the gradual increase of g EE , the constant steady state loses stability via the supercritical Hopf bifurcation (Izhikevich, 2007), red dot. After passing this point the neural populations start developing seizure oscillations. This activity regime is stable for large g EE variations, implying that seizure dynamics are possible for a large range of recurrent excitation. When g EE becomes higher than a critical value (Ͼ4.1 mS/cm 2 ) and the system loses stability via the subcritical Hopf bifurcation, green dot. It corresponds to the high activity state with no oscillations. This happens due to the sigmoid approximation of the population rate (Johannesma, 1968), when E/I reaches the saturation level (Fig. 1C,D). Second, we considered the excitatory to inhibitory conductance g EI (Fig. 4B). In this case, seizure activity is blocked when g EI is larger than 0.3 mS/cm 2 . If g EI is smaller than 0.3 mS/cm 2 , it leads to seizure activity via a subcritical Hopf bifurcation, green dot. Similar to the g EE bifurcation diagram, seizure dynamics are possible for a large range of g EI . These results show that a decrease in the excitatory conductance from excitatory to inhibitory populations is sufficient to provoke seizure activity. Note that even if g EI ϭ 0 mS/cm 2 , the excitatory population still receives the input from the inhibitory one because potassium reversal potential is elevated. These changes in potassium reversal potential drive both excitatory and inhibitory population even if synaptic drive is not present. For example, when g EI ϭ 0 mS/cm 2 , the increased potassium reversal potential still drives the inhibitory population, providing the inhibitory input to the excitatory population. It happens because it decreases the leak current thus depolarizing the membrane potential of excitatory and inhibitory neurons. Therefore, seizure oscillations are still present because inhibition is still present.
Seizure frequency in this case is near 1.25 Hz (Fig. 3B) and U E oscillates between -61 and -25 mV.
Third, we considered inhibitory-to-excitatory conductance g IE (Fig. 4C). When g IE ϭ 0, the model shows resting state activity. This corresponds to the condition when the inhibitory population does not have any influence on the excitatory one. Experimentally this scenario is achieved when inhibitory neurotransmission is completely blocked. Therefore, in the complete absence of inhibition, seizure activity could not be generated. In turn, pre-ictal oscillations are not possible without the contribution of the external synaptic noise I E ͑t͒ when g IE ϭ 0 mS/cm 2 . When there is stochastic synaptic input, it occasionally brings the system to the oscillatory regime associated with seizures (Fig. 2C). Then oscillations are promoted due to recurrent excitation and terminated via AHP adaptation mechanism. Thus, without synaptic stimulatin of the inhibitory population, the model is incapable of seizure generation. In turn, pre-ictal oscillations do not require inhibition, but strongly depend on the recurrent excitatory-toexcitatory connections g EE , adaptation g AHP , and the external Seizure Rest g EE (mS/cm 2 ) g EI (mS/cm 2 ) g IE (mS/cm 2 ) g II (mS/cm 2 ) g AHP (mS/cm 2 ) V GABA (mV) Rest Rest Rest Rest Figure 4. Analysis of the population model. A-D, Bifurcation diagrams for the variations of the maximal synaptic conductances, including recurrent excitation g EE , excitation from excitatory to inhibitory population g EI , inhibition from inhibitory to excitatory population g IE , and the recurrent inhibition in the inhibitory population g II , respectively. E, F, Bifurcation diagrams for adaptation in the excitatory population g AHP and GABA reversal potential V GABA from the inhibitory-to-excitatory current, g IE i͑U E Ϫ V GABA ͒. Diagrams A-D were calculated for g IE ϭ 2 mS/cm 2 ; E, g IE ϭ 0.5 mS/cm 2 ; and F, g IE 1 mS/cm 2 . The value of U E characterizes the average membrane potential in the resting state and maximal/minimal values of U E during the oscillations. Red and green dots correspond to the supercritical and subcritical Andronov-Hopf bifurcations. Solid and dotted lines depict the stable and unstable solutions.
synaptic input I E ͑t͒. When inhibitory to excitatory conductance g IE becomes strong enough, around 0.65 mS/cm 2 , seizure oscillations become truncated and the system moves back to the resting state via subcritical Andronov-Hopf bifurcation, green dot. Fourth, we evaluated the role of recurrent inhibitory conductance g II for seizure dynamics (Fig. 4D). When there is substantial amount of self-inhibition in the inhibitory population, it leads to an increase of excitation in the whole system because of synaptic coupling. If g II is above 2.1 mS/cm 2 , it leads to the development of seizure oscillations via a supercritical Hopf bifurcation, red dot. Seizure activity in this case persists for the large variations in g II variations, from 2.2 to Ͼ10 mS/cm 2 .
We then analyzed the effect of adaptation in the excitatory population. We found the regime in the parameter space of the model for which g AHP becomes the critical parameter for seizure oscillations. To find this regime we slightly modified the parameter set, where g IE ϭ 0.5 mS/ cm 2 instead of 2 mS/cm 2 . In this case, g AHP could substantially affect seizure oscillations. When g AHP is in the range of 1-3 mS/cm 2 , there is a large region in the parameter space that produces seizure oscillations. If g AHP is larger than 3 mS/cm 2 , the seizure dynamics becomes truncated due to the inhibitory effect of adaptation via subcritical Andronov-Hopf bifurcation, green dot. Yet when adaptation is not strong enough, g AHP is lower, the model demonstrates seizure oscillations. If g AHP is lower than 1 mS/cm 2 , seizure oscillations become impossible and the model moves to the high activity state without oscillations via supercritical Andronov-Hopf bifurcation. Additionally, we found that in the complete absence of adaptation, seizure oscillations are still possible in the model (results not shown), but pre-ictal oscillations could not be generated because of the absence of negative feedback provided by adaptation. To be able to reproduce seizure oscillations together with pre-ictal oscillations induced by GABA A blockade, adaptation in the excitatory population is required.
We further studied the critical role of V GABA for seizure generation. It has been recently found that changes in V GABA are associated with the rhythm generation in the hippocampus (Cohen et al., 2002;Huberfeld et al., 2007). The analysis was performed for slightly modified parameter set, where g IE ϭ1 mS/cm 2 , such that V GABA becomes the bifurcation parameter. The other parameters remained the same. We have changed the initial parameter set to find the region of the parameter space where V GABA could play the crucial role for oscillations. We found that when V GABA is higher than -59 mV, it leads to ictal oscillations (Fig. 4F). If V GABA drops below -48 mV, the oscillations stop. These transitions take place due to supercritical and subscritical Andronov-Hopf bifurcations. Thus, there is substantial range of V GABA where its increase leads to the development of seizures, which might take place due to chloride accumulation before an ictal discharge (Huberfeld et al., 2007;Lillis et al., 2012).
In summary, using bifurcation analysis, we characterized the parameter regions of the model where seizure oscillations could take place. We found that transitions from seizure to rest and from rest to seizure take place via supercritical and subcritical Andronov-Hopf bifurcations. In all studied cases we found that resting and oscillatory solutions exist for large parameter variations, implying the stability of found solutions (Prinz et al., 2004;Marder and Taylor, 2011). We showed that variations of synaptic g EE , g EI , g IE , g II , and intrinsic conductances g AHP could bring the system toward seizure and move it back to the resting state. It implies that combination of recurrent synaptic currents and spike-frequency adaptation in the excitatory population accounts for the transitions between seizure and resting states.

Discussion
The objective of this study was to investigate the role of intrinsic excitability and inhibition as mechanisms of seizure dynamics. We constructed a novel neural mass model, consisting of interacting excitatory and inhibitory neural populations driven by external synaptic input. By comparing the model with the LFP data from human hippocampal/subicular slices, we found that it could accurately represent resting states, ictal discharges, and pre-ictal oscillations after the blockade of inhibition . Analysis of the model showed that synaptic and intrinsic conductances play the most crucial role for transitions between resting and seizure activity. By analyzing the parameter space of the model, we found the oscillatory regimes specific for the resting state and seizure dynamics, and found that transitions between these regimes take place via subcritical and supercritical Andronov-Hopf bifurcations.
Starting with the pioneering work of Wilson and Cowan (Wilson and Cowan, 1972), neural mass models have traditionally aimed to reduce the complexity of neural dynamics toward interactions between excitation and inhibition. This approach has been validated in multiple studies for describing the large-scale brain activity patterns (Jirsa et al., 2010). Additionally, it has been shown that intrinsic properties of single neurons such as spikefrequency adaptation (Fröhlich et al., 2008) substantially change spiking patterns and thus neural dynamics (Kager et al., 2000;Bazhenov et al., 2004;Buchin et al., 2016b). So far these types of interactions have not been explicitly taken into account in neural mass models.
In this work, we developed a novel mass model by adding AHP currents (Buchin and Chizhov, 2010a) to the excitatory population. This allowed to efficiently take into account not only seizure and resting state dynamics  but also pre-ictal oscillations. In our model seizure activity takes place due to imbalance between self-excitation, adaptation and inhibition. We found that reducing the amount of inhibition to the excitatory population provokes seizure activity. Nonetheless, inhibition plays an important role in orchestrating seizures as well (Fig. 2B). We found that the complete lack of inhibition leads to the development of slow oscillations with significantly different frequency content than seizures (Fig. 2C). Thus, we propose that inhibition, together with single neuron intrinsic properties provided by adaptation, plays an important role controlling the seizure dynamics. only in the absence of inhibition (Fig. 2C), as in the experimental data when the GABA A synaptic activity is completely blocked.
The primary advantage of our model compared to more abstract ones such as  is that it provides more firm biophysical explanations linking single neuron properties to population dynamics (Johannesma, 1968;Gerstner and Kistler, 2002;Chizhov and Graham, 2007). Our approach could be extended to take into account the shunting effect of inhibition by adjusting the firing rate transfer function (Chizhov et al., 2014). To describe the additional mechanisms of seizure transition, the present model could include slow activity-dependent parameter changes similar to Ullah et al., 2009;Proix et al., 2014;Chizhov et al., 2018). There are multiple biophysical mechanisms that could play the role of slow variable bringing the network toward seizure (Naze et al., 2015), including dynamic ion concentration of extracellular potassium (Bazhenov et al., 2004), intracellular chloride (Jedlicka et al., 2011;Buchin et al., 2016b), and intracellular sodium (Krishnan and Bazhenov, 2011;Karus et al., 2015), in pyramidal cells. The population model could be further modified to incorporate these slow mechanisms to describe seizure initiation.
A common problem with neural mass models in general is their limited ability to generate the experimentally measurable signals (Lytton, 2008). In this work, we used the average voltage of the excitatory neural population as the approximation of the LFP signal near the neurons' somas (Ratnadurai-Giridharan et al., 2014). Given the distant dependence of the LFP signal, the current model should be considered only as a first approximation (Buzsáki et al., 2012). More detailed approaches describing populations of two-compartmental neurons (Chizhov, 2014;Chizhov et al., 2015) could also provide better approximation for the LFP.
Epilepsy is a complex phenomenon involving the dynamic interactions between multiple components of the nervous system (Lytton, 2008). In this work, we have investigated the particular role of inhibition and adaptation and their implications for seizure dynamics. Reconciling modeling results with experimental data, we have shown that seizure activity cannot be generated in the complete absence of the inhibitory population and adaption in the excitatory population. Further development of theoretical and experimental approaches in epilepsy research may lead to a better understanding of its mechanisms and the development of new therapeutic targets.