Slow ion concentration oscillations and multiple states in neuron-glia interaction-insights gained from reduced mathematical models

Abstract

When potassium in the extracellular space separating neurons and glia reaches sufficient levels, neurons may fire spontaneous action potentials or even become inactivated due to membrane depolarisation, which, in turn, may lead to increased extracellular potassium levels. Under certain circumstances, this chain of events may trigger periodic bursts of neuronal activity. In the present study, reduced neuron-glia models are applied to explore the relationship between bursting behaviour and ion concentration dynamics. These reduced models are built based on a previously developed neuron-glia model, in which channel-mediated neuronal sodium and potassium currents are replaced by a function of neuronal sodium and extracellular potassium concentrations. Simulated dynamics of the resulting two reduced models display features that are qualitatively similar to those of the existing neuron-glia model. Bifurcation analyses of the reduced models show rich and interesting dynamics that include the existence of Hopf bifurcations between which the models exhibit slow ion concentration oscillations for a wide range of parameter values. The study demonstrates that even very simple models can provide insights of possible relevance to complex phenomena.

Keywords: bifurcation analysis; bursting; ion concentration dynamics; mathematical model; neuron–glia interplay.

FIGURE 1

Properties of the neuron–ECS–glia system and of neuron and neuron–glia models. (A) Schematic representation of the neuron–glia system with channels, pumps, and cotransporters in the neuron and glia membranes (depicted sodium and potassium channels in the neuron membrane comprise several channels). (B) Bifurcation diagram of the model of Kager et al. (2000) using [K⁺]_o as the bifurcation parameter [(Na⁺)_n is fixed at 10 mM, and model parameter values are given in Supplementary Material]. “HB” and “HC” indicate locations of a Hopf bifurcation and a heteroclinic bifurcation, respectively. (C) Two-parameter bifurcation diagram of the model of Kager et al. (2000) using [K⁺]_o and [Na⁺]_n as bifurcation parameters showing how curves that correspond to the HB and HC bifurcations in (B) separate the parameter plane into regions corresponding to the resting state (RS), spontaneous discharge (SD), and depolarisation block (DB) behaviour. Curves displayed are solution orbits corresponding to the dynamics displayed in (D) and (E) using the same colour coding. (D)–(E): Dynamics of neuronal membrane potential (D) and of [K⁺]_o (E) obtained from simulations using the full neuron–glia model (parameter values are given in Supplementary Material). The glial sodium–potassium pump rate $J_{\text{NaKATPase,max}}^{\left(\text{g}\right)}$ is multiplied by factors of 0.8 (top panels in D and E), 0.65 (middle), and 0.55 (bottom) to generate different types of behaviour.

FIGURE 2

Calculated and fitted sodium transmembrane neuronal current. (A) Magnitude of the sodium current into the neuron as a function of [K⁺]_o and [Na⁺]_n computed from the neuronal model of Kager et al. (2000). (B) Best fit to the current plotted in (A) to the function given in Eq. 7a obtained by the parameters y ₁ = 0.367, y ₂ = 25.52, y ₃ = 6.39, y ₄ = 35.38, y ₅ = 1.10, y ₆ = 2.51, y ₇ = 3.33, y ₈ = 16.38, and y ₉ = 137.9 (the unit for y ₁ and y ₃ is pmol ms⁻¹ cm⁻², for y ₂: pmol ms⁻¹ cm⁻² mM⁻¹, for y ₄, y ₅, y ₇, and y ₉: mM. y ₆ and y ₈ are unit-free). (C) Contour plot of the current in (A) as a function of [K⁺]_o and [Na⁺]_n with indications of the expected behaviour in three regions of the ([K⁺]_o,[Na⁺]_n)-plane.

FIGURE 3

Comparison of [K⁺]_o dynamics obtained in simulations using the full and reduced models. Columns (A–F) display [K⁺]_o dynamics for increasing values of the sodium–potassium pump rate obtained by multiplication by a factor f _NaK in the range 0.58–0.74, indicated at the top of each column. Full model dynamics are displayed in the top panel, RM1 in the middle, and RM2 in the bottom panel.

FIGURE 4

Bifurcation and phase plane analysis for RM1 and RM2. (A1-A2): Bifurcation diagram of RM1 with the sodium–potassium pump rate reduction factor f _NaK as the bifurcation parameter. (A1) shows the full diagram, and (A2) shows details of the diagram within the rectangle depicted in (A1). The black curve indicates stable steady states; red curve, unstable steady states; black filled circles, stable oscillatory solutions; and red filled circles, unstable oscillatory solutions. “HB” and “HC” indicate locations of Hopf and heteroclinic bifurcations, respectively. (B1–B2): Same as (A1–A2) for RM2. (C1–C4): Solution orbits and phase plane analysis for RM2 in the $({[{\mathrm{K}}^{+}]}_{\mathrm{o}},{[\mathrm{N}{\mathrm{a}}^{+}]}_{\mathrm{n}})$ -plane. In all plots, red and blue curves indicate nullclines of $d{N}_{{\text{Na}}^{+},\text{n}}/dt$ and $d{N}_{{\text{K}}^{+},\text{o}}/dt$ , respectively, and black curves are solution orbits. The depicted dynamics correspond to different values of f _NaK; (C1) f _NaK = 0.74, (C2) f _NaK = 0.64, (C3) f _NaK = 0.62, and (C4) f _NaK = 0.58.

FIGURE 5

Bifurcation analysis and [K⁺]_o dynamics showing the effect of pump rate and NKCC1. (A) Two-parameter bifurcation diagram showing how dynamics of the RM1 depends on f _NaK and f _NKCC1. The three regions indicated by RS, SD, and DB correspond to regions of the (f _NaK, f _NKCC1)–plane where the resting state, spontaneous discharge, and depolarisation block behaviour are observed, respectively. “A1”–“A4” in the figure indicate the locations in the parameter plane associated with simulations whose results are shown in (A1–A4). (B) Same as (A) for RM2. “B1”–“B4” in the figure indicate the locations in the parameter plane associated with simulations whose results are shown in (B1–B4).

FIGURE 6

Bifurcation analysis and [K⁺]_o dynamics showing the effect of pump rate and [K⁺]_g. (A) Two-parameter bifurcation diagram showing how dynamics of the RM1 depends on f _NaK and [K⁺]_g. The three regions indicated by RS, SD, and DB correspond to regions of the $(f_{\text{NaK}},{[{\text{K}}^{+}]}_{\text{g}})$ -plane where the resting state, spontaneous discharge, and depolarisation block behaviour, respectively, are observed. “A1”–“A4” in the figure indicate the locations in the parameter plane associated with simulations whose results are shown in (A1–A4). (B) Same as (A) for RM2. “B1”–“B4” in the figure indicate the locations in the parameter plane associated with simulations whose results are shown in (B1–B4).

FIGURE 7

Explanation for how burst cycles are maintained. (A) [K⁺]_o (blue) and [Na⁺]_n (red) dynamics during one cycle of bursting obtained by numerically solving the RM2 model Eqs 10a, 10b using the default parameter set except the sodium–potassium pump rate, which is multiplied by 0.70 to generate periodic bursting solutions. The times of entering and exiting from the SD region are indicated (black dashed vertical lines) and labelled “1” and “2,” respectively. (B) Magnitude of glial transmembrane K⁺ flux (blue) and of neuronal transmembrane K⁺ and Na⁺ fluxes (red) during one cycle. The black vertical dashed lines and the numbering have the same meaning as in (A). (C) Solution orbit during one cycle (blue) in the $({[{\mathrm{K}}^{+}]}_{\mathrm{o}},{[\mathrm{N}{\mathrm{a}}^{+}]}_{\mathrm{n}})$ -plane and the line that separates the RS and SD regions (black). Labels “1” and “2” refer to the times indicated in (A) and (B), respectively.