Mathematical model of Na-K-Cl homeostasis in ictal and interictal discharges

Despite big experimental data on the phenomena and mechanisms of the generation of ictal and interictal discharges (IDs and IIDs), mathematical models that can describe the synaptic interactions of neurons and the ionic dynamics in biophysical detail are not well-established. Based on experimental recordings of combined hippocampal-entorhinal cortex slices from rats in a high-potassium and a low-magnesium solution containing 4-aminopyridine as well as previous observations of similar experimental models, this type of mathematical model has been developed. The model describes neuronal excitation through the application of the conductance-based refractory density approach for three neuronal populations: two populations of glutamatergic neurons with hyperpolarizing and depolarizing GABAergic synapses and one GABAergic population. The ionic dynamics account for the contributions of voltage-gated and synaptic channels, active and passive transporters, and diffusion. The relatively slow dynamics of potassium, chloride, and sodium ion concentrations determine the transitions from pure GABAergic IIDs to IDs and GABA-glutamatergic IIDs. The model reproduces different types of IIDs, including those initiated by interneurons; repetitive IDs; tonic and bursting modes of an ID composed of clustered IID-like events. The simulations revealed contributions from different ionic channels to the ion concentration dynamics before and during ID generation. The proposed model is a step forward to an optimal mathematical description of the mechanisms of epileptic discharges.

The model of synaptically interacting neuronal populations is based on our previous study [1].We consider two excitatory and one inhibitory neuronal populations, denoted by indexes E1, E2 and I, correspondingly.The types of synapses are denoted by the types of mediator (AMPA, GABA or NMDA), and pre-and postsynaptic neurons.
A mathematical description of each single population is based on the probability density approach [2], namely, the conductance-based refractory density (CBRD) approach [3].The approach considers a population of an infinite number of Hodgkin-Huxley-like neurons receiving both a common input and an individual for each neuron noise.In any arbitrary case of transient or steady-state stimulation the firing rate of such population can be quite precisely and computationally effectively calculated by solving a system of equations in partial derivatives, 1-d transport equations.The equations govern an evolution of neuronal states distributed in a phase space of the time elapsed since last spikes, * t .They contain the Hodgkin-Huxley equations for the membrane voltage and gating variables, parameterized by * t , as well as the equation for the neuronal density in * t -space, ) , ( * t t E ρ , where the index E substitutes for E1 or E2.The output characteristic of the population's activity is the firing rate , which is equal to E ρ in the state of a spike, . The equations written below describe an excitatory population of adaptive regular spiking pyramidal cells according to our previous works [3,4].Basic  .In comparison with one-compartment model, the extra parameters is the ratio of dendritic to somatic conductances γ and the dendritic length. The inhibitory synapses are located at soma, contributing into the somatic synaptic current soma I , whereas the excitatory synapses are at dendrites, determining the dendritic synaptic current dendr I .Due to the construction of the 2-compartment model [5], both type synaptic conductances are imposed to be somatic, in spite of the localization, in order to be compared with experimental whole-cell somatic registrations.Approximations of voltage-gated ionic currents are based on the CA1 pyramidal cell model from [6], where instead of full description of calcium dynamics and calciumdependent potassium currents a cumulative after-spike hyperpolarization (AHP) current, that provides an effect of slow spike timing adaptation [7].Parameterized by * t , the governing equations are as follows: PDF Created with deskPDF PDF Writer -Trial :: http://www.docudesk.com l is the square ratio of the dendritic length to the characteristic length.The somatic and dendritic synaptic currents soma I and dendr I are calculated as where the differential operator represents the solution of the reverse problem of dendritic current estimation from somatically registered-like conductances [5].The synaptic conductance kinetics is estimated from somatic responses to stimulation of presynaptic neuronal population, thus it implicitly accounts not only the kinetics of synaptic channels but also the dendritic and axonal propagation delays.For the dendritic compartment, the differential operator sharpens the transient effect of the channels, thus providing better agreement between somatic postsynaptic currents and potentials.This sharpening affects only glutamatergic channels located on the dendritic compartment.

Hazard function
The source term in the eq.(A1) is the hazard function H which is defined as the probability for a single neuron to generate a spike, if known actual neuron state variables.The approximation of the hazard function H has been obtained for the case of white noise [3] and color noise [4] where T is the membrane potential relative to the threshold, scaled by the noise amplitude V σ which increases with the synaptic conductance as The term A is the hazard for a neuron to cross the threshold because of noise, derived analytically [3] and approximated by exponential and polynomial for convenience; B is the hazard for a neuron to fire because of depolarization due to deterministic drive, i.e. the hazard due to drift in the voltage phase space.Note that the H -function is independent of the basic neuron model and does not contain any free parameters or functions for fitting to any particular case.Thus, H -function is the same for excitatory and inhibitory populations.

Voltage-dependent channels
PDF Created with deskPDF PDF Writer -Trial :: http://www.docudesk.com The set of ionic currents includes the voltage-dependent potassium currents DR I and A I responsible for spike repolarization, the slow potassium current M I that contributes to spike frequency adaptation and the potassium current AHP I , implicitly dependent on calcium dynamics and contributing to spike frequency adaptation.Approximating formulas for the currents Na I , DR I , A I and M I are taken from [6]; the approximation for AHP I is given in [7].The voltage-dependent potassium current DR I : The voltage-dependent potassium current A I : PDF Created with deskPDF PDF Writer -Trial :: http://www.docudesk.com The adaptation current AHP I :

Boundary conditions.
According to the conservation of the number of neurons in a population, the firing rate is calculated as a sink of neurons from their state (A21) The reset values for the fast gating variables in eqs.(A18,A19) were obtained with the basic single neuron model.With a rather arbitrary input providing a spike, these values were measured at the moment of a voltage maximum at the spike.The reset level for each slow conductance in the CBRD model was calculated as a sum of its value at a peak of spike-release distribution in the * t -space and an increment at spike: Here 0 tot g is the total somatic conductance at rest, and syn g is the total synaptic conductance; S is the membrane area.The dependence of is taken from a full single neuron model [3], allowing to take into account the effect of sodium channel inactivation on the threshold dynamics [8].V σ is the noise amplitude meaning the dispersion of individual neuron's voltage fluctuations in a stationary state.Its scaling with syn g approximately reflects the fact of the synaptic noise increase with the increase of mean synaptic drive [9].Stochastic input to E2-neurons noise I was modeled as Ornstein-Uhlenbeck process with the time correlation 10 ms and the dispersion 20 pA for the ID-regime simulation and 40 pA for the IID-regime simulation.
The equations for the input synaptic conductances are given below, as well as the values of the reversal potentials.When calculating the dynamics of a neural population, the integration of eqs.(A2,A3,A5-A16) determines the evolution of the distribution of voltage E U across * t .Then, the effect of crossing the threshold and the diffusion due to noise are taken into account by H -function, eq.(A4), substituted into the equation for neuronal density, eq.(A1).The integral eq.(A17) results in the output firing rate

CBRD-approach for a population of interneurons
The model for the fast spiking single-compartment interneurons is the reduction of the model of the adaptive regular spiking two-compartment neurons: )

Lognormal distribution of synaptic weights within each population
The CBRD-approach is generalized to the case of lognormal distribution of synaptic weights within each j -population [10].In this case, instead of equal total synaptic current, neurons receive lognormally distributed current.For the current scaled by its mean across the distribution, x , the distribution is The membrane potential of neurons parameterized with x , j x U , can be found as .The output firing rate is defined as In numerical simulations, we set the parameter of the lognormal distribution 0.75 = LN σ and discretized the x -space by 15 intervals.Mean somatic and dendritic membrane potentials are calculated as follows:

Connections
The synaptic conductances are described with the second-order differential equations [11] with introduced synaptic plasticity factors  mS/cm 2 ,

Representative neurons
PDF Created with deskPDF PDF Writer -Trial :: http://www.docudesk.com Representative neurons of each of the populations were modeled by the basic single neuron model with the same synaptic inputs as for the populations.The model is described by the equations for the membrane voltage, eq.(A2,A3,A5-A16) for E1 and E2-populations, and eq.(A25,A5-A10) for I-population, where the sum of partial derivatives were substituted by the total derivative in time t, and the sodium current was explicitly present in the right-hand part of eq.(A2) and eq.(A25).The sodium current dependent on voltage V was approximated by the 4-state Markov model [6]:

LFP model
Two compartment model of principal neurons (E1-population) allows calculation of the local field potential (LFP) originating from dipole-like configuration of membrane currents [5] ( )

2 .
per area of homogeneously distributed neurons p , the specific intracellular resistivity i r , the conductivity of extracellular medium σ , the maximum sodium current at a spike max Na I , the sipke duration AP τ , and the dendritic length L .These coefficients were estimated as PDF Created with deskPDF PDF Writer -Trial :: http://www.docudesk.com

+∞
The increment values for the slow gating variables in eqs.(A22-A23) were also measured at a single spike of the single neuron model.
p t * is such that PDF Created with deskPDF PDF Writer -Trial :: http://www.docudesk.com ,Here i ϕ is the presynaptic firing rate.In neglect of spatial propagation and temporal delays the presynaptic firing rate is equivalent to the somatic firing rate, i.e. τ are the rise and decay time constants.We imply that the synaptic time constants are estimated from the somatic responses to the stimulation of a presynaptic neuronal population, thus these time constants characterize not only synaptic channel kinetics but the dendritic and axonal propagation delays as well.The time scale . The LFP signal measured at the somatic depth,