Ictal wavefront propagation in slices and simulations with conductance-based refractory density model

The mechanisms determining ictal discharge (ID) propagation are still not clear. In the present study, we aimed to examine these mechanisms in animal and mathematical models of epileptiform activity. Using double-patch and extracellular potassium ion concentration recordings in rat hippocampal-cortical slices, we observed that IDs moved at a speed of about 1 mm/s or less. The mechanisms of such slow propagation have been studied with a mathematical, conductance-based refractory density (CBRD) model that describes the GABA- and glutamatergic neuronal populations’ interactions and ion dynamics in brain tissue. The modeling study reveals two main factors triggerring IDs: (i) increased interneuronal activity leading to chloride ion accumulation and a consequent depolarizing GABAergic effect and (ii) the elevation of extracellular potassium ion concentration. The local synaptic transmission followed by local potassium ion extrusion and GABA receptor-mediated chloride ion accumulation underlies the ID wavefront’s propagation. In contrast, potassium ion diffusion in the extracellular space is slower and does not affect ID’s speed. The short discharges, constituting the ID, propagate much faster than the ID front. The accumulation of sodium ions inside neurons due to their hyperactivity and glutamatergic currents boosts the Na+/K+ pump, which terminates the ID. Knowledge of the mechanism of ID generation and propagation contributes to the development of new treatments against epilepsy.


Introduction
The exact nature of spatiotemporal dynamics of brain activity during seizures is of great importance. However, the mechanisms of epileptic discharge propagation are not yet well understood. Epileptic activity in the brain is divided into several forms, including long-lasting (with a duration of tens of seconds) ictal discharges (IDs) and different types of short-lasting discharges (SDs) of less than one or a few seconds [1]. Several forms of SDs are distinguished: (i) predominantly, GABAergic interictal discharges (IIDs) emerging between IDs; (ii) IID-like preictal discharges (PIDs) that may take place just before an ID; and (iii) GABA-glutamatergic late short discharges (LSDs) that constitute either an ID in its late phase (known as intraictal bursts in [2]) or continuously repeating in status epilepticus. These discharges propagate through the brain tissue at different speeds. The typical speed of SD propagation is about tens of millimeters per second [2][3][4][5][6]. ID propagation is slower; the rate is less than 1 mm/s in the cerebral cortex of patients [7][8][9] and animals in vivo [10,11] and in vitro [2,3,[12][13][14]. The difference in speeds of IDs and SDs suggests different mechanisms of propagation. The consensus among previous studies is that SD generation and propagation mechanisms are not related to ion dynamics, which implies only neuronal excitation and synaptic interactions [6,[15][16][17]. In contrast, the generation of IDs critically depends on ion dynamics [18][19][20][21], which is often mistakenly omitted in many other modeling studies. Whether the contribution of the ion dynamics into the ID propagation is direct-due to ion diffusion-or indirect-as a consequence of neuronal excitation-is still debated [22].
A slowly propagating wave of ID is accompanied by a wave of increasing extracellular potassium ion concentration [1]. However, it is still unclear whether ictal activity causes an increase in potassium ion concentration or, conversely, diffusion of potassium ions from sites of high concentration causes ictal activity. Therefore, at least two alternative mechanisms of ID propagation have been suggested: (i) synaptic transmission between neurons, including spike propagation through axons and passive conduction of postsynaptic signals through dendritic neuronal branches [23,24]; (ii) the diffusion of potassium ions, resulting from substantial transient increases in extracellular potassium ion concentrations that are observed during IDs [1,25,26].
Recently, using the proposed spatially distributed extension of the simple model of IDs and SDs Epileptor-2 [27], we compared the impacts of synaptic transmission and potassium ion diffusion in 2-D simulations [28]. We showed that potassium ion diffusion can contribute to ictal discharge propagation; however, it leads to an ID propagation speed of much less than 0.1 mm/s. Instead, the synaptic connectivity-based mechanism provides a speed comparable to experimental data. However, due to the reductions made in the Epileptor-2, it cannot be compared to electrophysiological experiments in detail and, in particular, does not exclude any indirect effects of potassium ion diffusion. Therefore, here we use our detailed mathematical model based on the conductance-based refractory density (CBRD) approach and extend it to the 1-D spatially distributed case. We simulate the propagation of repeating IDs and compare the simulations with our recordings of neuronal electrical activity and extracellular potassium ion concentration in two points of the entorhinal cortex (EC) at a distance of 2 mm in rat brain slices in the 4-aminopyridine model. The obtained results indicate that the accumulation of extracellular potassium ions due to synaptic activity but not potassium ion diffusion promotes ID propagation.

Experimental observations of ictal front propagation
We induced the epileptiform activity in rat brain slices and performed the simultaneous recordings of the synaptic currents and the [K + ] o fluctuations (Fig 1A). The epileptiform activity in the implemented in vitro model was described in detail in our previous works [21,29,30]. Briefly, two main types of epileptiform discharges were observed: GABA-mediated IIDs (marked with blue arrows in Fig 1A) and subsequent IDs (marked with brown arrows in Fig 1A), which are mediated by both GABA and glutamate. IIDs had a duration of about 1 s and, in some cases, transitioned to the IDs, which had the duration of 25-100 s and emerged every 5-10 min. The IIDs, which directly preceded IDs, will be referred to as preictal discharges (PIDs).
During each ID, we observed a substantial increase in extracellular potassium ion concentration ( Fig 1B). In most cases, the rise in potassium ion concentration came in two stages ( Fig  1C). The first minor rise of [K + ] o corresponded to the GABA-mediated PID. The second more significant rise of [K + ] o occurred during the tonic phase of ID and coincided with the emergence of glutamate-mediated currents.
[K + ] o peaked during the tonic phase of the ID (the average peak [K + ] o was 9.5 ± 0.3 mM, n = 27 IDs from nine slices), after which it gradually decreased. The late-stage clonic discharges produced transient increases of [K + ] o (Fig 1D), which did not change the overall decreasing trend. In half of the IDs included in the analysis, the K + concentration transiently decreased below the level observed before the ID.
Next, we investigated the ID propagation in the neocortex performing simultaneous wholecell patch-clamp recordings in the entorhinal (ERC) and perirhinal cortices (PRC) at a distance of 2 mm between the recording electrodes. In 85% of cases included in the analysis (with the delays longer than 1 s), we detected the IDs in both cortical areas with a delay in the initiation of several seconds (Fig 2), thus obtaining a wide variety of speed values of the order of magnitude of tenths of a millimeter per second. The delay between the IDs was observed in both voltage-clamp and current-clamp modes (Fig 2A and 2B, respectively). Some IDs (15%) failed to propagate and could be detected in only one cortical area ( Fig 2C).
As the shape of the K + transients during each ID is relatively standard with its rising phase being related to the recruitment of glutamatergic neurons into the activity, we measured the time delay between rises of K + transients in the ERC and PRC to estimate the propagation speed (Fig 3). The half-maximal [K + ] o was used as a reference point to measure the delay between the events (Fig 3B), thus ignoring the initial weak component, presumably determined by the potassium ion outflux from synapses activated by long-range connections. The delays between K + transients were comparable to those observed using the patch-clamp method, and the median speed was estimated to be 0.52 mm/s (IQR: 0.34-0.81, n = 17 pairs of IDs from eight slices). No significant difference between propagation speed from ERC to PRC and from PRC to ERC was detected (p = 0.42, Mann-Whitney U test). The recordings, where the ID-related K + transient failed to propagate from one cortical area to another, were not used for speed estimations (Fig 3C). Concluding the presentation of experimental results, we The representative recordings of synaptic currents in the voltage-clamp mode (black trace) and extracellular potassium ion concentration (green trace). V hold = -27 mV, which is between the reversal potentials of GABA A and glutamate receptor-mediated currents. Outward GABA-mediated synaptic current corresponds to IIDs and PIDs and is marked with blue arrows. Inward, predominantly glutamate-mediated currents correspond to IDs and are marked with brown arrows. Note that all IIDs, PIDs, and IDs induce K + transients. The blue bar marks a fragment of the recording, which contains PID and ID and is extended in (B). (B) K + transient has a two-stage rise, corresponding to PID and ID initiation, respectively. The ID initiation phase is extended in (C). (C) Red dashed lines indicate the start of the PID and the emergence of glutamate-mediated components at the beginning of the ID. The late stage of the ID consisted of short discharges, which are extended in (D). (D) Each of these late-stage discharges produces a transient increase of [ PLOS COMPUTATIONAL BIOLOGY point out that our experiments and similar previous studies mentioned in the Introduction reveal a challenging task to explain the mechanism of slow propagation of ID wavefront with the speed of the order of magnitude of tenths of a millimeter per second.

Simulation of the ictal front wave
Aiming to reproduce and explain the experimentally observed epileptiform activity, we have constructed a mathematical model (see the Methods section) that describes the main physiological mechanisms underlying interactions between neuronal populations in the cortical tissue. Below, we begin with the description of the simulated phenomenon, then pass to the mechanisms of ID generation and end up with the analysis of ID propagation.

Phenomenon
We simulated neuronal activity in a spatially extended cortical domain (Fig 4A, yellow  domain), thus obtaining the spatial-temporal distributions of the main variable, the membrane voltage, ionic concentrations, etc. (Fig 4B-4D) and the signals from one (Fig 5) or two spatially remote sites S1 and S2 (Figs 4E and 5C). The model reproduces the spontaneous generation of

PLOS COMPUTATIONAL BIOLOGY
repeating IDs and their spread. Similar to the experimental data, the model generates IDs every few minutes (mean ID frequency was 0.59 s -1 ; Fig 4F) that last for tens of seconds (Figs 4B-4E and 5). At the beginning of each ID, the outward GABAergic current predominates, and the simulated current at the holding voltage of -27 mV exposes short positive, i.e., GABAergic, bursts (Fig 4E, next to bottom; blue arrows in Fig 5A; compare to the experiment in Fig 1A). Then, the inward glutamatergic current prevails, exposing long-lasting negative, i.e., glutamatergic, components (Fig 4E, next to bottom; brown arrows in Fig 5A), similar to what is observed in the experiment (Figs 1 and 2A). These GABAergic and glutamatergic components constitute an ID. Each ID is characterized by the rapid increase in the concentration of the extracellular potassium ions, which peaked in simulations up to 14 mM (Fig 4F). The traces of [K + ] o at the sites S1 and S2, shown in Fig 4E, have similar shapes. The rapid growth of [K + ] o is delayed at one of the electrodes. It is in contrast to such traces obtained in the model  (E) signals from two sites, S2 (black) and S1 (blue): the extracellular potassium ion concentration, the membrane potential of two remote representative neurons, the current calculated at the "holding" voltage -27mV, and the reversal potential of that takes into account long-range (all-to-all) connections, which reveal a rather simultaneous initiation with weak components (Fig 5C), similar to that in experiments (Fig 3). As discussed below in the paragraph of the Results section "ID propagation in the networks with additional all-to-all connections," the long-range connections do not affect the bulk of IDs, and thus they are omitted in most of the simulations.

Mechanism of ID generation
We observe that in our conditions, interneurons do not inhibit but instead stimulate the pyramidal cells, whose activity then constitutes each new ID. Between IDs, the interictal and preictal GABAergic components seen as positive current events in  Fig 7J). The interneurons start to fire spontaneously before the ID (see small initial events in the firing rate trace for the I-population in Figs 6E and 7J). The I-neurons' initial firing is determined by the noise that mimics spontaneous synaptic activity, observed in the experiment and caused by increased synaptic transmitter release after 4-aminopyridine application. The enhanced I-neurons' firing increases the frequency of IPSCs (positive current events between IDs in Fig 5A and orange curve for the GABAergic current in Fig 6D).
The chloride movement into and out of cells is described as a sum of the passive conductive flux [31] and fluxes mediated by three proteins: the two cation-chloride-transporters, KCC2 and NKCC1, and the GABA A receptors (see I I GABA in Fig 6D). Their effect leads to a gradual accumulation of chloride ions inside neurons. The GABA reversal potential (Fig 4E, bottom, and Fig 6F) and the leak reversal potential (Fig 6F) increase due to chloride ion accumulation inside the cells. At every GABA-mediated event, the membrane potential tends to approach the GABA reversal potential, which is always above the chloride ion reversal potential ( Fig 6F) because of the bicarbonate ion permeability of GABA A receptors. In the case of membrane potential being above the chloride ion reversal potential, the activation of GABA A receptors increases the chloride ion entry into the neurons, even in the absence of depolarization by glutamatergic currents.
The KCC2 transporters diminish the difference between the reversal potentials for chloride and potassium ions. Thus, the elevation of [Cl -] i results in the elevation of [K + ] o and depolarization. The gradual depolarization of neuronal membranes and the depolarizing effect of spontaneous GABA-and glutamate-mediated events lead to excitation and synchronization of some I-neurons. The frequency of short bursts of I-cell excitation increases (see more significant events in the firing rate trace for the I-population in Fig 6E after t = 180 s).
The GABAergic activity leads to chloride ion accumulation in E-neurons (see the intracellular chloride ion concentration in Fig 6G), raising the reversal potential of GABA A receptormediated current and membrane potential. Therefore, I-neurons excite E-neurons, triggering glutamatergic events (see the E-population firing rate in Fig 6E and Fig 7I).
Glutamatergic excitation and elevation of [K + ] o act as positive feedback, thus maintaining the neural network's overexcitation. In opposition to glutamatergic drive, the synaptic resource, i.e., the pool of synaptic vesicles ready to release, is depleted ( Fig 7C) and interrupts firing, thus splitting ID into a series of bursts, LSDs ( Fig 7E). Sodium and potassium ion currents through voltage-gated and glutamatergic channels elevate [K + ] o and [Na + ] E i (Fig 6A and 6B). High [Na + ] i activates the Na + /K + pump. The pump provides a hyperpolarizing current and recovers the potassium transmembrane gradient, thus at two sites S1 and S2, located at one-fourth of the length from both ends of the simulated "slice" (Fig 4A), in the model that takes into account longrange (all-to-all) connections. The blue arrow marks the initiation of potassium response simultaneously at the two sites. The red notch shows the delay (about 2ms) used to estimate the propagation speed (0.6mm/s for the distance between the electrodes 1.2mm).
https://doi.org/10.1371/journal.pcbi.1009782.g005 This mechanism of ID generation is consistent with previous studies [21,32,33]. In particular, the dynamics of the neuronal membrane potential and the ionic concentrations is similar to the "best approximations based on the available literature," presented in [32], as seen from a comparison of the panels of Fig 8. [K + ] o peaks at the initial phase of ID with an undershoot accompanying the interval between IDs, as in Fig 1A. [Na + ] E i is likely to peak at the end of the ID [18]. [Cl -] E i reflects the recordings from [33,34]. The aforementioned activation of the Na + / K + pump explains the earlier peak of [K + ] o , the later peak of [Na + ] E i and the undershoot of

ID propagation
In simulations, as in our experiments, different IDs emerge in distinct parts of the slice. For the simulation shown in Fig 4, the sites of ID origination are marked by arrows in bottom panel B. For instance, the last two IDs ( Fig 4B) appear on opposite sides of the "slice." Except for the duration, the signals recorded in spatially remote sites are similar (Figs 4E and 7). The ID at the leading site lasts longer, such that the ID ceases almost simultaneously at different places, as often observed in experiments [2,35]. At the onset of ID, [K + ] o and the membrane potential in all types of neurons rise more gradually in the site-follower (Fig 7, blue curves) than in the leading site (black curves). In the later phase of ID, the LSDs appeared almost simultaneously in both sites. Any PIDs that were sometimes observed in the experiments were not present in the simulations.
While some of the IDs appear almost simultaneously in different locations, i.e., their apparent speed is large, the others propagate with limited speed. Our study focuses on the propagation of IDs at a minimum velocity. A few-second delay is observed between some of the IDs recorded at a 2-mm distance ( Fig 4E). For instance, the delay is about 5 s for the third ID shown in Fig 7A. Respectively, the ID wavefront speed is estimated to be 0.4 mm/s, which is a typical minimal speed among different IDs in different realizations of simulations with the fixed control settings (mean value 0.26 mm/s; Figs 7F and 9H). This value within an order of magnitude is consistent with our experimental estimates and the experimental data from other laboratories [3,10,11].

ID propagation in conditions of no potassium ion diffusion
Our mathematical model includes two mechanisms that may contribute to ID propagation: (i) the diffusion of potassium ions through extracellular space, described by Eq 27 and (ii) the propagation of action potentials through axons, followed by the transmission through synapses and propagation of postsynaptic signals by dendrites, described by the connection profile, Eq 33, with the synaptic conductance equations, Eqs 14-18. To separate the effects of the two mechanisms, we blocked potassium ion diffusion (last term in Eq 27). Still, we observed repeating ictal discharges ( Fig 9B) with a pattern that differed slightly but was qualitatively similar to the control case (compare Fig 9B to 9A). The time difference in the ID onsets was similar to that in the control case, suggesting that the ID propagation speed was about 1 mm/s or less. This evidence indicates that the diffusion of potassium ions is insignificant. It indicates potential and the reversal potentials of the total leak, chloride ion and GABAergic currents, and the reversal potential of the potassium ion current, for the E-and I-cells, respectively; (H) the ECS volume. https://doi.org/10.1371/journal.pcbi.1009782.g006

PLOS COMPUTATIONAL BIOLOGY
that the primary mechanism of ID propagation is based mainly on spatially distributed synaptic connections.

ID propagation in the networks with non-spreading connections
Next, we tested whether artificial elimination of synaptic connections' spatial spread has a decisive effect on discharge propagation. For this purpose, we changed the characteristic length of the Gaussian profile of the strengths of the connections λ to 0 in Eq 33, thus obtaining φ i,j (t,x) = ν i (t,x), i.e., the presynaptic firing rate φ i,j (t,x) at any position x for the population j is now determined by the firing activity ν i (t,x) of the corresponding neuronal population i in the same position x, no longer integrating the firing rates from nearby or remote neurons. Simulation has shown a dramatic change in the activity pattern ( Fig 9C). Still, at every position x at the same time interval 400 s, we observe four or five IDs; however, the discharges are not so synchronized as in the previous simulations. However, there are aligned patterns with a slope corresponding to the extremely slow propagation (between 0.1-0.01 mm/s). These effects seem similar to the spreading depression, which may spread through K diffusion, with almost no synaptic effects.
The dramatic difference between activity patterns in Figs 8A and 9 confirms the primary role of synaptic connections and their spread in wavefront propagation.
ID propagation in conditions of excessive potassium ion diffusion. Since the potassium ion diffusion may hypothetically provide the propagation of discharges, in the conditions of eliminated spatial spread of synaptic connections, we increased the diffusion coefficient 1000 times (390 instead of 0.39 μm 2 /ms) but maintaining λ!0. The simulation gave the spatially structured pattern shown in Fig 9E. It is consistent with the modeling results by Martinet et al. (2017) [22], where the diffusion coefficient was set to be 10 5 μm 2 /ms. Nevertheless, such a significant discrepancy in the numbers argues against the potassium ion hypothesis.
ID propagation in the networks with long connections. If the potassium ion diffusion's role is negligible, then the speed of the ID wavefront must depend on the spatial parameter λ in the governing equations. We increased λ fivefold and obtained the pattern with fast-propagating   (Fig 9G). The speed was 1.1 mm/s instead of 0.28 mm/s in the control case (Fig 9A) or 0.35 mm/s in a lack of potassium ion diffusion and short connections (Fig 9B). Hence, the connection length is the major factor that determines the speed of IDs. Neglecting potassium ion diffusion, the length of connections is the only spatial scale in the model, which enters through Eq 33. Hence the ID speed is just proportional to λ. Estimating the speed of the slowest waves in simulations with different values of λ and different realizations of noise (Fig 9H), we obtained the coefficient of the proportionality to be 4.2 ± 0.7 s -1 .

ID propagation in the networks with additional all-to-all connections
In the previous simulations, we considered Gaussian-like profiles of spatially-extended connections without long-range connections, which may affect ID propagation. We also suggest that the effect of long-range connections explains the appearance of PIDs observed in the experiments but missed in the simulations mentioned above. To test these hypotheses, we added homogeneous global connections to the Gaussian profile. To this end, a fraction of the presynaptic rate (0.2) was set to be proportional to the firing rate averaged across the entire domain. As a result, the activity pattern remained approximately the same as in the control case ( Fig 9D), but weak pre-ictal bursts became evident (see the arrow in Fig 9D). The PIDs are accompanied by the weak raise of [K] o , initiated almost immediately in the whole "slice," as seen from the traces at two sites ( Fig 5C). The weak [K] o growth is determined by the potassium outflux from the synapses activated by the long-range connections. Since the fraction of these connections is small, the weak raise of [K] o does not affect the bulk of ID. Thus, the longrange connections do not crucially affect the ID propagation but determine PIDs.

The effect of volume dynamics
The volume of extracellular space (ECS) changes significantly during IDs [36], which affects the discharges [19] and hypothetically might affect the ID propagation. Volume dynamics provide feedback for the modulation of ionic concentrations. A change in ionic balance results in a change in osmolarity and, consequently, volume change, which in turn primarily affects the extracellular concentrations. We describe the volume dynamics with Eq 32, and its effect is considered for the extracellular ion concentrations in Eqs 27-29. The ECS volume is much smaller than the intracellular volume, which is characterized by the ratio β. This ratio is close to the reverse ECS volume fraction estimated from the experiments. The ECS volume fraction ranges from 5% to 36% [37,38]. An effective ECS volume fraction is even less if we exclude the contribution of large reservoirs of ECS into the estimated value, so we set β = 10 in contrast to 7 from the previous modeling study [39]. Since the ECS volume is much smaller than the intracellular volume, any water flux out or inside cells leads to a more significant change of the ECS volume than of the intracellular volume. Therefore, it has a stronger effect on extracellular than intracellular ionic concentrations. The ECS volume reduces after a single ID at almost 50% and gradually recovers between IDs (Fig 6H), which is consistent with the experimental observations [36]. Such a volume change affects the extracellular ionic concentrations even if the intracellular concentrations remain constant. To check the ECS volume change effect on ID generation, we performed a fixed volume simulation (Fig 9F). fixed ECS volume, except for the [Cl -] o changes, which are bigger in the case of fixed volume. In this case, ID generation is more irregular, observed as a larger difference in interdischarge intervals (compare Fig 9F to 9A).

The effect of glia
To consider the effect of glia, we took into account a glial buffer in the form from [40]. Simulations revealed a quite expected increase of the intervals between IDs but no evident effects on the ID propagation (compare Fig 9I to 9A).

The phenomenon of epileptiform discharge propagation
We studied seizure propagation properties using an in vitro model of epileptiform activity and mathematical modeling. With electrode recordings made at two distant sites of the cortex, we found an ID propagation speed of less than 0.5 mm/s. In simulations, we reproduced the regime of repeating ID generation and propagation with speed close to the experimental estimates within an order of magnitude. Our experimental observations and simulations are consistent with data in patients and results obtained from in vivo and in vitro animal models. For instance, Martinet et al. (2015) [41] used invasive subdural electrode arrays covering a broad area of the cortex (8 × 8 cm) of patients with pharmacoresistant epilepsy to visualize seizure spread and observe waves of cortical region recruitment. They measured a mean recruitment speed of about 4 mm/s. Our estimations (Figs 3B and 4F) are also consistent with the previous EEG and fMRI studies in human patients, where propagation speeds varied from 0.2 to 10 mm/s [42][43][44][45]. The multiunit ECoG recordings showed the propagation to be 0.83 mm/s [7]. In in vitro models of epileptic activity with the zero magnesium (0 Mg 2+ ) solutions, the speed varied from 0.1 to 1.34 mm/s [12,46,47].
The speed of the IDs was crucially slower than that of the SDs, including LSDs (referred to as the afterdischarges in other studies [35]) and interictal discharges, which were measured to be > 30 mm/s [2,3,12]. This difference is also evident from our experiments and simulations and corresponds to our previous study on SDs [6]. The ID propagation is much lower than the axonal conduction velocities found in the brain, e.g., for mossy fibers (300 mm/s) [48]. From the other extreme, the ID speed is much faster than that of cortical spreading depression observed in migraine (3 mm/min) [49]. Our simulation with non-spreading connections ( Fig 9C) seems to be close to the case of cortical spreading depression. It shows similar speeds of discharges propagating purely by means of potassium diffusion. However, this issue requires further investigation.
Typically, in the late stage of an ID, the LSDs recorded at spatially remote sites are highly correlated [35]. They look synchronized on the time-scale of an entire ID because of their high speed. So, naturally, the IDs recorded at different sites typically cease simultaneously. That is why the ID spread studies reveal the tendency of IDs to be shorter at the locations where the ID front appears later [41], as also seen in our examples.
In our experiments, IDs were preceded by PIDs, which were not seen in the control simulation. However, after consideration of the long-range connections, such PIDs have been revealed in simulations, which clarifies the origin of the discharges: the PIDs were reflected through the long connections from the leading zone of ID generation.

Mechanism of ID generation
The ID origination scenario in our model supports the GABAergic hypothesis of seizure initiation [50,51]. This hypothesis proposes that synchronous activation of inhibitory interneurons underlies the inciting events that result in a seizure. In our case, the involvement of pyramidal cells goes through chloride accumulation [33,52] rather than through post-inhibitory rebound excitation [51,53]. The synchronous activation of interneurons is also supported by the chloride ion accumulation in these cells. That is why the ionic dynamics play a crucial role in the initiation of ID generation and maintenance of this activity. Therefore, it was essential to describe it in detail in our model. Briefly, intracellular chloride and extracellular potassium ion concentrations elevation provides positive feedback during a single ID generation. In contrast, the intracellular sodium ion accumulation does negative feedback terminating each ID through the activation of the sodiumpotassium pump [54]. Considering the fast negative feedback through the short-term synaptic depression and the calcium-dependent potassium channels, we obtain IDs in the form of clustered bursts of spikes. At the peak of depolarization during ID, neurons may show tonic firing or fall into the depolarization block. In total, the scenario is consistent with classical observations [1].

Mechanisms of ictal wavefront propagation
In the literature, several mechanisms of ID and SD propagation were suggested, the first two of which we considered in the present paper (Fig 10): (i) synaptic transmission between neurons, including spike propagation through axons and passive conduction of postsynaptic signals through dendritic neuronal branches [23,24]; (ii) the diffusion of potassium ions, as a substantial transient increase in extracellular potassium ion concentration is observed during IDs [1,22,25,26]; (iii) the ephaptic interactions of neurons through an electric field [55]; and (iv) the electrodiffusion of potassium or glutamate ions in the extracellular space [56]. As shown with modeling [6], IID and LSD propagation is determined by neuronal branching and synaptic transmission (the first mechanism) only. However, regarding IDs, with only this fast propagation mechanism being considered, it was unclear how the ictal wavefront travels within cortical tissue at such a slow speed as less than 1 mm/s. Moreover, the wavefront of potassium ion elevation usually accompanies an ID [1], as observed in our experiments. This observation points to the critical role of potassium ion elevation and its depolarizing effect. However, whether the potassium ions spread by diffusion is the cause of the ictal front's propagation (the second mechanism) was an open question. The speed of the potassium ion waves due to ionic diffusion is a constraint made by the effective diffusion coefficient that is affected by the tortuosity and estimated to be less than 1 mm 2 /ms [57][58][59]. It is uncertain if the potassium ion diffusion-based hypothesis and simulations [22] are consistent with the experiments. An artificially increased diffusion coefficient may provide speeds comparable to those of IDs, as in [22], where the coefficient was five orders of magnitude higher than the actual values (1 cm 2 /s). Recently, direct evidence concerning the lack of potassium ion diffusion effect on ID propagation from the hippocampus to the neocortex has been provided [60]. The authors reasoned that if activity in the follower territory is triggered by a rise in the number of potassium ions diffused from the other site, then for those events, the rise would appear to occur significantly earlier relative to the local firing. In fact, they have found no significant difference in the latencies, suggesting that the entrainment of neocortical events by the hippocampal activity in their preparation does not happen via the diffusion of potassium ions.
The third mechanism based on ephaptic interactions is known to determine the propagation of slow periodic epileptiform activity in the conditions of absent synaptic interactions [55]. However, it is generally too weak: the extracellular electrical potential with a typical magnitude (about 1 mV) can hardly increase the transmembrane potential up to its threshold value, which is usually much more than a few millivolts above the resting level. Nevertheless, the ephaptic interactions are supposed to trigger excitation in the case of reduced extracellular space [61]. The fourth mechanism of electrodiffusion is also too weak. As estimated earlier [28], the electric field gradients at an order of 2 mV per 100 μm may provide the speed of potassium ions of a few μm/s, which is negligibly small. Therefore, the last two mechanisms can be excluded from consideration. Thus, the present study was aimed to distinguish between the first two mentioned mechanisms of ID propagation through electrophysiological recordings and modeling.
Our simulations have helped to distinguish between the first two mechanisms (Fig 10). We have found that the potassium diffusion-based mechanism (Fig 10A) leads to a much lower speed of propagation, which refutes this hypothesis. On the contrary, the neuronal branchingbased mechanism (Fig 10B) determines the ictal front rate; thus, it plays a more important role. In this mechanism, the elevation of [K + ] o is also essential for the excitation; however, it might be localized, i.e., it follows the propagation of synaptic activity. These findings are consistent with our previous simplified modeling study [28]. We suggest that this mechanism is also the major one that explains the above-mentioned observations made in the in vitro models and patients; though, the large-scale recordings with EEG or fMRI electrodes may also reflect the contribution of the cortico-thalamo-cortical loops [41,[62][63][64][65][66].

Mechanism of propagation of SDs, which are much faster than IDs
As described in our previous study of SDs with modeling and electrophysiological tools [6], the mechanism of SD propagation includes only synaptic connections. It does not necessitate the ionic dynamics, though the pathologically impaired, quasi-constant level of intracellular chloride ion concentration is required for the generation of SDs. In our present simulations, the ionic concentrations were dynamic, but the SD mechanism was the same. Thus, it validates the conclusion of our previous study, stating that the SD propagation is independent of the ionic dynamics. In contrast, the ID propagation necessarily involves at least the dynamic change in extracellular potassium ion concentration. Consequently, the ID front propagates much slower than SDs, at the speed of a few tens of mm/s, so on the time scale of an entire ID, SDs are observed as almost synchronous events [2,3].

Factors affecting ID propagation
The present detailed model has revealed the main factors affecting propagation and some of the inefficient ones. The most influential parameter was the length of the connections. In our simulations, it increased the ID speed almost proportionally. It may also explain the bulk difference between the studies in humans and rodents: the typical length of collaterals in human brains is much longer. On the contrary, modulation of the potassium ion diffusion coefficient within its physiological range was inefficient. The volume dynamics and the glial buffering also did not contribute significantly to the ID propagation.
In conclusion, our study distinguishes between two primary hypotheses on the mechanisms of ictal wavefront propagation. It highlights the role of conventional signal propagation through neuronal branches. This axo-dendritic propagation is accompanied by positive feedback excitation provided by the elevation in the level of extracellular potassium ions that are released from excited neurons through voltage-gated and active glutamatergic receptors. These results are of certain importance to novel treatments against epilepsy.

Ethics statement
Three-week-old male Wistar rats (n = 15) were used in this study. The animals were kept under standard conditions with free access to food and water. All animal procedures followed the guidelines of the European Community Council Directive 86/609/EEC and were approved by the Sechenov Institute of Evolutionary Physiology and Biochemistry Bioethics Committee.

Brain slice preparation
The rats were sacrificed via decapitation, and their brains were removed rapidly. The brain slice preparation method has been described previously (Amakhin et al., 2016;. A vibrating microtome (Microm HM 650 V; Microm, Germany) was used to cut horizontal 350-μm-thick slices that contained the hippocampus and the adjacent cortical regions (including the entorhinal cortex (ERC) and the perirhinal cortex (PRC)). Artificial cerebrospinal fluid with the following composition (in mM) was used: 126 NaCl, 24 NaHCO 3 , 2.5 KCl, 2 CaCl 2 , 1.25 NaH 2 PO 4 , 1 MgSO 4 , and 10 dextrose. The artificial cerebrospinal fluid was aerated with a gas mixture of 95% O 2 and 5% CO 2 . All chemicals used to prepare the solutions were purchased from Sigma-Aldrich (St. Louis, MO, USA) unless stated otherwise. 1-2 slices per rat were used for the experiments.

In vitro model of epileptiform activity
Epileptiform activity in rat brain slices was induced using a pro-epileptic solution (

Whole-cell patch-clamp recordings
The recordings were performed at 30˚C. Pyramidal neurons in the deep layers of the medial ERC and layer 3 of PRC were visualized using a Zeiss Axioscop 2 microscope (Zeiss, Oberkochen, Germany) equipped with differential interference contrast optics and a video camera (Grasshopper 3 GS3-U3-23S6M-C; FLIR Integrated Imaging Solutions Inc., Wilsonville, OR, USA). Patch electrodes (3)(4)(5) were pulled from borosilicate glass capillaries (Sutter Instrument, Novato, CA, USA) using a P-1000 pipette puller (Sutter Instrument, Novato, CA, USA). A cesium methanesulfonate-based pipette solution (composition in mM: 127 CsMeSO3, 10 NaCl, 5 EGTA, 10 HEPES, 6 QX314, 4 ATP-Mg, and 0.3 GTP; pH adjusted to 7.25 with CsOH) was used for voltage-clamp recordings. For current-clamp recordings, a potassium gluconate-based pipette solution was used, in mM: 136 K-Gluconate, 10 NaCl, 5 EGTA, 10 HEPES, 4 ATP-Mg, and 0.3 GTP; pH adjusted to 7.25 with KOH. Whole-cell recordings were performed using a Multiclamp 700B (Molecular Devices, Sunnyvale, CA, USA) patch-clamp amplifier and an NI USB-6343 A/D converter (National Instruments, Austin, TX, USA) using WinWCP 5 software (University of Strathclyde, Glasgow, U.K.). The data were filtered at 10 kHz and sampled at 20 kHz. In all cells included in the sample, access resistance was less than 15 MO and remained stable (�20% increase) across the experiment. The liquid junction potential was compensated offline for the voltage-clamp recordings by subtracting 7 mV.

Extracellular potassium ion concentration recordings
Recordings of the extracellular potassium ion concentration were performed using the monopolar K + -selective microelectrodes [60]. The pipettes were pulled from borosilicate glass (Sutter Instrument). The pipettes' interior surface was exposed to hexamethyldisilazane vapor (Sigma-Aldrich) at 220˚C for 90 min. The pipettes were then backfilled with 100 mM KCl solution. A small volume of the K + sensor (Potassium Ionophore I, cocktail A; Sigma-Aldrich, cat. no. 99311) was taken into the salinized pipette's tip using slight suction. The recording of electrode voltage was performed using the Multiclamp 700B patch-clamp amplifier in currentclamp mode. We checked the stability of the electrodes at the start and end of each recording. Data from unstable electrode recordings were discarded. The extracellular K + concentration ([K + ] o ) at a given moment (t) was calculated from the electrode voltage, (V(t)), as follows: where S is the scaling factor, which was estimated by applying solutions with different [K + ] o at the tips of ionophore-filled electrodes using a fast application system (HSSE-2/3, ALA Scientific Instruments Inc., USA).
In all electrodes tested, the scaling factor was within a small range (0.043-0.045), so for all obtained recordings, S was set equal to the average value of 0.044mV -1 .

Experimental data analysis and statistics
The data analysis was performed using custom software written in Wolfram Mathematica 12 (Wolfram Research, Champaign, IL, USA). To estimate SLE propagation speed, a time lag between the corresponding K + transients in ERC and PRC was utilized. The time lag between the [K + ] o rises was measured at half the peak amplitudes. The characteristic distance was about 2 mm. The upper limit of ID speed propagation was calculated as a ratio of the distance between recorded neurons and the time lag. 2-4 IDs from each slice were recorded for speed estimation.
Sigmaplot 14 (Systat Software Inc., San Jose, CA, USA) was used for the statistical analysis of the results. Dixon's Q-test (at the 95% confidence level) was used to reject outliners. The Kolmogorov-Smirnov test was employed for the evaluation of the normality of sample data. For data that passed the normality test, the results were expressed as mean ± standard error of the mean. Otherwise, the results are expressed as the median and 25%-75% interquartile range (IQR).

A CBRD approach for populations of pyramidal neurons and interneurons
The model for connected populations uses a one-population model as a building block. Such a population receives input signals as synaptic conductances, which are determined by the presynaptic firing rates. The population's output signal is the firing rate, which affects other populations' presynaptic firing rates, as described by the equation of neuronal activity propagation. The population activity is also affected by the gradients of ionic concentrations, determining the ionic channels' driving forces.
The synaptically interacting neuronal populations model is based on our previous study [21]. We consider excitatory and inhibitory neuronal populations, denoted by indices E and I, respectfully, which are distributed in the 1-D space of the x-coordinate and connected by the synaptic AMPA, GABA, and NMDA receptors (Fig 11).

CBRD approach for a single population of neurons
The mathematical description of every single population is based on the probability density approach [67], namely, the CBRD approach [68]. This approach considers a population of an infinite number of Hodgkin-Huxley-like neurons receiving both a common input and an individual noise input for each neuron. In any arbitrary case of transient or steady-state stimulation, such a population's firing rate is approximated with a system of equations in partial derivatives, 1-D transport equations. The equations govern an evolution of neuronal states distributed in the phase space of the time elapsed since the last spike, t � . The equations include the Hodgkin-Huxley equations for the membrane voltage and gating variables, parameterized by t � , and the equation for the neuronal density in t � -space, ρ p (t,t � ), where the index p substitutes for E or I. The output characteristic of the population's activity is the firing rate ν p (t), which is equal to ρ p in the state of a spike, t � = 0. According to our previous work, the particular form of the equations written below is for regular-spiking pyramidal cells [69,70]. Basic neurons are single-compartment ones with the membrane voltage U p (t,t � ). Approximations of voltage-gated ionic currents are based on the neocortical pyramidal cell model from [71], including for E-neurons the calcium dynamics and calcium-dependent potassium current that provides an effect of slow spike timing adaptation. Parameterized by t � , the governing equations are as follows: @r p @t þ @r p @t � ¼ À r p HðU p ;g p tot Þ; ð1Þ where g p tot ðt;t � Þ is the total conductance, including the leak, voltage-gated and synaptic conductance g p syn ðt; U p Þ ¼ g AMPA;p ðtÞ þ g NMDA;p ðt; U E Þ þ g GABA;p ðtÞ.

Hazard function
The source term in the Eq 1 is the hazard function H. This function is defined as the probability to generate a spike for a single neuron that is actually in the state characterized by known neuronal state variables such as the mean membrane potential and gating variables and noise. The approximation of the hazard function H has been obtained for the case of white noise [69] and color noise [70] as a function of U(t) and g E tot ðt;t � Þ, and parameterized by the noise amplitude in the resting state s 0 V , the spike threshold voltage V th , and the ratio of the membrane time constant τ m = C/g tot to the noise time constant τ Noise , i.e. k = τ m /τ Noise : A ¼ 1 t m e 0:0061À 1:12TÀ 0:257T 2 À 0:072T 3 À 0:0117T 4 ð1 À ð1 þ kÞ À 0:71þ0:0825ðTþ3Þ Þ; where T is the membrane potential relative to the threshold, scaled by the noise amplitude σ V , which increases with the synaptic conductance: s V ¼ s 0 V ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 þ g syn =g L q . The term A is the hazard for a neuron to cross the threshold because of noise, derived analytically [69] and approximated by exponential and polynomial for convenience; B is the hazard for a neuron to fire because of depolarization due to deterministic drive, i.e., the hazard due to drift in the voltage phase space. Note that the H-function is independent of the basic neuron model and does not contain any free parameters or functions for fitting to any particular case. Thus, H-function is the same for excitatory and inhibitory populations.

Voltage-dependent channels
The set of ionic currents includes the voltage-dependent potassium currents I DR and I A responsible for spike repolarization, the slow potassium current I M that contributes to spike frequency adaptation and the potassium current I K−Ca , dependent on calcium dynamics and contributing to slow spike frequency adaptation. Approximating formulas for the currents I DR and I K−Ca are taken from [71].

Boundary conditions
According to the conservation of the number of neurons in a population, the firing rate is calculated as a sink of neurons from their state t � due to spiking, ρ p (t,t � )H(U p (t,t � ), integrated over the whole phase space, i.e.
n p ðtÞ � r p ðt; 0Þ ¼ It is the boundary condition for Eq 1.
The spike duration is taken into account by introducing the time interval 0<t � <Δt AP during which the voltage and the gating variables are fixed to their reset values. It defines the boundary conditions for Eqs 2-7 at t � = Δt AP which are as follows: The reset values for the fast gating variables in Eqs 10-12 were obtained with the basic single neuron model. With a rather arbitrary input providing a spike, these values were measured at the moment of a voltage maximum at the spike. The reset level for the slow conductance in the CBRD model was calculated as its value at the peak of spike-release distribution in the t � -space: where t �p is such that rðt; t �p ÞHðt; t �p Þ ¼ max Here g 0 tot is the total conductance at rest, and g syn is the total synaptic conductance; S is the membrane area; σ V is the noise amplitude meaning the dispersion of individual neuron's voltage fluctuations in a stationary state. Its scaling with g syn approximately reflects the synaptic noise increase with the increase of mean synaptic drive [72].
The noise term I noise in Eq 2 is nonzero only for the interneurons (Fig 11). It is calculated as an Ornstein-Uhlenbeck process with the amplitude 25pA and the correlation time 4ms.
The depolarisation block is modeled through the dynamic threshold as V E; 100ms . When calculating the dynamics of a neural population, the integration of Eqs 2-9 determines the evolution of the distribution of voltage U E across t � . Then, the effect of crossing the threshold and the diffusion due to noise are taken into account by H-function, Eq 4, substituted into the equation for neuronal density, Eq 1. The integral Eq 10 results in the output firing rate ν E (t).

Short-term synaptic plasticity
The synaptic depression was modeled with the Tsodyks-Markram model [74]: with τ glu = τ GABA = 500 ms, u glu = 0.2, and u GABA = 0.1. A crucial role of short-term synaptic depression in discharge termination was found in experiments and modeling [75]. Simulations with the absence of synaptic depression show that the neuronal populations spontaneously switch to a high-activity state without termination.

Representative neurons
Representative neurons of each of the populations were modeled by the basic single neuron model with the same synaptic inputs as for the populations. This model is described by the equations for the membrane voltage, Eqs 2,3,5-13, where the sum of partial derivatives was substituted by the total derivative in time t, and the sodium current was explicitly present in the right-hand part of Eq 2. The sodium current dependent on voltage V was approximated as in [71]: a m ¼ À 0:32ðU þ 50Þ=ðexpðÀ ðU þ 50Þ=4Þ À 1Þ; b m ¼ 0:28ðU þ 23Þ=ðexpðÀ ðU þ 23Þ=5Þ À 1Þ; with � g Na ¼ 7 mS=cm 2 and V Na = 50mV.

Ionic dynamics
The ionic concentrations that strongly affect reversal potentials are the extracellular potassium concentration [K + ] o , the intracellular chloride ion concentrations ½Cl À � E i and ½Cl À � I i for E-and Ineurons, respectively, and the intracellular sodium concentrations ½Na þ � E i and ½Na þ � I i (Fig 12). The potassium balance equation was modified after Wei et al. [76,77]. The contribution of inhibitory neurons α I was equal to 1/4 in accordance with the fraction of inhibitory neurons in the cortical tissue. For each of the populations, E and I, the intracellular concentrations of the potassium, chloride, and sodium ions were calculated according to the balance equations: d½Na þ � i dt ¼ g À I Na;leak À I Na;glu þ q Na nðtÞ À 3I pump þ I NKCC1 þ I NCX where β is the ratio of intra/extracellular volume; I E=I K;leak and I E=I K;active are the potassium currents through the leak and active voltage-gated channels; I E=I K;glu is the potassium current through glutamatergic channels; I E=I pump is the current of potassium dublets, provided by the Na + /K + pump; I E=I KCC 2 is the flux of potassium due to the KCC2-cotransporter; I I NKCC1 is the potassium flux due to the NKCC1-cotransporter taken into account only for the I-population; γ is the surface-tovolume and charge-to-concentration translating parameter; I Cl,leak and I Na,leak are the chloride and sodium ion leak currents, respectively; I GABA is the chloride ion current through GABA-A-controlled receptors; and I Na,glu is the sodium current through the glutamatergic receptors. Reversal potentials were obtained from the Nernst equations: The leak currents are as follows: is the average potassium current through the voltage-gated channels. The M-and K-Ca-components are absent for interneurons for simplicity. The potassium and sodium currents through the glutamatergic channels [81] were approximated as linear dependent on voltage with the fractions of the total glutamatergic conductance. The fraction was estimated as 0.2 for the potassium and 0.4 for the sodium ions [82]: The glial buffer is [40]: The parameters are as follows: g ¼ S=ðF volÞ ¼ 10 À 3 ½ðmM=msÞ=ðmA=cm 2 Þ�; β = 10; α I = 0.25. The Na + /K + pump, KCC2, and NKCC1 transporter current amplitudes are as follows: I pump,max = 0.8 μA/cm 2 The leak conductances are the same for both populations: g KL = 20 μS/cm 2 , g CIL = 20 μS/cm 2 , and g NaL = 7 μS/ cm 2 . The membrane area is 3�10 −5 cm 2 . The sodium charge transferred by a single spike is q Na = 0.1 μC/cm 2 . The potassium concentration mixture with the bath solution was characterized with the parameter D bath = 0.25s −1 . The coefficient of the spatial diffusion of potassium ions is D 1d = 0.39 μm 2 /ms. The initial concentrations are [K + ] o = [K + ] bath = 3.5mM, ½Cl À � The calcium concentration inside neurons increases mainly due to the influx through the NMDA receptors and decreases due to the outflow through the sodium-calcium exchangers. The effect of the exchanger is approximated as a relaxation term I E NCX ¼ ½Ca 2þ � E i =ðgt Ca Þ in the equation for the calcium concentration inside E-neurons where the time constant of the relaxation due to Na + /Ca 2+ exchanger τ Ca = 200 ms, the contribution of the calcium ions into glutamatergic ionic transport χ = 0.001 and the reversal potential V Ca = 20mV. The calcium dynamics was neglected for interneurons.

Volume equation
Extracellular space (ECS) volume dynamics has been described phenomenologically, as in [19]: where v ¼ v o =v 0 o is the ratio of the ECS volume to its initial value; β 0 = 10 is the initial ratio of intracellular versus ECS volumes; τ v = 250 ms is the characteristic time constant of the volume change because of the change of osmolarity Δπ: Dp ¼ ½Na þ � i þ ½Cl À � i þ ½K þ � i À ½Na þ � o À ½Cl À � o À ½K þ � o À Dp 0 , where

Equation of neuronal activity propagation
The horizontal cortical connections are supposed to be local and isotropic. They are determined by the relationship between the somatic rates v i and presynaptic firing rates φ i,j , where i and j are the indexes of the pre-and postsynaptic populations, respectively. In the case of 1-d geometry, all variables depend on the spatial coordinate x, oriented along with the layers of the cortical tissue. A Gaussian profile of the strengths of the connections is assumed, i.e.
where λ is the characteristic length, assumed to be equal for all types of connections, λ = 50 μm, which roughly corresponds to electrophysiological estimations from paired recordings [83]. In 1-d representation, the cortical area was considered as a segment of the length 2.5mm. The spatial inhomogeneity of parameters was set with the linear distribution of � g AMPA;E , which was equal to zero at x = 0, and maximum at x = L.

Software and source code
The model was implemented in the software "Brain" that is available from the public repository with DOI: https://doi.org/10.6084/m9.figshare.15113370.v1.

Estimation of activity pattern characteristics
To obtain statistical data from simulations, we varied realizations of the noise term I noise in Eq 2. From 5 realizations of 400s simulations, we used the second, third and fouth IDs, thus obtaining n = 15 values for each characteristic. We estimated the following characteristics of the activity: the speed of the wavefronts, the interdischarge frequency, and the maximum value of the extracellular potassium ion concentration reached during IDs. The first IDs after the onset of the simulation were excluded. For each ID, the speed was estimated from the slope of the slowest fragment of the wavefront, as shown with the black lines in Fig 4B-4D. Similar values were obtained in estimations of a time lag between the K+ transients in remote sites at half the peak amplitudes as indicated in Fig 7A; in this case, the ID speed was calculated as a ratio of the distance between recorded neurons and the time lag. The ID frequency was measured as the reverse time between the fronts at one location S1, so as the maximum extracellular potassium ion concentration. The results are expressed as the median and 25%-75% interquartile range.