epilepsy-potassium-calculation/article/article.tex

\documentclass[10pt,letterpaper]{article}
\usepackage[top=0.85in,left=0.75in,footskip=0.75in]{geometry}
\usepackage[section]{placeins}
\usepackage{amsmath,amssymb}

% Use adjustwidth environment to exceed column width (see example table in text)
\usepackage{changepage}

% Use Unicode characters when possible
\usepackage[utf8x]{inputenc}

% textcomp package and marvosym package for additional characters
\usepackage{textcomp,marvosym}

% cite package, to clean up citations in the main text. Do not remove.
\usepackage{cite}

% Use nameref to cite supporting information files (see Supporting Information section for more info)
\usepackage{nameref,hyperref}

% line numbers
\usepackage[right]{lineno}

% ligatures disabled
\usepackage{microtype}
\DisableLigatures[f]{encoding = *, family = * }

% color can be used to apply background shading to table cells only
\usepackage[table]{xcolor}

% array package and thick rules for tables
\usepackage{array}

% package to indent first paragraphs
\usepackage{indentfirst}

% create \thickcline for thick horizontal lines of variable length
\newlength\savedwidth
\newcommand\thickcline[1]{%
\noalign{\global\savedwidth\arrayrulewidth\global\arrayrulewidth 2pt}%
\cline{#1}%
\noalign{\vskip\arrayrulewidth}%
\noalign{\global\arrayrulewidth\savedwidth}%
}

% \pdv command for partial differential notation
\newcommand*{\pdv}[3][]{\ensuremath{\frac{\partial^{#1} #2}{\partial {#3}^{#1}}}}

\def\del#1{{\textbf{}}}

% Text layout
%\raggedright
\setlength{\parindent}{0.5cm}
\textwidth 6.25in
\textheight 8.75in
% don't indent paragraphs and make a command for indentation when needed
%\newlength\tindent
%\setlength{\tindent}{\parindent}
%\setlength{\parindent}{0pt}
%\renewcommand{\indent}{\hspace*{\tindent}}


% Bold the 'Figure #' in the caption and separate it from the title/caption with a period
% Captions will be left justified
\usepackage[aboveskip=1pt,labelfont=bf,labelsep=period,justification=raggedright,singlelinecheck=off]{caption}
\renewcommand{\figurename}{Fig}

% Use the PLoS provided BiBTeX style
\bibliographystyle{plos2015}

% Remove brackets from numbering in List of References
\makeatletter
\renewcommand{\@biblabel}[1]{\quad#1.}
\makeatother



\usepackage{graphicx,caption}
\graphicspath{ {./images/} }
% Header and Footer with logo
\usepackage{lastpage,fancyhdr,graphicx}
\usepackage{epstopdf}
%\pagestyle{myheadings}
\pagestyle{fancy}
\fancyhf{}
%\setlength{\headheight}{27.023pt}
%\lhead{\includegraphics[width=2.0in]{PLOS-submission.eps}}
\rfoot{\thepage/\pageref{LastPage}}
\renewcommand{\headrulewidth}{0pt}
\renewcommand{\footrule}{\hrule height 2pt \vspace{2mm}}
\fancyheadoffset[L]{2.25in}
\fancyfootoffset[L]{2.25in}
\lfoot{\today}


\begin{document}
  \vspace*{0.2in}

  % Title must be 250 characters or less.
  \begin{flushleft}
  {\Large
  \textbf\newline{
  %Epileptic Seizure Propagation Across Cortical Tissue: Simple Model Based On Potassium Diffusion
  A Simple Model of Epileptic Seizure Propagation: Potassium Diffusion Versus Axo-Dendritic Spread}
  }
    \newline
    \\
    Anton V. Chizhov\textsuperscript{1,2*},
    Aleksei E. Sanin\textsuperscript{2}
    \\
    \bigskip
    \textbf{1} Computational Physics Laboratory, Ioffe Institute, Saint Petersburg, Russia
    \\
    \textbf{2} Laboratory of Molecular Mechanisms of Neural Interactions, Sechenov Institute of Evolutionary Physiology and Biochemistry of the Russian Academy of Sciences, Saint Petersburg, Russia
    %\\
    %\textbf{3} Affiliation Dept/Program/Center, Institution Name, City, State, Country

    \bigskip

    * anton.chizhov@mail.ioffe.ru

  \end{flushleft}

  \section*{Abstract}
  %Should be in Abstract: First the motivation. What task is actual and unsolved. Further methods, results and conclusion.


  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 theextracellular potassium diffusion, and the second is the propagation of spikes and postsynaptic signals along axons and dendrites.
  %Modeling show that both mechanisms produce similar in many aspects activity except of the spatial velocity and the length of epileptic discharge wave.
  As shown, thepotassium 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 it 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.


  %********************************************************************
  \section{Introduction}
  %********************************************************************

  % epilepsy for medicine and modeling
  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.

  % phenomenon of spread
  The spread of activity through cortical circuits has been studied in experiments by means of electrical registrations and optical imaging \cite{Trevelyan2006, Trevelyan2007, Smith2016}, and high-density microelectrode arrays \cite{Schevon2012}.
  Experiments show slow propagation of an ictal wavefront and fast spread of discharges behind the front \cite{Smith2016} \cite{Schevon2019}.
  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 \cite{Wong1990, Weissinger2000, Weissinger2005, Trevelyan2006, Trevelyan2007, Cammarota2013} and in vivo (0.6 mm/s in \cite{Wenzel2017} with two-photon microscope and 0.5 mm/s in \cite{Rossi2017} 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 \cite{Schevon2019}.

  % mechanisms
  The mechanism of the propagation is still an open question \cite{Smith2016}.
  The main hypotheses consider the diffusion of potassium ions, the synaptic interactions, the electric field interactions, or the potassium electro-diffusion.  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 \cite{Muller2018}. 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 electro-diffusion 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 \cite{Lian2001}, 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 \cite{Vern1977, Sykova1997}, and it is uncertain if the potassium diffusion-based hypothesis and simulations \cite{Martinet2017} are consistent with the experimental estimations. Moreover, the mentioned experiments \cite{Lian2001} describe not the behaviour of the ictal front but only the bursts inside the domain of spontaneous activity generation.


  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 \cite{Rossi2017}, and brain slice data shows that ahead of the ictal wavefront there are very large feedforward excitatory and inhibitory synaptic conductances \cite{Trevelyan2006, Trevelyan2007, Schevon2012}. 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 \cite{Trevelyan2006}. 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 \cite{Schevon2012}. 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.


  % models
  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 \cite{Wang2017}. Whereas the generation of interictal discharges is modeled in the conditions of impaired-but-fixed ionic concentrations \cite{Chizhov2017}, the dynamics of ictal discharges and the excitability of the cortical tissue is hypothesized to be governed by the ionic dynamics \cite{Gentiletti2017}, \cite{Chizhov2019}. 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 \cite{Wei2014a} or a network without a spatial structure  \cite{Bazhenov2004, Krishnan2011}. 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 \cite{Chizhov2018}, called Epileptor-2 after the known abstract model Epileptor \cite{Jirsa2014} that was further extended to large-scale simulations \cite{Proix2018}, \cite{Olmi2019}. Our model might also be extended to the spatially distributed case.
  % present model
  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 \cite{Jirsa1996}.


  %********************************************************************
  \section{Materials and methods}
  %********************************************************************

  The previous implementation of Epileptor-2 \cite{Chizhov2018} 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 \cite{Bazhenov2004, Kager2000, Cressman2009} 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 \cite{Chizhov2018}. 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.

  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. These equations will be discussed after the following introduction of the full model. The parameters are introduced in Table \ref{table1}.
  %The diffusion term was included into potassium ion dynamics with the diffusion coefficient of the same value as in \cite{Bazhenov2004}. And for the spike propagation a separation of the general population firing rate on soma's and presynaptic's firing rates was considered. In that concern presynaptic's firing rate depends on soma's one and represents its value after axial propagation.
  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.
  The equations for the balance of the extracellular potassium and intracellular sodium concentrations, $[K]_o$ and $[Na]_i$, are as follows:
  \begin{eqnarray}
    \label{eqn:K}
    \pdv{[K]_o}{t} &=& D_K\biggl(\pdv[2]{[K]_o}{x} + \pdv[2]{[K]_o}{y}\biggr) + \frac{[K]_{bath} - [K]_o}{\tau_K} - 2\gamma ~I_{pump}(t) + \delta_{K}\theta(t),\\
    \label{eqn:Na}
    \frac{d[Na]_i}{dt} &=& \frac{[Na]^0_i - [Na]_i}{\tau_{Na}} - 3I_{pump}(t) + \delta_{Na}\theta(t).
  \end{eqnarray}
  The equation for membrane potential $V$ is
  \begin{equation}
    \label{eqn:V}
    C\frac{dV}{dt} = -g_L V + u(t)
  \end{equation}
  The equation for the synaptic resource $x^D$ is
  \begin{equation}
    \label{eqn:xD}
    \frac{dx^D}{dt} = \frac{1-x^D}{\tau_D} - \delta_{x} ~x^D \varphi(t)
  \end{equation}
  The equation of spatial communications of neuronal populations via firing and postsynaptic propagation along axonal and dendritic trees is written for the presynaptic rate $\varphi(t)$ governed by the somatic rate $\nu(t)$:
  \begin{equation}
    \label{eqn:phi}
    \pdv[2]{\varphi}{x} + \pdv[2]{\varphi}{y} =  \frac{\varphi-\nu}{\lambda^2} %\tag{5.1}
  \end{equation}
  In the above equations,
  $\theta(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 $\nu(t)$ or $\varphi(t)$, where $\nu(t)$ is the firing rate on somas, and $\varphi(t)$ is the firing rate on  presynapses:
  \begin{equation}
    \label{eqn:theta}
    \theta(t) = \nu(t) \text{ (for Model 1)~~~ or ~~~} \theta(t) = \varphi(t) \text{ (for Models 2 and 3)} %\tag{5.2}
  \end{equation}
  %The choice of $\theta$ in Eq. (\ref{eqn:theta}) depends on a model of the mechanism of spatial propagation.
  %For the case of diffusion it is $\theta = \nu$. For the case of presynaptic firing rate it is $\theta = \varphi$.
  %Inhibitory population firing rate is assumed to be proportional to $\nu$.
  The somatic firing rate $\nu$ is calculated with a sigmoidal input-output function, where $[x]_+$ is equal to $x$ for the positive argument and $0$ otherwise:
  \begin{equation}
    \label{eqn:nu}
    \nu(t) = \nu_{max}\Bigg[\frac{2}{1+\exp\Big[-2(V(t)-V_{th})/k_v\Big]}-1\Bigg]_{+}
  \end{equation}
  %$[K]_o$ and $[Na]_i$ are extracellular potassium concentration and intraneuronal sodium concentration respectively. $V$ is the membrane depolarization. $x^D$ is the synaptic resource.
  The input current $u(t)$ includes the potassium depolarizing current, the synaptic drive, and the noise $\xi$, respectively:
  \begin{equation}
    \label{eqn:uu_diff}
    u(t) = g_{K,leak}(V_K(t) - V^0_K) + G_{syn}\varphi(x^D(t)-c_{IE}) + \sigma \xi(t)
  \end{equation}
  The potassium reversal potential is obtained from the ion concentrations via the Nernst equation:
  \begin{equation}
    \label{eqn:V_K}
    V_K(t) = 26.6mV\ln\biggl(\frac{[K]_o}{130mM}\biggr)
  \end{equation}
  The $Na^+/K^+$ pump current is taken from \cite{Cressman2009} in the form:
  \begin{equation}
    \label{eqn:I_pump}
    I_{pump}(t) = \frac{\rho}{\big(1+\exp(3.5-[K]_o)\big)\big(1+\exp\big(({25-[Na]_i})/{3}\big)\big)}
  \end{equation}
  The representative neuron is modeled with an adaptive quadratic integrate-and-fire neuron \cite{Izhikevich2003}. The equations for the membrane potential $U$ and the adaptation current $w$ are as follows:
  \begin{equation}
    \label{eqn:U}
    C_U\frac{dU}{dt} = g_U(U-U_1)(U-U_2) - w + u + I_{a}
  \end{equation}
  \begin{equation}
    \label{eqn:w}
    \tau_{w}\frac{dw}{dt} = - w
  \end{equation}
  \begin{equation}
    \label{eqn:U_reset}
    \text{if } U > V^T \text{ then } U = V_{reset}, w = w + \delta w
  \end{equation}

  %\subsubsection{The extensions to the original Epileptor-2 model}
  This new spatially extended version of the model includes the diffusion term $D_K\big(\partial^2{[K]_o}/\partial{x^2} + \partial^2{[K]_o}/\partial{y^2}\big)$ in the equation for $[K]_o$, eq.(\ref{eqn:K}), and the equation for the presynaptic firing rate, eq.(\ref{eqn:phi}).
  Eq.(\ref{eqn:phi}) 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 $\varphi$ is obtained as a convolution of the somatic firing rate $\nu$ with the exponentially decaying kernel $\exp(-r/\lambda)$, where $r$ is the distance between pre- and postsynaptic neurons and $\lambda$ is the characteristic length of the connectivity profile.
  Alternatively, eq.(\ref{eqn:phi}) can be derived from the mean field equations proposed in \cite{Jirsa1996}, \cite{Robinson1997}, assuming that the axonal delay is negligible.
  %The equation of spatial communications of neuronal populations via firing and postsynaptic propagation along axonal and dendritic trees has been proposed in  \cite{Jirsa1996}.
  %Spatial communications of neuronal populations is provided by two factors, the spread of action potentials through axons and 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, one obtains an equation that determines the presynaptic rate $\varphi$ as governed by the rate $\nu$, which is as follows:
  %According to \cite{Jirsa1996} presynaptic $\varphi$ and soma $\nu$ firing rates can be linked as in the equation below:
  \del{
  \begin{equation}
    \label{eqn:phi_wave}
    \pdv[2]{\varphi}{t} + 2\omega\pdv{\varphi}{t} + \omega^2\varphi - c^2(\pdv[2]{\varphi}{x} + \pdv[2]{\varphi}{y}) = (\omega^2 + \omega\pdv{}{t})\nu
  \end{equation}
  \begin{equation}
    \label{eqn:omega}
    \omega = \frac{c}{\lambda}
  \end{equation}
  Where $c$ is the mean velocity of cortical pulse propagation from a soma to a presynapse, $\lambda$ is the average spatial range of cortical fibers connectivity. Assuming that $\frac{\lambda}{c} << \tau_m$ where $\tau_m$ is synaptic constant we can propose that when $\omega \rightarrow \infty$ and $c \rightarrow \infty$. Thus we can keep only quadratic terms of $\omega$ in (\ref{eqn:phi_wave}) and simplify it to the below form:
  \begin{equation}
    \label{eqn:phi_puass}
    \omega^2\varphi - c^2(\pdv[2]{\varphi}{x} + \pdv[2]{\varphi}{y}) = \omega^2\nu  \tag{\ref{eqn:phi_wave}.1}
  \end{equation}
  By dividing (\ref{eqn:phi_puass}) on $\omega^2$ and applying (\ref{eqn:omega}) $\frac{c}{\omega} = \lambda$ we receive the following equation which is just another form of (\ref{eqn:phi}).
  \begin{equation}
    \varphi - \lambda^2(\pdv[2]{\varphi}{x} + \pdv[2]{\varphi}{y}) = \nu \tag{\ref{eqn:phi_wave}.2}
  \end{equation}
  }
  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.(\ref{eqn:phi}).

  We used the Neumann boundary conditions for eq.(\ref{eqn:K}) and eq.(\ref{eqn:phi}), i.e. the fluxes of $[K]_o$ and $\varphi$ were zeroed at all the boundaries. Most of the parameters were set as in the previous version of Epileptor-2 \cite{Chizhov2018}. They are given in Table \ref{table1}. A few values were modified (given in bold font in Table \ref{table1}).
  %to provide a in bold are those that differ from the original ones or of the new parameters.

  \begin{table}
    \newlength\q
    \setlength\q{\dimexpr .25\textwidth - 3\tabcolsep}
    \noindent\begin{tabular}{p{\q}p{\q}p{2\q}}
               \hline
               Parameter & Value with units & Description \\
               \hline
               $\tau_K$ & 100s & Potassium time constant \\
               $\tau_{Na}$ & 20s & Sodium time constant \\
               ${C}/{g_L}$ or $\tau_m$ & 10ms & Membrane time constant \\
               $\tau_D$ & 2s & Synaptic depression time constant \\
               $D_K$ & $\pmb{4\cdot10^{-6}}$cm\textsuperscript{2}/s & Potassium diffusion coefficient \\
               $\delta_K$ & \textbf{0.04}mM & Potassium concentration increment at spike \\
               $\delta_{Na}$ & 0.03mM & Sodium concentration increment at spike \\
               $\delta_x$ & 0.01 & Synaptic concentration increment at spike \\
               $\xi$ & $\mathcal{N}(0,1)$ & Gaussian white noise \\
               ${\sigma}/{g_L}$ & 25mV & Noise amplitude \\
               $\rho$ & 0.2mM/s & Maximum pump flux \\
               $\gamma$ & \textbf{20} & Volume ratio \\
               ${G_{syn}}/{g_L}$ & 5mV$\cdot$s & Postsynaptic charge \\
               $c_{IE}$ & 0.5 & Inhibitory vs. excitatory conductance ratio \\
               ${g_{K,leak}}/{g_L}$ & \textbf{1.0} & Potassium leak conductance \\
               $[K]_o^0$ & 3mM & Initial extracellular potassium concentration \\
               $[K]_{bath}$ & \textbf{7}mM & Extracellular potassium concentration \\
               %$V_K^0$ & -100mV & Initial potassium reversal potential \\
               $[Na]_i^0$ & 10mM & Initial intracellular sodium concentration \\
               $\nu_{max}$ & 100Hz & Soma maximal firing rate \\
               $\lambda$ & \textbf{0.385}mm & Spatial scale of cortical fiber connectivity \\
               $V_{th}$ & 25mV & Threshold potential \\
               $k_{\nu}$ & 20mV & Gain \\
               $g_U$ & 1.5nS/mV & Representative neuron conductivity \\
               $C_U$ & \textbf{1050}pF & Representative neuron capacity \\
               $V^T$ & 25mV & Representative neuron threshold potential \\
               $V_{reset}$ & -40mV & Representative neuron reset potential \\
               $U_1$ & -60mV & Representative neuron reversal potential 1 \\
               $U_2$ & -40mV & Representative neuron reversal potential 2 \\
               $U^0$ & -70mV & Representative neuron initial potential \\
               $I_{a}$ & 116pA & Representative neuron tonic current \\
               $\tau w$ & 200ms & Adaptation current time constant \\
               $\delta w$ & 100pA & Adaptation current increment ap spike \\

               \hline
    \end{tabular}
    \caption{Parameters of the models. Basic parameter values are from \cite{Chizhov2018}, except the modified values that are given in bold.}
    \label{table1}
  \end{table}

  \subsection*{Numerical scheme and computer model}
  The simulations were performed in the Python 3.6 environment. The spatial derivatives in eqs. (\ref{eqn:K}) and (\ref{eqn:phi}) 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. (\ref{eqn:K}). The linear algebraic system corresponding to the discretized eq. (\ref{eqn:phi}) 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.\ref{fig:scheme}). The domain was discretized with a grid of 80x80 cells. The results were not significantly dependent on the numerical parameters; they were quantitatively dependent  on realizations of noise.
  %To preserve the results reproduction the noise seed was fixed at the value 66.
  The numerical realization of the model is available from the website \href{https://bitbucket.org/vogdb/epilepsy-potassium-calculation/}{bitbucket.org/vogdb/epilepsy-potassium-calculation/}.

  \begin{figure}
    \centering
    \includegraphics[width=0.4\textwidth]{images/Fig0_scheme.png}
    \caption{Computational domain includes a 6mm by 6mm piece of the cortical tissue. The central circular domain with $R=0.3~$mm (red) is a center of excitation due to high conductance $G_{syn}$, in contrast to the low conductance on the periphery (blue).}
    %Fig. 0
    \label{fig:scheme}
  \end{figure}

  %********************************************************************
  \section{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 \cite{Chizhov2018}. The first model supplies Epileptor-2 with the diffusion term in the equation for the extracellular potassium concentration, eq.\ref{eqn:K}. 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.\ref{eqn:phi}, 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.

  % Электродиффузия
  We do not introduce into the models an electro-diffusion of potassium ions. 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~v'~\partial [K]_o/\partial x $
  in the equation for the extracellular potassium concentration, eq.\ref{eqn:K}, where $u$ is the ion mobility and $v'$ is the voltage gradient \cite{Hodgkin1953}. Taking the value for $u$ from \cite{Hodgkin1953} to be equal to $5\cdot 10^{-4}cm^2/(s\cdot V)$
  and estimating the voltage gradient to be $2mV$ per $100\mu m$, based on local field recordings \cite{Schevon2012}, 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 electro-diffusion is too weak to contribute significantly to the potassium spread and can be omitted.

  % Постановка задачи для расчётов
  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 \cite{Muller2018}. 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.\ref{fig:scheme}), the excitability $G_{syn}$ is elevated by setting $G_{syn}/g_L=5$mV$\cdot$s in contrast to $G_{syn}/g_L=1$mV$\cdot$s at the periphery.
  Both models are stochastic due to the introduced spatially homogeneous noise in Eq.\ref{eqn:uu_diff}. %The nature of this parameter can be homogeneous or heterogeneous across the spatial area. For the simulations we use the former homogeneous type of noise as the latter type does not affect resulting epileptic waves except the frequency of their appearance. With the homogeneous noise epileptic waves appear two times often than with the heterogeneous noise thus requiring less simulation time. All other wave characteristics stay the same.

  % 0-мерная модель
  \subsection{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 \cite{Chizhov2018}.
  Ictal (ID) and interictal discharges (IID) are reproduced.
  In the case of IDs, these quasi-periodic spontaneous events occur with an interburst interval of an order of minutes (Fig.\ref{fig:diff_K_point_center_period}). 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 (Fig \ref{fig:diff_K_point_center}A, \ref{fig:syn_K_point_center}A); i.e. ID is a cluster of IID-like bursts.
  The membrane potential of the representative neuron (Fig.\ref{fig:diff_K_point_center_period}A,\ref{fig:diff_K_point_center}A,\ref{fig:syn_K_point_center_period}A,\ref{fig:syn_K_point_center}A) and the concentrations of potassium and sodium ions (Fig. \ref{fig:diff_K_point_center_period}C,\ref{fig:syn_K_point_center_period}C) 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 (Fig.\ref{fig:diff_K_point_center}C,\ref{fig:syn_K_point_center}C). 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 (Fig.\ref{fig:diff_K_point_center}D,\ref{fig:syn_K_point_center}D). $[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 \cite{Chizhov2019}.
  When a certain high level of the intracellular sodium concentration is reached, the potassium-sodium pump activates (Fig.\ref{fig:diff_K_point_center}C,\ref{fig:syn_K_point_center}C).
  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 (Fig. \ref{fig:diff_K_point_center_period}C,\ref{fig:syn_K_point_center_period}C).


  \begin{figure}
    \centering
    \includegraphics[width=1.0\textwidth]{diff/x40y40_period.png}
    \caption{Diffusion mechanism (Model 1). Repeated IDs recorded in the center of the cortical domain (Fig.\ref{fig:diff_K_board}). \textbf{A}, The representative neuron membrane depolarization $U$. \textbf{B}, The total input current $u$. \textbf{C}, The ionic concentrations $[K]_o$ and $[Na]_i$, and the $Na^+/K^+$ pump current $I_{pump}$. \textbf{D}, The somatic firing rate $\nu$ and the synaptic resource $x^D$.}
    %Fig. 1
    \label{fig:diff_K_point_center_period}
  \end{figure}

  \begin{figure}
    \centering
    \includegraphics[width=1.0\textwidth]{diff/x40y40.png}
    \caption{Diffusion mechanism (Model 1). An ictal discharge recorded in the center of the cortical domain (Fig.\ref{fig:diff_K_board}). \textbf{A}, The representative neuron membrane depolarization $U$. \textbf{B}, The total input current $u$. \textbf{C}, The ionic concentrations $[K]_o$ and $[Na]_i$, and the $Na^+/K^+$ pump current $I_{pump}$. \textbf{D}, The somatic firing rate $\nu$ and the synaptic resource $x^D$.}
    %Fig. 2
    \label{fig:diff_K_point_center}
  \end{figure}

  \begin{figure}
    \centering
    \begin{minipage}{0.5\linewidth}
      \centering
      \captionsetup{width=0.8\linewidth}
      \includegraphics[width=1.0\linewidth]{diff/K_points.png}
      \caption{Diffusion mechanism. A comparative plot of $[K]_o$ at two points S1 and S2 shown in Fig. \ref{fig:diff_K_board}.}
      \label{fig:diff_K_points}
    \end{minipage}\hfill
    \begin{minipage}{0.5\linewidth}
      \centering
      \captionsetup{width=0.8\linewidth}
      \includegraphics[width=1.0\linewidth]{diff/U_points.png}
      \caption{Diffusion mechanism (Model 1). A comparative plot of $U$ at the sites S1 and S2.}
      %Fig. 3
      \label{fig:diff_V_points}
    \end{minipage}
  \end{figure}

  %The proposed here spatially extended models show the same characteristic features of activity in the center of epileptic discharge generation as the spatially homogeneous model.
  From a mathematical point of view,  the dynamics of IDs is determined by the quasi-periodic oscillations of the extracellular potassium and intracellular sodium ionic concentrations, as shown in \cite{Chizhov2018}. These concentrations are governed by  a slow subsystem of the full system of equations, based on Eq.(\ref{eqn:K}) without the diffusion term and Eq.(\ref{eqn:Na}). On the contrary, the IID-like events that constitute IDs are governed by a fast subsystem, based on Eqs.(\ref{eqn:V}, \ref{eqn:xD}). These events represent spontaneous large-amplitude oscillations (Fig. \ref{fig:diff_K_point_center}A,\ref{fig:syn_K_point_center}A).


  %*************************************************************************************************
  % Модель 1
  \subsection{Spatial aspects in Model 1: Diffusion}
  \label{Results_Model1}

  In the consideration of the extracellular potassium diffusion-based mechanism of activity spread, simulations were conducted  with the diffusion equation eq.(\ref{eqn:K}) and eq.(\ref{eqn:theta}) set to be $\theta = \nu$, thus omitting eq.(\ref{eqn:phi}). The diffusion coefficient was taken from \cite{Vern1977}.

  The model's behavior at the central point is qualitatively similar to the spatially homogeneous case \cite{Chizhov2018} and differs  only quantitatively.
  %The observations from multiple points allows us to calculate additional characteristics of the model such as wave's velocity and length. We start from the point's description of the model, continue with a comparison of measurements at multiple points and finally look at the spatial picture of the resulted potassium wave.
  The quasi-periodic spontaneous IDs occur at an interval of about 110-130 seconds (Fig \ref{fig:diff_K_point_center_period}). Each ID is characterized by a high rate of activity lasting about 17 seconds (Fig \ref{fig:diff_K_point_center}) and consisting of short bursts resembling IIDs (Fig \ref{fig:diff_K_point_center}A). 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. \ref{fig:diff_K_board}). $[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.

  %TODO! The similar picture can be seen in \cite{Whalen2018}.

  \begin{figure}
    \centering
    \makebox[\textwidth][c]{\includegraphics[width=1\textwidth]{diff/K_board.png}}
    \caption{Diffusion mechanism (Model 1). Potassium concentration spatial-temporal patterns during generation of two IDs. Two sites S1 and S2 are remote on 2mm.}
    \label{fig:diff_K_board}
  \end{figure}

  The wave profile remains approximately the same during its propagation, as seen from a comparison of $[K]_o$ evolution at two sites (Fig. \ref{fig:diff_K_points}), located in the center and the periphery and marked in Fig. \ref{fig:diff_K_board}. 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. \ref{fig:diff_V_points}. The ID duration is approximately the same at the two sites. The amplitude of the $[K]_o$ pulses varies within 10$\%$ (Fig. \ref{fig:diff_K_points}).
  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. \ref{fig:diff_K_board}). Its velocity is roughly the same, about $0.03mm/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. \ref{fig:diff_K_points}). Still, the velocities are much smaller than the experimental measures.

  %*************************************************************************************************% Модель 2
  \subsection{Spatial aspects in Model 2: Axo-dendritic spread}
  \label{Results_Model2}

  In the consideration of the synaptic mechanism of activity spread, simulations were conducted with the equation for the axo-dendritic propagation of spiking activity, eq.(\ref{eqn:phi}) and eq.(\ref{eqn:theta}) set to be $\theta = \varphi$, thus without eq.(\ref{eqn:K}).

  The model's behavior at the central point is also qualitatively similar to the spatially homogeneous case \cite{Chizhov2018} and the previous scenario (Section \ref{Results_Model1}), with only quantitative differences. The quasi-periodic spontaneous IDs occur at an interval of about 220 seconds (Fig \ref{fig:syn_K_point_center_period}). Each ID is characterized by a high rate of activity lasting about 40 seconds (Fig \ref{fig:syn_K_point_center}) and consisting of short bursts resembling IIDs (Fig \ref{fig:syn_K_point_center}A). 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. \ref{fig:syn_K_board}).
  $[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 \ref{fig:syn_K_point_center}C), 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 \ref{fig:syn_K_board}).
  After reaching its minimum, the potassium concentration slowly increases toward the threshold of another ID initiation after about 381s (Fig. \ref{fig:syn_K_point_center_period},\ref{fig:syn_K_board}).


  %Velocity of the first K wave is about $0.15 mm/s$. The second wave is faster with velocity about $0.21mm/s$.


  \begin{figure}
    \centering
    \includegraphics[width=1.0\textwidth]{syn/x40y40_period.png}
    \caption{Synaptic mechanism  (Model 2). Repeated IDs recorded in the center of the cortical domain(Fig.\ref{fig:syn_K_board}). \textbf{A}, The representative neuron membrane depolarization $U$. \textbf{B}, The total input current $u$. \textbf{C}, The ionic concentrations $[K]_o$ and $[Na]_i$, and the $Na^+/K^+$ pump current $I_{pump}$. \textbf{D}, The somatic firing rate $\nu$ and the synaptic resource $x^D$.}
    %Fig. 6
    \label{fig:syn_K_point_center_period}
  \end{figure}

  \begin{figure}
    \centering
    \includegraphics[width=1.0\textwidth]{syn/x40y40.png}
    \caption{Synaptic mechanism (Model 2). An ictal discharge recorded in the center of the cortical domain (Fig.\ref{fig:syn_K_board}). \textbf{A}, The representative neuron membrane depolarization $U$. \textbf{B}, The total input current $u$. \textbf{C}, The ionic concentrations $[K]_o$ and $[Na]_i$, and the $Na^+/K^+$ pump current $I_{pump}$. \textbf{D}, The somatic firing rate $\nu$ and the synaptic resource $x^D$.}
    %Fig. 7
    \label{fig:syn_K_point_center}
  \end{figure}

  \begin{figure}
    \centering
    \begin{minipage}{0.5\linewidth}
      \centering
      \captionsetup{width=0.8\linewidth}
      \includegraphics[width=1.0\linewidth]{syn/K_points.png}
      \caption{Synaptic mechanism (Model 2). A comparative plot of $[K]_o$ at two different points.}
      \label{fig:syn_K_points}
    \end{minipage}\hfill
    \begin{minipage}{0.5\linewidth}
      \centering
      \captionsetup{width=0.8\linewidth}
      \includegraphics[width=1.0\linewidth]{syn/U_points.png}
      \caption{Synaptic mechanism (Model 2). A comparative plot of $U$ at two different points.}
      \label{fig:syn_V_points}
    \end{minipage}
  \end{figure}

  \begin{figure}
    \centering
    \makebox[\textwidth][c]{\includegraphics[width=1\textwidth]{syn/K_board.png}}
    \caption{Synaptic mechanism (Model 2). Potassium concentration spatial-temporal patterns during generation of two IDs.}
    \label{fig:syn_K_board}
  \end{figure}

  The velocity of the first potassium wave is about $0.11 mm/s$. The second moves with roughly the same velocity (Fig. \ref{fig:syn_K_board}). The speed is proportional to the spatial scale of spatial connections $\lambda$. The speed value 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 \cite{Cammarota2013}, correspondingly.

  The wave profile remains approximately the same during the propagation, as seen from a comparison of $[K]_o$ evolution at two sites (Fig. \ref{fig:syn_K_points}), located in the center and the periphery, and marked in Fig. \ref{fig:syn_K_board}. The moving pulse of $[K]_o$, shown in Fig. \ref{fig:syn_K_points}, 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. \ref{fig:diff_V_points}. 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. \ref{fig:syn_V_cut}A), similar to the analysis performed in \cite{Ma2012} (see Fig.6 there). A diagram of the activity versus the vertical coordinate and time  (Fig. \ref{fig:syn_V_cut}B) demonstrates the effect of synchronization of the voltage bursts between  the center of the ictal discharge generation zone and its periphery. The propagation of the voltage bursts is almost instantaneous in comparison with that of the potassium wave (Fig.\ref{fig:syn_K_board}). This effect of essentially different speeds of the waves of slow and fast variables is quite consistent with the experimental observations from \cite{Ma2012}. 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$).

  \begin{figure}
    \centering
    \makebox[\textwidth][c]{\includegraphics[width=1.2\textwidth]{syn/cut_V_at163s.png}}
    \caption{Synaptic mechanism (Model 2). A.  The membrane potential $V$ at the 163rd second of the simulation. The red bar marks the region of interest (ROI) for spatial-temporal analysis. B. Propagation pattern of the ROI within the interval of 400 ms since the time moment at 163s.}
    \label{fig:syn_V_cut}
  \end{figure}

  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 (Fig. \ref{fig:syn_K_points}, \ref{fig:phi_correlation}). By contrast, the afterdischarges and preictal discharges are much more precisely correlated. Small  firing rate bursts precede an ID (Fig. \ref{fig:phi_correlation}). 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 from the IIDs. 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.

  \begin{figure}
    \centering
    \includegraphics[width=1.0\textwidth]{images/syn/phi_correlation.png}
    \caption{Preictal events at different points are highly correlated, whereas the onsets of ID are delayed. The firing rate $\varphi$ is plotted at different points with the distance 1 mm between them. }
    \label{fig:phi_correlation}
  \end{figure}

  Model 2 shows essentially larger waves than Model 1, as seen from a comparison of the domains with high $[K]_o$ in Fig.\ref{fig:syn_K_board} to those in Fig.\ref{fig:diff_K_board}. 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
  \cite{Kutsy1999, Blume2001, Trevelyan2006}, 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. However, 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\div200$mm/s \cite{Albowitz1993, Schevon2012}. 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 \div 0.8$ preserve the wave-like solution and change the speed in $0.3 \div 20$ times, as obtained in simulations with a narrow computational domain.)

  %\subsection{Velocity}
  The crucial role of the depolarization provided by the elevation of $[K]_o$ in the wavefront propagation is confirmed by the increase of the speed with the increase of $g_{K,leak}$. The potassium-mediated positive feedback depends on the sensitivity of the membrane polarization to $[K]_o$, evaluated by $g_{K,leak}$. Hence the stronger the feedback, the faster the ictal wavefront.

  %*************************************************************************************************
  \subsection{Model 2: Domain with a lesion}
  \label{Results_Model2_lesion}

  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 (Fig.\ref{fig:Model2_lesion}). The activity cannot propagate across the lesion. Instead, it passes around the lesion and spreads through the entire area behind.

  \begin{figure}
    \centering
    \makebox[\textwidth][c]{\includegraphics[width=1\textwidth]{images/Model2_lesion_K_board_2019-07-12_10_14.png}}
    \caption{Model 2: domain with a lesion. Potassium concentration spatial-temporal patterns during generation of two IDs.}
    \label{fig:Model2_lesion}
  \end{figure}

  %*************************************************************************************************% Комбинированная модель
  \subsection{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 equations (\ref{eqn:K}-\ref{eqn:U_reset}). As suggested, simulations with the combined model have shown results quite similar to those from Model 2, which were described in the subsections \ref{Results_Model2} and \ref{Results_Model2_lesion}. 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.\ref{fig:Model2_lesion} 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.\ref{fig:Model2_lesion} and \ref{fig:Model3_lesion}. Consequently, the potassium diffusion does not affect the spread of IDs.

  \begin{figure}
    \centering
    \makebox[\textwidth][c]{\includegraphics[width=1\textwidth]{images/Model3_lesion_K_board_2019-07-12_12_02.png}}
    \caption{Model 3: domain with a lesion. Potassium concentration spatial-temporal patterns during generation of a single ID.}
    \label{fig:Model3_lesion}
  \end{figure}




  %********************************************************************
  \section*{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 \cite{Martinet2017}, 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 \cite{Wong1990, Trevelyan2006, Trevelyan2007} as well as in vivo in animal models \cite{Wenzel2017, Rossi2017} and in human patients \cite{Schevon2012, Schevon2019}. 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.

  {\sl Role of feedforward inhibition in seizure propagation}

  \noindent Current experimental studies of seizure propagation reveal the role of feedforward inhibition  \cite{Trevelyan2006, Trevelyan2007, Schevon2012, Wenzel2017, Schevon2019, Cammarota2013}; however, a certain mechanism is still debatable \cite{Schevon2019}. 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 \cite{Trevelyan2006}. 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 \cite{Trevelyan2006}. 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 \cite{Albowitz1993}.

  {\sl Role of potassium-mediated positive feedback in seizure propagation}

  \noindent 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 \cite{Weissinger2000}.

  {\sl Speed of preictal bursts and afterdischarges}

  \noindent 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 \cite{Schevon2012, Schevon2019} (for instance, see Fig. 4a in \cite{Schevon2012}). By contrast, the afterdischarges and preictal discharges are much more precisely correlated.
  In electrophysiological recordings by multiple electrodes \cite{Schevon2012}, 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 \cite{Schevon2012}. 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 \cite{Ma2012} or intracellular signals \cite{Trevelyan2007b}. 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 \cite{Schevon2012,Schevon2019}.

  {\sl A lesion prevents the ID propagation }

  \noindent 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 \cite{Lian2001}, where epileptiform discharges were found to be synchronized on both sides of a complete cut. Noting that the observed discharges were not ictal dischages but spontaneously repeating short epileptiform events, we suggest that their synchronization occurs on time scales longer than that of a single ID propagation. In the case of repeating events, the correlation of activity in two parts of the slice may be mediated by the potassium diffusion, which is otherwise too slow to affect an ID.

  {\sl Analogy to spreading depression waves}

  \noindent 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) \cite{Whalen2018}. 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 \cite{Ayata2015, Whalen2018}. As is known \cite{Ayata2015}, 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.


  %********************************************************************
  \section*{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.

  %\section*{Supporting information}

  \section*{Acknowledgments}
  This work was supported by the Russian Science Foundation (project 16-15-10201).

  \nolinenumbers


  \bibliography{article}


\end{document}