A simple model of epileptic seizure propagation: Potassium diffusion versus axo-dendritic spread

The mechanisms of epileptic discharge generation and spread are not yet fully known. A recently proposed simple biophysical model of interictal and ictal discharges, Epileptor-2, reproduces well the main features of neuronal excitation and ionic dynamics during discharge generation. In order to distinguish between two hypothesized mechanisms of discharge propagation, we extend the model to the case of two-dimensional propagation along the cortical neural tissue. The first mechanism is based on extracellular potassium diffusion, and the second is the propagation of spikes and postsynaptic signals along axons and dendrites. Our simulations show that potassium diffusion is too slow to reproduce an experimentally observed speed of ictal wavefront propagation (tenths of mm/s). By contrast, the synaptic mechanism predicts well the speed and synchronization of the pre-ictal bursts before the ictal front and the afterdischarges in the ictal core. Though this fact diminishes the role of diffusion and electrodiffusion, the model nevertheless highlights the role of potassium extrusion during neuronal excitation, which provides a positive feedback that changes at the ictal wavefront the balance of excitation versus inhibition in favor of excitation. This finding may help to find a target for a treatment to prevent seizure propagation.


Introduction
Epilepsy is characterized by repeated seizures associated with abnormal intense electrical neural discharges. Despite ongoing research, the mechanisms of the generation and propagation of these discharges are not yet fully understood. Understanding these mechanisms is important for medical treatment development and helpful for mathematical modeling as an explanatory example of neuronal synchronization, which is a simpler regime of activity than normal functioning. On the other hand, in contrast to normal brain simulations, epileptic discharges involve the dynamics of ionic concentrations, thus requiring a more complex mathematical description.

PLOS ONE
PLOS ONE | https://doi.org/10.1371/journal.pone.0230787 April 10, 2020 1 / 21 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 The spread of activity through cortical circuits has been studied in experiments by means of electrical registrations and optical imaging [1][2][3], and high-density microelectrode arrays [4]. Experiments show slow propagation of an ictal wavefront and fast spread of discharges behind the front [3] [5]. The ictal wavefront progresses through the cortical area at a pace of < 1 mm/ s, which is consistent with propagation speeds measured with electrodes and imaging in brain slice models [1,2,[6][7][8][9] and in vivo (0.6 mm/s in [10] with two-photon microscope and 0.5 mm/s in [11] with widefield imaging in mouse neocortex). In the wake of the ictal wavefront, firing bursts are observed that are highly coherent across microelectrode recording sites within the same region, and correlated with ictal discharges recorded from adjacent macroelectrodes. The area demonstrating this hypersynchronized bursting is termed the ictal core [5].
The mechanism of the propagation is still an open question [3]. The main hypotheses consider the diffusion of potassium ions, the synaptic interactions, the electric field interactions, or the potassium electrodiffusion. As shown, the major role in excitability belongs to the extracellular potassium concentration. Recently, the spatial patterns of the extracellular potassium distribution have been registered by means of a nanoparticle-based technique [12]. The wavefront of potassium elevation from a seizure in the 4-AP (4-aminopyridyne, which strengthens synaptic connections) based model of cortical epilepsy spreads with a speed of about tenths of millimeters per second, which is similar to the typical speed of the ictal front. However, whether the potassium ion spread only accompanies the ictal front or explains the spread of the ictal front by diffusion or electrodiffusion is still a contested question. On one hand, observations of interictal discharges in slices show that correlated discharges exist even in the case of a lesion across a slice and are blocked after prevention of diffusion [13], thus pointing to non-synaptic, presumably potassium diffusion-based origin of the activity propagation or synchronization. On the other hand, the speed of potassium waves due to diffusion is limited [14,15], and it is uncertain if the potassium diffusion-based hypothesis and simulations [16] are consistent with the experimental estimations. Moreover, the mentioned experiments [13] describe not the behaviour of the ictal front but only the bursts inside the domain of spontaneous activity generation. In the present paper, by means of simulations we show that the potassium diffusion may cause the propagation of the excitation that leads to the discharge generation but do not determine the speed of ictal wavefronts.
Synaptic communications seem to play a dominant role in ictal front propagation. That is, the seizure propagation respects the excitatory network of connectivity underlying normal processing [11], and brain slice data shows that ahead of the ictal wavefront there are very large feedforward excitatory and inhibitory synaptic conductances [1,2,4]. The ictal wavefront generates huge feedforward excitation, yet a rapid feedforward inhibition provides a powerful restraint. Seizures in vivo propagate as slowly as the ictal wavefronts in the in vitro zero-magnesium model. However, epileptiform events recorded in disinhibited slices propagate much faster, between 10 and 200 mm/s, presumably because of little or no effective feedforward inhibition to slow propagation [1]. The exact mechanism by which the restraint is overcome during recruitment of cortical territories to a seizure remains to be determined, but may involve activity-dependent mechanisms that either boost excitatory neurotransmission or compromise inhibitory neurotransmission, or both [4]. Here, we suppose that an activity-evoked elevation of the extracellular potassium level might be a key factor in the change of the balance of excitation versus inhibition at the ictal front.
Some mathematical models that consider spatial propagation suggest that the excitability of the tissue surrounding the seizure core may play a determining role in the seizure onset pattern [17]. Whereas the generation of interictal discharges is modeled in the conditions of impaired-but-fixed ionic concentrations [18], the dynamics of ictal discharges and the excitability of the cortical tissue is hypothesized to be governed by the ionic dynamics [19], [20], [21]. Generally, a computational approach to this issue requires a biophysical consideration of the neuronal population interactions in the conditions of changing ionic concentrations of sodium, potassium, chloride, and calcium ions inside and outside the neurons and glial cells. This problem is quite complex and computationally expensive. The most well-elaborated biophysical models consider either a single neuron [22], or a network without a spatial structure [23,24], or a spatially structured but synaptically uncoupled neuro-glial network [19]. Taking into account the spatial propagation adds extra independent coordinates and thus essentially increases the complexity of models. For this reason, a consideration of spatial propagation requires a reduced but biophysically plausible model able to reproduce ictal discharges. Recently, we have proposed a spatially concentrated biophysical model of ictal and interictal discharges [25], called Epileptor-2 after the known abstract model Epileptor [26] that was further extended to large-scale simulations [27], [28]. Our model might also be extended to the spatially distributed case. In the present work, the Epileptor-2 model is extended by introducing the diffusion equation for the potassium concentration and the equation of spatial communications of neuronal populations via firing and postsynaptic propagation along axonal and dendritic trees [29].

Materials and methods
The previous implementation of Epileptor-2 [25] was restricted to the consideration of a spatially homogeneous activity. A primary advantage of that model is a biophysical interpretation of its governing variables that describe the pathological states of brain activity. For that purpose, the ionic dynamics equations used in earlier studies [23,30,31] are incorporated into a rate-based model for recurrently connected excitatory and inhibitory neuronal populations. The ionic dynamics description comprises equations for the concentrations of extracellular potassium and intracellular sodium. The firing rate of the inhibitory population is assumed to be proportional to that of the excitatory population, and thus, the inhibitory population is taken into account implicitly. The firing rate is described as a rectified sigmoid function of a membrane potential. The membrane potential is described by Kirchoff's current conservation law, which is written for a one-compartment neuron. The expressions for the excitatory and inhibitory synaptic currents, the input-output function, the rate-based equations for the ionic dynamics, etc., are justified in [25]. The short-term synaptic depression is described according to the Tsodyks-Markram model. An adaptive quadratic integrate-and-fire model is used to reveal a spiking activity of a representative neuron. Behavior of such representative neuron is quite similar to that of a real neuron during generation of ictal discharges (Fig 1).
For the spatially inhomogeneous case, the model is supplied with the equations for the extracellular potassium ion diffusion and the spike and synaptic current propagation along neuronal axons and dendrites. As in the previous implementation, the full system of equations is split into three subsystems that describe (i) the ionic dynamics, (ii) the neuronal excitability, and (iii) a neuron-observer. To be explained below, the equations for the balance of the extracellular potassium and intracellular sodium concentrations, [K] o and [Na] i , are as follows: The equation for the membrane potential V is The equation for the synaptic resource x D is The equation of spatial communications of neuronal populations via firing and postsynaptic propagation along axonal and dendritic trees is written for the presynaptic rate φ(t) governed Two ictal discharges as bursts of clustered interictal-like short bursts are seen in the membrane voltage. In the experiment, the ictal discharges were recorded in a pyramidal neuron from a rat entorhinal cortex using an in vitro 4-AP model of epileptic activity. Modified from [21]. https://doi.org/10.1371/journal.pone.0230787.g001 by the somatic rate ν(t): In the above equations, θ(t) represents the firing rate affecting potassium and sodium concentrations. It is different for each model of the mechanism of propagation (Models 1, 2, and 3; see the Results section). It is equal to either ν(t) or φ(t), where ν(t) is the firing rate on somas, and φ(t) is the firing rate on presynapses: yðtÞ ¼ nðtÞ ðfor Model 1Þ or yðtÞ ¼ φðtÞ ðfor Models 2 and 3Þ ð6Þ The somatic firing rate ν is calculated with a sigmoidal input-output function, where [x] + is equal to x for the positive argument and 0 otherwise: The input current u(t) includes the potassium depolarizing current, the synaptic drive, and the noise ξ, respectively: The potassium reversal potential is obtained from the ion concentrations via the Nernst equation: The functional form of Na + /K + pump current is taken from [31] (see also [32]) with modifications from [22]: Eq (1) includes the following terms: (i) the lateral diffusion term, (ii) the relaxation of [K] o to the bath concentration due to diffusion and glial buffering, (iii) the potassium influx due to ATP-dependent Na-K pump and (iv) the potassium elevation proportional to the firing rate, which cumulatively approximates the effects of the potassium-chloride cotransporters, the voltage-gated channels, and glutamatergic synapses, whose activation is roughly proportional to the firing rate. Eq (2) includes (i) the [Na] i relaxation term describing leakage and intracellular diffusion, (ii) the outflow due to the ATP-dependent Na-K pump, and (iii) the [Na] i elevation during the firing activity by means of the voltage-dependent and glutamatergic channels. Eq (3) describes a mean membrane depolarization due to the input u(t) using a single-compartment leaky neuron model. Not taking into account a spike generation, the voltage V(t) reflects a nominal, extreme level of membrane polarization. Together with Eq (8), it is a stochastic ordinary differential equation for the Ornstein-Uhlenbeck process. Eq (8) takes into account (i) the depolarizing effect of the potassium concentration increase relative to its initial concentration, (ii) the total synaptic current, and (iii) the noise. The current reflects the indirect effects of the ionic concentrations on membrane polarization through the voltage-gated and leak channels. The first term is a linearization of the dependence of all potassium channels (mainly, the leak) on the reversal potential deflection V K ðtÞ À V 0 K . The synaptic term consists of two parts, positive G syn φ(t)x D (t) and negative −c IE G syn φ(t). The positive component is an excitatory, short-term depressive current proportional to the presynaptic firing, thus implying an instantaneous synaptic kinetics. The short-term plasticity in the negative, inhibitory term is neglected; the coefficient c IE sets the balance of inhibition versus excitation. Eq (4) describes the short-term synaptic depression according to the Tsodyks-Markram model written in the rate-dependent form [33].
This new spatially extended version of the model includes the diffusion term (1), and the equation for the presynaptic firing rate, Eq (5). Eq (5) implies that the spatial communications of neuronal populations is provided by two factors: (i) the spread of action potentials through axons and (ii) the spatial integration of postsynaptic currents at the dendritic branches. Taking into account only local isotropic connections and implying that the dendritic propagation affects the presynaptic firing rate attributed to somatic location, the presynaptic firing rate φ is obtained as a convolution of the somatic firing rate ν with the exponentially decaying kernel exp(−r/λ), where r is the distance between pre-and postsynaptic neurons and λ is the characteristic length of the connectivity profile. Alternatively, Eq (5) can be derived from the mean field equations proposed in [29,34], assuming that the axonal delay is negligible. In our consideration, the axonal delay is much smaller than the membrane time constant which is the smallest time scale in the Epileptor-2 model, and thus it is to be neglected, which leads to Eq (5). We used the Neumann boundary conditions for Eqs (1) and (5), i.e. the fluxes of [K] o and φ were zeroed at all the boundaries.
The representative neuron is modeled with an adaptive quadratic integrate-and-fire neuron [35]. The equations for its membrane potential U and adaptation current w are as follows: Most of the parameters were set as in the previous version of Epileptor-2 [25]. They are given in Table 1. A few values were modified (given in bold font in Table 1).

Justification of main assumptions
The model describes only intracellular concentration of sodium and extracellular of potassium. First, the concentrations of intracellular potassium and extracellular chloride and sodium ions are assumed to be constant, because their effect on the network is much smaller than that for their counterparts. As in [21], we estimate the sensitivity of the network to changes of the concentrations. We evaluate the sensitivity with the absolute values of the derivatives of the reversal potentials on the concentrations, i.e. |@V ion /@ behaviour of these variables [36]. This issue and other assumptions are discussed in more details in [25].
We do not introduce into the models an electrodiffusion of potassium ions. In fact, the electric field produced by membrane currents on excited neurons at the front of the ictal discharge may accelerate the potassium ions in the extracellular space. This effect might be taken into account by the extra term u K V 0 @[K] o /@x in the equation for the extracellular potassium concentration, Eq (1), where u K is the ion mobility and V 0 is the voltage gradient [37]. Not considering the term with the second order derivatives, the equation would be a transfer equation of the form where "r.h.s." is the right hand side of Eq (1); the factor u K V 0 determines the speed of propagation of [K] o perturbations. Taking the value for u K from [37] to be equal to 5 � 10 −4 cm 2 /(s � V) and estimating the voltage gradient to be 2mV per 100μm, based on local field recordings [4], we obtain the speed to be of an order of a few micrometers per second, which is negligibly small in comparison to the seizure propagation speed. This result means that the electrodiffusion is too weak to contribute significantly to the potassium spread and can be omitted.

Numerical scheme and computer model
The simulations were performed in the Python 3.6 environment. The spatial derivatives in Eqs (1) and (5) have been approximated with the central difference scheme on a regular mesh with Neumann boundary conditions. The Euler-Maruyama explicit numerical scheme was applied for the integration in time of the stochastic ordinary differential equations and Eq (1). The linear algebraic system corresponding to the discretized Eq (5) was solved by the Jacobi method. The typical value of a time step was 1 ms. The computational domain was represented as a square with a side 6 mm and the circular center of excitation with R = 0.3 mm (Fig 2). The domain was discretized with a grid of 80x80 cells. The results were not significantly dependent on the numerical parameters and realizations of noise. The numerical realization of the model is available from the website gin.g-node.org/asanin/epilepsy-potassium-calculation.

Results
We have hypothesized two different mechanisms of the spatial propagation of epileptic activity. The first mechanism is based on the diffusion of potassium in the extracellular medium. The second is determined by the spread of spikes and synaptic currents through axons and dendrites isotropically distributed within the cortex. In order to distinguish between the two mechanisms, we have considered two different models of spatial propagation. Both models generalize Epileptor-2, the recently proposed spatially homogeneous model of epileptic activity [25]. The first model supplies Epileptor-2 with the diffusion term in the equation for the extracellular potassium concentration, Eq (1). Below we refer to this model as "the diffusion model," or Model 1. The second model adds to the system the equation of spiking activity spread, Eq (5), which connects presynaptic to somatic firing rates by assuming an exponentially decaying profile of connectivity with some characteristic length of the connections.
Below we refer to this model as "the synaptic model," or Model 2. Our simulations are aimed to reproduce the spatial-temporal patterns of activity of the cortical neural tissue during the generation of epileptic ictal discharges after local application of the proepileptic agent 4-AP [12]. In our models, we consider a square-shaped domain of nervous tissue with a small circular central zone being a source of epileptic discharges. In the central zone (Fig 2), the excitability G syn is elevated by setting G syn /g L = 5mV�s in contrast to G syn / g L = 1mV�s at the periphery. Both models are stochastic due to the introduced spatially homogeneous noise in Eq (8).

Temporal aspects of activity in the center of epileptic discharge generation
The two spatially distributed models show patterns of activity in the center of epileptic discharge generation similar to those in the original spatially homogeneous model Epileptor-2 [25]. Ictal (ID) and interictal discharges (IID) are reproduced. In the case of IDs, these quasiperiodic spontaneous events occur with an interburst interval of an order of minutes (Fig 3A). Each ID is characterized by a high-rate spiking activity, lasts about a few tens of seconds and consists of short bursts that resemble IIDs (Figs 3Ba and 4Ba); i.e. ID is a cluster of IID-like bursts. The membrane potential of the representative neuron (Figs 3Aa, 3Ba, 4Aa and 4Ba) and the concentrations of potassium and sodium ions (Fig 3Ac and 3Bc) reflect the spontaneous occurrence of the discharges.
[K] o dynamics determines the onset and the duration of an ID. The ID begins as soon as the slowly increasing [K] o reaches a certain threshold level (Figs  3Bc and 4Bc). Then, [K] o increases rapidly, because of intensive potassium extrusion through potassium voltage-gated and glutamatergic channels that are active during the ID; the activity of these channels is reflected in the firing rate (Figs 3Bd and 4Bd).
[K] o grows until it is balanced by the Na-K pump. The peak of [K] o takes place at the middle of the ID. After that, [K] o begins to decrease, finally returning to its baseline and even lower. The phase of the ID in which [K] o concentration approaches the baseline determines the termination of the ID. The Na-K pump is activated by the elevated intracellular sodium concentration. The sodium concentration increases because of high spiking and glutamatergic synaptic activity during IDs [21]. When a certain high level of the intracellular sodium concentration is reached, the potassium-sodium pump activates (Figs 3Bc and 4Bc). The Na + /K + pump peaks at the end of the ID. Its activity remains high until the baseline potassium concentration is restored. Then, the burst terminates, and the sodium concentration slowly decays to the original concentration before the next ID (Figs 3Ac and 4Ac).
From a mathematical point of view, the dynamics of IDs is determined by the extracellular potassium and intracellular sodium ionic concentrations which dynamics is close to slow deterministic oscillations, as shown in [25]. It is governed by a slow subsystem of the full system of equations, based on Eqs (2) and (1) without the diffusion term, both dependent on a time averaged firing rate that is a function of [K] o . On the contrary, the IID-like events that constitute IDs are governed by a fast subsystem, based on Eqs (3) and (4). These events represent spontaneous large-amplitude oscillations (Figs 3Ba and 4Ba).

Spatial aspects in Model 1: Diffusion
In the consideration of the extracellular potassium diffusion-based mechanism of activity spread, simulations were conducted with the diffusion equation Eqs (1) and (6) set to be θ = ν, thus omitting Eq (5). The diffusion coefficient was taken from [14].
The model's behavior at the central point is qualitatively similar to the spatially homogeneous case [25] and differs only quantitatively. The quasi-periodic spontaneous IDs occur at an interval of about 110-130s (Fig 3A). Each ID is characterized by a high rate of activity lasting  (Fig 2). B: An ictal discharge recorded in the center of the cortical domain. a: The representative neuron membrane depolarization U. The inset is a magnification of the first burst. b: The total input current u. c: The ionic concentrations [K] o and [Na] i , and the Na + /K + pump current I pump . d: The somatic firing rate ν and the synaptic resource x D . C: Potassium concentration spatialtemporal patterns during generation of two IDs. Two sites of "recordings" are marked as "S1" and "S2", remote at a distance of 2mm. D: A comparative plot of [K] o at two points S1 and S2. E: A comparative plot of U at the sites S1 and S2. F: The dependence of the front wave speed on the diffusion coefficient (obtained with a narrow computational domain).

PLOS ONE
about 17s (Fig 3B) and consisting of short bursts resembling IIDs (Fig 3Ba). The potassium threshold for ID initiation is about 4mM. Initiated at the center, the ID forms a radial wave, which spreads across the entire cortical domain, as seen from patterns of the extracellular potassium concentration (Fig 3C). [K] o rapidly increases at the front of the wave and more gradually decreases after the activation of the Na-K pump due to increasing sodium concentration, thus constituting the rear phase of the wave.
[K] o finally returns to its baseline and even lower. After reaching its minimum, the potassium concentration slowly increases toward the threshold of another ID initiation.
The wave profile remains approximately the same during its propagation, as seen from a comparison of [K] o evolution at two sites (Fig 3D), located in the center and the periphery and marked in Fig 3C. The moving pulse of [K] o corresponds to a single ID, formed as a cluster of IID-like discharges seen in the voltage plot in Fig 3E. The ID duration is approximately the same at the two sites. The amplitude of the [K] o pulses varies within 10% (Fig 3D). At the same time, the voltage patterns of IIDs are different, reflecting the spontaneous bursting character of the activity.
The velocity of the first potassium wave is about 0.035 mm/s. The second wave initiates at about 149s (Fig 3D). Its velocity is roughly the same, about 0.03 mm/s. The acceleration is due to the higher level of [K] o in front of the second wave, as seen from the comparison of the shapes of the first and second IDs detected at the same site S1 (Fig 3D). Still, the velocities obtained in physiological range of the diffusion coefficient values (Fig 3F) are much smaller than the experimental measures.

Spatial aspects in Model 2: Axo-dendritic spread
In the consideration of the synaptic mechanism of activity spread, simulations were conducted with the equation for the axo-dendritic propagation of spiking activity, Eqs (5) and (6) set to be θ = φ, thus without the diffusion term in Eq (1).
The model's behavior at the central point is also qualitatively similar to the spatially homogeneous case [25] and the previous scenario (Section 3.2), with only quantitative differences. The quasi-periodic spontaneous IDs occur at an interval of about 220s (Fig 4A). Each ID is characterized by a high rate of activity lasting about 40s ( Fig 4B) and consisting of short bursts resembling IIDs (Fig 4Ba). The potassium threshold for ID initiation is about 4mM. The spatial propagation is qualitatively similar to that in simulations with Model 1. Initiated at the center, the ID forms a radial wave, which spreads across the entire cortical domain (Fig 4C). [K] o rapidly increases at the front of the wave and more gradually decreases after the activation of the Na-K pump due to increasing sodium concentration (Fig 4Bc), thus constituting the rear phase of the wave.
[K] o finally returns to its baseline and even lower (shots at 170s and later in Fig 4C). After reaching its minimum, the potassium concentration slowly increases toward the threshold of another ID initiation after about 381s (Fig 4A and 4C).
The velocity of the first potassium wave is about 0.11 mm/s. The second moves with roughly the same velocity (Fig 4C). The speed is roughly proportional to the spatial scale of spatial connections λ (Fig 4F), being also dependent on the geometrical scale of the domain of excitation. The speed value is quite consistent with the range of experimental measures mentioned in the Introduction section. For instance, both the speed and the spatial scale are comparable with the speed and the size of neuronal clusters active during the discharges registered in [9].
The wave profile remains approximately the same during the propagation, as seen from a comparison of [K] o evolution at two sites (Fig 4D), located in the center and the periphery, and marked in Fig 4C. The moving pulse of [K] o , shown in Fig 4D, corresponds to a single ID.
At different points of time, the same ID forms different clusters of IID-like discharges seen in the voltage plot in Fig 4E. These clusters of spike bursts are different at the two sites, because of the spontaneous generation of IID-like events.
To clarify the spatial character of IID-like event generation and propagation, we have outlined a region of interest (ROI) in the form of a vertical strip aligned across the center of the simulation area (Fig 5A). The diagram of the activity in ROI versus the vertical coordinate y and time t (Fig 5B) demonstrates the effect of synchronization of the voltage bursts between the center of the ictal discharge generation zone and its periphery. During the time period of 0.4s, which is short in comparison with the time scales of the potassium wave propagation during an ID (Fig 4C), the voltage profile across y frequently fluctuates (Fig 5B). The fluctuations are seen as thin vertical stripes on the y−t diagram. Each such fluctuation is a burst that constitutes an ID. The outer envelope of the fluctuations, being averaged in time, corresponds to the front of the potassium wave. The fact that the bursts-stripes are vertical means that the activity within each stripe occurs simultaneously within some range of y. Alternatively, an event with a finite speed would have an inclined front. Thus the propagation of the voltage bursts is almost instantaneous in comparison with that of the potassium wave ( Fig 4C).
This effect of essentially different speeds of the waves of slow variables and fast (the potassium concentration and the membrane potential, correspondingly) is quite consistent with the experimental observations from [38] (see Fig 6 there). In their experiment, the wave of an intrinsic optical signal that reflected the cerebral blood volume was much slower than the wave of a voltage-sensitive dye signal. Suggesting that the cerebral blood volume follows the extracellular potassium concentration, our simulation explains that the slow wave originates as an envelope of fast voltage bursts and travels with the speed of IDs. The speed of the slow waves in our simulation (0.15 mm/s) was also comparable to that in their experiment (0.45 mm/s).
Because of slow propagation, the delays between the onsets of IDs at the points remote on scale of a millimeter are as large as a few seconds (Figs 4D and 6). By contrast, the afterdischarges and preictal discharges are much more precisely correlated. Small firing rate bursts precede an ID (Fig 6). They are either independent IIDs (for example, see the bursts at about t = 95s) or those that correspond to early bursts belonging to an ID that begins at some distance. In the latter case, the events simply reflect the excitation happening within the core of an ID. In both cases, the events do not recruit neurons into active spike generation that constitutes an ID. The real recruitment, characterized by an intense, sharp increase in the firing rate, occurs later.
Model 2 shows essentially larger waves than Model 1, as seen from a comparison of the domains with high [K] o in Fig 4C to those in Fig 3C. This difference is because of the higher speed and longer duration of IDs in Model 2. The large waves are more consistent with experiments. Also, the speed of the wave in Model 1 is too slow in comparison to measurements [1,39,40], which estimate a range of about 0.2-10 mm/s. These observations favor Model 2. The potassium diffusion mechanism is too slow to explain the ictal discharge propagation. In contrast, the propagation of spikes and synaptic currents by the axons and dendrites is a sufficient mechanism to explain the ictal discharge propagation.
As experiments have shown, the disinhibition dramatically increases the ictal wavefront speed up to 10 � 200 mm/s [4,41]. In our model, the contribution of the inhibition is determined by the factor c IE . In simulations, the decrease of this factor accelerates the propagation up to two orders of magnitude. (Variations of c IE in the range 0.25 � 0.8 preserve the wave- If the propagation depends on the synaptic connections, then damage to the connections is able to prevent the activity spread. We demonstrate this with a simulation of activity spread in a domain with a straight lesion set with the condition φ = 0 on a line segment (Fig 7). The activity cannot propagate across the lesion. Instead, it passes around the lesion and spreads through the entire area behind.

Model 3: Combination of the two mechanisms of activity spread
Our finding that the axo-dendritic connections provide faster propagation of epileptic discharges than potassium diffusion does prompt us to check whether the same fast propagation takes place in the situation where the both mechanisms are functioning together. To this end, we have considered a model that incorporates the entire system of Eqs (1)- (13). As suggested, simulations with the combined model have shown results quite similar to those from Model 2, which were described in the subsection 3.3. Therefore, taking the potassium diffusion into account in Model 2 does not change the solution, thus proving that potassium diffusion does not affect the spread of IDs. Axo-dendritic spread prevails over potassium diffusion.
In the case of the propagation in the domain with a lesion, as in the case shown in Fig 7 for the simulation with Model 2, the diffusion is potentially able to cross through the lesion, in contrast with the spiking activity which is stopped at the damage of the connections. However, within the time period of a single ID spread, the effects of the potassium diffusion are negligibly small, which is seen from a comparison of Figs 7 and 8 (compare the narrow regions of increased potassium near the lesion in two figures). Consequently, the potassium diffusion does not affect the spread of IDs.
The observed minor role of the potassium diffusion seems to be contradictory to the experiment performed by Lian et al. [13]. With a complete lesion of a slice by a blade they have detected epileptiform discharges in the two halves of the slice. The discharges were time-correlated, if the potassium diffusion was not impaired. The discharges were interictal-like ones, weaker and more frequent than ictal discharges. To reproduce the effect, we have set a complete lesion of synaptic connections on a line crossing the computational domain ( Fig 9A). With the basic set of parameters, the diffusion could not affect the discharge propagation across the lesion. We have modified some parameters of sensitivity to potassium, excluded any possible spatial correlations due to noise by supplying a spatially non-homogeneous noise, and increased the potassium diffusion coefficient to its characteristic value in pure electrolytes. As a result, we observed discharges in the entire domain (Fig 9C and 9E), which were weaker and more frequent than those in the previous simulations (compare Fig 9C with Figs 3D and 4D). They originate in the central, more excitable region and spread rapidly within the left subdomain. The elevated potassium diffuses through the lesion and evokes discharges in the subdomain behind the lesion. With the diffusion blocked, the right sub-domain behind the lesion keeps silent (Fig 9B and 9D). Therefore, the diffusion of the potassium ions is able to excite nearby regions, but with a significant delay, and moreover this effect is observed in a narrow range of simulation parameters. Summarizing, the potassium diffusion is able to transfer excitation but do not determine the speed of ictal wavefronts.

Discussion
With the help of a simple but biophysically meaningful model, we have considered two different hypotheses of ictal discharge propagation. The first one, Model 1, is based on the extracellular potassium diffusion. The second hypothesis, Model 2, is based on spatial propagation due to synaptic interactions. Both models reproduced ictal discharges as traveling waves that are similar to experimental ictal waves by many characteristics including the frequency, the duration, the composition of interictal-like bursts, the spontaneous character of bursts, the ionic dynamics and the dynamics of spiking. The main difference between the models is the speed of the waves. The diffusion model led to inconsistently thin and overly slow ictal waves, if the diffusion coefficient is constrained by its experimental estimations. Only an artificially increased diffusion coefficient may give comparable speeds of waves, as in [16], where the coefficient was 5 orders of magnitude higher than the real values (1 cm 2/ s). These observations are not in favor of the first hypothesis. On the contrary, the synaptic mechanism of propagation resulted in simulated ictal waves with a speed similar to experimental measurements in slices [1,2,6] as well as in vivo in animal models [10,11] and in human patients [4,5]. Further, we have shown that an addition of the diffusion mechanism to the synaptic one (Model 3) does not change the solution; i.e., the diffusion is too slow and cannot significantly contribute to the propagation. This comparison justifies the hypothesis that the synaptic transmission may fully determine the speed of the ictal wavefront.  The evolution of the extracellular potassium concentration at the three sites of the domain, S1, S2 and S3, in simulations with (C) and without (B) the diffusion. D and E. The spatial-temporal patterns of the extracellular potassium concentration during generation of two discharges in the left domain. The discharge in the right domain originates due to the potassium diffusion. Specific parameters were: τ K = 10s, g K,leak /g L = 3, K bath = 4mM and G syn /g L = 1mV�s in the periphery, K bath = 7mM and G syn /g L = 5mV�s in the central zone, a spatially nonhomogeneous noise, D K = 2 � 10 −5 cm 2/ s. https://doi.org/10.1371/journal.pone.0230787.g009

Role of feedforward inhibition in seizure propagation
Current experimental studies of seizure propagation reveal the role of feedforward inhibition [1,2,4,5,9,10]; however, a certain mechanism is still debatable [5]. During recruitment of new territories to an existing ictal event, the driving force is provided by the glutamatergic output of the seizing, "core" territories. Interneurons contribute to the feedforward inhibitory response. The ictal wavefront phenomenon coincides with a shift from predominately inhibitory to predominately excitatory currents. Ahead of this event, strong inhibitory currents effectively block pyramidal firing [1]. This "inhibitory veto" operates in the penumbra, limits the ictal wavefront's advance, and thus determines the speed of the wavefront. This mechanism is consistent with the behavior of our Models 2 and 3, but there is one detail that is different in our explanation. It was hypothesized that a collapse of inhibition takes place at the ictal wavefront, which is, however, difficult to reveal directly from experiments [1]. Our Model 2 reveals that it is not necessary to assume any additive factor such as the collapse of inhibition in order to explain the changes at the ictal wavefront. At the same time, the model supported high sensitivity of the speed of the wavefront to the intensity of inhibition, such that disinhibition is able to significantly speed up the front, as in [41].
Role of potassium-mediated positive feedback in seizure propagation Instead of feedforward inhibition, an increase of potassium concentration leads to the change in excitation-versus-inhibition balance. The potassium increases because of its outflux through voltage-gated and glutamatergic channels opened in the excited state at the ictal wavefront. This localized elevation of potassium concentration depolarizes the neuronal membranes, thus providing a positive feedback of excitation. This extra excitation changes the balance and forces the glutamatergic cells to fire and contribute to the common excitation underlying the ictal discharge. This scenario is consistent with the observations of highly correlated and fast onsets of electric activity and potassium increase [7].

Speed of preictal bursts and afterdischarges
The correlation of the onset moments of an ictal discharge at distal points is determined by the speed of the discharge. It is slow, and the delay between the signals at the points remote at a distance of about one millimeter is typically as large as a few seconds [4,5] (for instance, see Fig. 4a in [4]). By contrast, the afterdischarges and preictal discharges are much more precisely correlated. In electrophysiological recordings by multiple electrodes [4], there was a marked discrepancy between the apparent onset time of the ictal rhythm and the actual time when neurons were recruited. The earliest low-frequency ictal rhythms propagated rapidly across all electrodes, however, there was only an irregular and relatively low level of unit activity, and the real recruitment, characterized by an intense sharp increase in unit activity, occurred only later [4]. Similar weak synchronous bursts were observed before the onset of an ictal discharge in Model 2, which were simply reflections of the bursts that are generated inside the ictal core and compose the ictal discharge. Also, the afterdischarges are well correlated if observed as local field potentials [38] or intracellular signals [42]. They are also synchronized inside the ictal core in Model 2. Therefore, Model 2 is consistent with the description of the spreading seizure in terms of the ictal core and the penumbra [4,5].
A lesion prevents the ID propagation Our simulation of a spread in a domain with a lesion has shown that an ID cannot directly propagate across lesions, but instead, it passes around them. This finding seems to be contradictory to some experiments in slices [13], where epileptiform discharges were found to be synchronized on both sides of a complete cut. Noting that the observed discharges were not ictal discharges but spontaneously repeating weaker epileptiform events, we have suggested that their synchronization occurs because of the potassium diffusion but with a delay longer than that of natural ID propagation. The propagation of activity between two parts of the slice may be mediated by the potassium diffusion, which is otherwise too slow to affect an ID. Our simulation with a complete lesion has confirmed this suggestion.
Analogy to spreading depression waves Though we have found that potassium diffusion cannot be a major factor in the propagation of ictal wavefronts, the simulated slow waves in Model 1 are quite similar to concomitant waves of potassium and neuronal excitation observed in experimental models of spreading depression (SD) [43]. Obtained with a realistic potassium diffusion coefficient, the speed of waves is close to characteristic values of SD waves, of an order of a few mm/min [43,44]. As is known [44], the mechanism of SD is also based on potassium diffusion, which explains the observed similarity. However, the details of our model should not be directly translated to the case of SD.

Conclusion
The proposed model of ictal discharge generation and spread, found to be consistent with a majority of experimental observations, diminishes the role of potassium diffusion but supports the synaptic, potassium-mediated mechanism of ictal wavefront propagation. This prediction is important for the development of a treatment preventing seizure spread, and is to be further experimentally verified.