Mixed-mode oscillations in pyramidal neurons under antiepileptic drug conditions.

Subthreshold oscillations in combination with large-amplitude oscillations generate mixed-mode oscillations (MMOs), which mediate various spatial and temporal cognition and memory processes and behavioral motor tasks. Although many studies have shown that canard theory is a reliable method to investigate the properties underlying the MMOs phenomena, the relationship between the results obtained by applying canard theory and conductance-based models of neurons and their electrophysiological mechanisms are still not well understood. The goal of this study was to apply canard theory to the conductance-based model of pyramidal neurons in layer V of the Entorhinal Cortex to investigate the properties of MMOs under antiepileptic drug conditions (i.e., when persistent sodium current is inhibited). We investigated not only the mathematical properties of MMOs in these neurons, but also the electrophysiological mechanisms that shape spike clustering. Our results show that pyramidal neurons can display two types of MMOs and the magnitude of the slow potassium current determines whether MMOs of type I or type II would emerge. Our results also indicate that slow potassium currents with large time constant have significant impact on generating the MMOs, as opposed to fast inward currents. Our results provide complete characterization of the subthreshold activities in MMOs in pyramidal neurons and provide explanation to experimental studies that showed MMOs of type I or type II in pyramidal neurons under antiepileptic drug conditions.

Introduction

The entorhinal cortex (EC) plays a pivotal role in the generation of theta and gamma rhythms during rapid-eye-movement sleep, memory function, learning, and exploring behaviors in awake animals [1-5]. The EC is a vital part of the temporal lobe that acts as an intermediary in the neocortical-hippocampal network. Incoming information to the hippocampus is mostly transmitted by stellate cells through the prefrontal pathway in layer II of the EC; whereas outgoing information from the hippocampus is transmitted through layer V of the EC.

Morphological studies of layer V have shown that the majority of cells in this layer are pyramidal neurons, which display a robust, rhythmic subthreshold activity similar to stellate cells [2]. It has been also shown that pyramidal neurons in layer V of the EC might contribute to the generation of spontaneous seizures after status epilepticus, to immobility in the awake animal, and to slow-wave sleep [6-9]. Importantly, spontaneous seizures have been shown to be affected considerably by subthreshold oscillatory activities of these neurons [10-12]. Studies of antiepileptic drugs, in which persistent sodium current (I) is inhibited with drugs such as riluzole and phenyton, have interestingly showed that pyramidal neurons still exhibit subthreshold oscillations [4, 11–13], which probably contribute to the abnormal firing patterns of pyramidal neurons under antiepileptic drug conditions. Therefore, the main goal of this study is to investigate the properties of STOs in pyramidal neurons when I is inhibited, and to examine the relationship between these subthreshold activities and the firing patterns in pyramidal neurons.

These STOs are small-amplitude oscillations (< 10 mV), and are intrinsic neuronal phenomena that persist during synaptic transmission block [14-16]. STOs were first reported in myelinated nerves [17], and have since been observed in many neuron types, such as inferior olive [18], squid giant axon [19], EC [2, 20], hippocampal CA1 [21], pyloric dilator [22], neocortex [23], thalamus [24], and spinal motoneurons [25, 26]. Many experimental studies have suggested that STOs influence a variety of neural behaviors, such as the size of bursting (or spike clustering) [27], synaptic plasticity [21–25, 27, 28], rhythmic activities at the network level [29, 30], spike discharge panel and spike clustering [2, 21–25, 27–29], perception and memory [31, 32], synchronous activities [33], and motor coordination [34]. The properties of STOs are usually investigated by injecting a chirp signal (i.e., a sinusoidal current with increasing frequencies) into neurons, while its synaptic receptors are blocked [29]. Under this condition, the membrane potential of resonant neurons displays the maximum response (or STO) for stimuli with frequency matching the neurons’ resonant frequency.

The resonance frequencies are usually examined through linearized models of resonant neurons [35, 36]. One of the drawbacks of linearized models, however, is the disregard of the firing region of neural oscillations by removing the Hodgkin-Huxley sodium and potassium currents that mediate the action potential [36]. This limitation has impeded efforts to examine the effects of STO phenomena on the firing patterns of neurons, thus inhibited a better understanding of how such phenomena might impact neural encoding. Therefore, the focus of this study is on the nonlinear properties of conductance-based (CB) models, in order to improve our understanding of the role of STO phenomenon in shaping the firing pattern, and also to investigate the electrophysiological mechanisms underlying those firing patterns of pyramidal neurons in layer V of the EC.

It is now widely accepted that the presence of a potassium channels with slow kinetics induces the STO phenomena in many neurons. However, the effect of potassium channels in regulating these phenomena in epilepsy remains poorly understood. There are over 70 potassium channel genes in the mammalian body, however only few of them have been linked to diseases. Among these genes, KCNQ1 and KCNE1 underlie diseases including epilepsy and cardiac arrhythmias. KCNQ1 in conjunction with KCNE1 encode slow non-inactivating potassium current (I) [37]. The mutations in this current leads to altered resting potential, disruption of neural firing pattern, increased neurotransmitter release and calcium influx, which in turns could generate hyperexcitability and epileptic phenotype. On the other hand, recent studies showed that mutation in I could greatly impact inherited epilepsies because of its role in sustaining epileptic discharges and long membrane depolarization [12]. Therefore, I is the target of many antiepileptic drugs such as riluzole, phenytoin and topiramate ([12] and references therein). Riluzole can inhibit I, which accounts for about 1% of the sodium current, without altering the transient sodium current (I) thereby sustain normal neuronal excitability and information processing. Inhibition of I causes the loss of neuronal repetitive firing and the emergence of STOs and spikes (action potentials), termed mixed-mode oscillations [38]. In this study, therefore, we used a model of pyramidal neurons in layer V of EC in which the persistent sodium channel is inhibited to simulate the effect of riluzole in order to investigate the separate effect of I on the STO phenomenon.

Many experimental and theoretical studies have shown that the combination of STO and spiking phenomena in resonant neurons result in complex oscillatory activity, called mixed-mode oscillations (MMOs) [2, 14, 20, 39–43]. The concept of MMOs was first proposed in the Belousove–Zhabotinsky reaction [44], and has been subsequently reported in several models and experiments in biological and chemical systems [20, 45–47]. In this study, we are dealing with two types of MMOs, the first of which we call type I MMOs. In this type, the depolarization of membrane potential results in STO phenomena around the steady-state potential; then further depolarization results in action potentials (spikes) at the peak of STOs generating MMOs [48]. The STO phenomenon in MMOs mainly play a timing role and control the spike clustering (groups of action potentials separated by silent STO periods) [49-53]. By considering this role for STOs, the authors suggest that there is another type of MMO in which small-amplitude oscillations (presumably STOs) emerge after the spike and during the repolarization period. In literature, this type of MMO is known as pseudo-plateau bursting [39, 40, 54]. During these small-amplitude oscillations, neurons cannot fire action potentials. The ability to delay the generation of the next action potential is similar to the timing role of STO phenomena in type I MMOs. Therefore, we consider the combination of spike followed by STO-type oscillations in these neurons as type II MMOs.

From the mathematical point of view, there are several methods for analysing the MMO phenomena in nonlinear systems with multiple timescales, such as neurons. These methods include break-up/loss of stability of a Shilnikov homoclinic orbit [46], break-up of an invariant torus [55], and a subcritical Hopf bifurcation (HB) with an appropriate return mechanism [56]. However, these mechanisms cannot completely explain all MMO features in a nonlinear system, such as in neurons [57]. Due to these limitations, the present study employs canard theory, which can examine transient dynamics in multiple timescale systems [58].

The classic canards are generic and were initially studied by Benoit et al. through nonstandard analysis [59]. These canards describe the fast transition from small-amplitude STOs to large-amplitude oscillation (i.e., spikes) during the variation of a parameter. Recently, Wechselberger et al. reported that a class of three-dimensional canards, called type-I folded nodes, produces small-amplitude oscillations in MMOs [60]. Afterwards, several studies applied canard theory to the CB model of neurons to predict the number of STOs per spike in MMOs [61-65]. However, further research has yet to explore how the mathematical results of such studies relate to the neuronal electrophysiological mechanisms.

The first goal of this study is therefore to apply canard theory to the CB model of pyramidal neurons in layer V of the EC to investigate the properties of their MMO behaviors. It has been shown that these neurons have non-inactivating potassium current (I) and persistent sodium current (I). These currents mediate subthreshold activities with amplitudes in the 3–5 mV range and frequencies in the 5–15 Hz range [2, 33]. It has also been shown that these neurons can display MMOs [2]. Based on this finding, Jalics et al. developed a compartmental model derived from that of Acker et al., in which they examined the number of STOs per spike in MMOs. They used three-time-scale, singular perturbation stability analysis to show that these neurons can display both MMOs and a seemingly chaotic firing pattern [33, 48]. We mainly investigate the range of parameters in which the STO phenomena shape the spike cluster. This information is expected to improve understanding of the effects which ionic currents in pyramidal neurons have on information encoding in the EC [29].

The second goal of this study is to bridge the gap between the theoretical results obtained by the canard theory and the electrophysiological mechanisms that could potentially underlie the mathematical predictions. To this end, we first identify the ranges of maximum conductances of ionic currents wherein the MMOs exist then discuss the electrophysiological behaviors resulting from these ionic currents in those ranges. Various factors are considered, such as the input current, the evolution of membrane potential, the states of activation and inactivation functions, the time constants of membrane potential and ionic currents, the magnitude of maximum conductances, and the dynamics of ionic currents. Accordingly, this paper is organized as follows: Section 2 presents the full and reduced compartmental biophysical single-neuron models of pyramidal cells. Section 3 analyses how varying specific parameters impacts MMOs by using canard theory. Moreover, Section 3 discusses how the results obtained by applying canard theory to the reduced model can be interpreted in terms of electrophysiological mechanisms. It should be noted that, in this study, we examine the entire range of each parameter within which it is physiologically feasible for MMOs to exist; but we do not examine the number of STOs per spike in MMOs [53, 54, 56, 60].

In this study, our overarching goal was to apply canard theory to the CB model of pyramidal neurons in layer V of EC to investigate the electrophysiological mechanisms underlying information coding efflux from hippocampus to neocortex. We develop a correlated mathematical and electrophysiological understanding of the impact of subthreshold activities on the firing behaviour (i.e. spike clustering) of pyramidal neurons. Our findings improve our understanding of the contributions of ionic currents to shaping spike clustering.

Method

Conductance-based model

In this study, we analyzed a single-compartment biophysical model of pyramidal cells in layer V of the EC. This model, which was introduced by Acker et al. and modified by Jalics et al. [33, 48], generates MMOs over a range of parameter values. The model consists of slow non-inactivating potassium (I) and persistent sodium (I) currents, with standard Hodgkin–Huxley potassium (I), sodium (I), and leak currents (I) [66]. As shown in Fig 1, the current balance equation is, where I is the applied bias (DC) current (μA/cm2), C is the membrane capacitance (μF/cm2), V is the membrane potential (mV), I = g m (V−E), I = g m (V–E), I = g n4 (V−E), I = g m3 h (V−E), and I = g (V–E). E and G (x = ks, k, nap, na, L) are the Nernst (reversal) potential (mV) and the maximal conductance (μS/cm2), respectively. All gating variables follow a first-order differential equation with the form of, where,

Fig 1

Conductance-based model of pyramidal neurons.

The model composed of a capacitor and five parallel resistors. The constant resistor represents the passive leakage current, g. Four other resistors correspond to the active conductances, g, g, g, g. The g and g are conductances of Hodgkin-Huxley model for generating fast action potentials. E, E and E are the reversal potentials of potassium, sodium and leak currents respectively.

The definition of α, β and the parameter values are defined in the S1 File. x(V) denotes the steady-state function of gating variables, and τ(V) represents the time constant functions. Eqs 1 and 2 define the full six-dimensional model of a pyramidal cell in the layer V of EC.

The behavior of this model varies in three particular regions: (A) at I < 0.8802 (μA/cm2), damped oscillation converges to a stable rest state; (B) at 0.8802 < I < 0.99406 (μA/cm2), a bi-stability region exists and includes both a stable resting state and spike solutions; and (C) at I > 0.99406 (μA/cm2), stable periodic solutions exist. States in (A) and (B) are separated by unstable, small-amplitude oscillation caused by a subcritical HB. Moving from the first region to the second region, the trajectory follows a repelling unstable manifold, which is reminiscent of two-dimensional canard phenomenon. Moreover, the periodic solution in region (C) is related to the spiking behavior of neurons.

Reduced model

The reduced model of pyramidal neurons in layer V of the EC was introduced by Jalics et al. [48], as follows:

They used steady-state (in)activation functions and time constant functions to reduce the full model. Fig 2 illustrates these functions. Given that the bifurcation diagram of the reduced model is similar to that of the full six-dimensional model (not shown), Jalics et al. [48] concluded that the behaviour dynamics of the full and reduced models are comparable. Generally, the main reduction principle is elimination of the Hodgkin–Huxley spiking currents [67, 68] because Hodgkin–Huxley currents do not affect the membrane potential in the subthreshold range [68]. However, the elimination of Hodgkin–Huxley currents leads to a qualitative change in neuronal dynamics. Therefore, we reduced the full model using a different approach as described below in order to avoid changing the neuronal dynamics.

Fig 2

The evolutions of dynamics of the pyramidal cells in layer V of EC.

a) Steady-state (in)activation functions. b) Voltage-dependent time constant functions (τ = 90 msec).

In the reduced pyramidal neuron model the persistent sodium current (I was inhibited, which causes the spiking threshold of the neuron to increase. Therefore, a higher input current I was needed in order to activate the cell and examine its behaviour [48]. This increase also helps to sustain both firing rate and the level of current while the dynamics are eliminated. Importantly, our results show that pyramidal neuron model with inhibited I exhibits two types of MMOs.

In the reduced model, a single equilibrium loses its stability by a subcritical HB at I = 16.93 (μA/cm2), and the neuronal behavior switches among four regions (Table 1).

Table 1

Neuronal oscillation types for different values of input current.

(μA/cm2)	DO*	STO	MMO	Spiking
Iinp< 16.90	Yes	─	─	─
16.90 < Iinp< 16.93	─	Yes	─	─
16.93< Iinp< 18.26	─	─	Yes	─
Iinp > 18.26	─	─	─	Yes

*DO stands for damping oscillations.

We investigate the mechanisms underlying this behavior in the reduced model with I = 0 and I ≈ 17 (μA/cm2), unless otherwise mentioned. The variables V, m, and n evolve on distinct timescales. For the membrane potential V, the time constant is calculated by τ = C / g, where g is obtained by,

The time constant of m is τ = 90 (msec), whereas that of n can be obtained from the equation of the activation variable (Eq 2). The results show that the time constant of V changes faster than those of m and n, thereby creating one fast and two slow variables in the system (this relationship is further discussed in the following section). Moreover, setting C to the smaller values (i.e. C → 0) decreases the time constant of V. Therefore, C is considered as the parameter of the dimensionless singular perturbation problem in this model [63]. Herein, at C ≈ 0, m and n are regarded as slow variables and V is the only fast variable in the model. By taking advantage of these assumptions, we show that the MMO patterns in the reduced model emerge through a canard mechanism, as a result of the existence of a folded node singularity [57, 60, 69]. The existence of these folded node singularities guarantees the existence of STO phenomena in these neurons [60].

Geometric singular perturbation theory (GSPT)

The reduced model (Eqs 4–6) is a singularly perturbed system and consists of one fast variable (V) and two slow variables (m and n). We apply GSPT to analyze the model and determine whether its MMOs are canard-induced oscillations [70, 71]. GSPT separates the analysis of the model into two low-dimensional limiting problems by reducing the singular perturbation parameter C to 0.

The first system is a one-dimensional system at the fast timescale (τ1 = τ / ε) and called the layer system, which describes the fast dynamics and is given by,

The second system is a two-dimensional system at the slow timescale (τ) and called the reduced system, which describes the evolution of slow variables (m and n) and is given by (4) with C = 0.

The trajectory of the reduced system lies on the critical manifold S defined by S = {(V, m, n) ∈ ℝ3: f (V, m, n) = 0}. Furthermore, the critical manifold equation can be solved explicitly for one of the slow variables (m and n). Therefore, the critical manifold can be given as,

Fig 3 illustrates that at I = 17.5 (μA/cm2), the critical manifold S is a folded surface, which characterizes a large number of physiological systems with multiple timescales, such as neurons [57]. This property allows these systems to change from one state to another, such as between the spiking and subthreshold regions. In Fig 3, the critical manifold S is composed of three surfaces divided by two fold-curves (L and L). The middle surface (light blue) is the space of repelling points (∂f / ∂V > 0), and the upper (red) and lower (dark blue) surfaces are the spaces of attracting points (∂f / ∂V < 0). Fold-curves are obtained by,

Fig 3

The critical manifold for reduced model of EC neurons.

The critical manifold with fold-curves (L+ and L-) and their projections [P(L+) and P(L-)]. L+ and L- are the upper and lower fold-curves respectively. P(L+) and P(L-) are the projections of L+ and L- onto the lower and upper sheets of the critical manifold respectively. The upper sheet (red color) and lower sheet (dark blue color) are attracting sheets. The middle sheet (light blue color) is repelling sheet.

This equation yields two values for the membrane potential (V) and two equations for m in the form of m = m(n). The fold-curves L and L are projected vertically (along one-dimensional sets; V, m, and n) onto the upper P(L) and lower P(L) surfaces, respectively; whereas m and n are regarded as constants. To analyze the behavior of the reduced flow (when C→ 0) on the critical manifold S, we project the reduced system onto the two-dimensional subspace (V and n). The reduced flow can be described by three equations: 1) A differential equation for V, which is given by implicitly differentiating f = 0; 2) the differential equation for n (Eqs 6 and 3) the explicit form of the critical manifold for m (Eq (8)). The implicit differentiation of f = 0 leads to, where dm/dt and dn/dt satisfy Eqs 5 and 6, respectively, and m satisfies Eq 10. Given that the critical manifold S is given as a two-dimensional graph Eq 10, we project the reduced model onto (V, n) as, where f, f, and f are df/dV, df/dm, and df/dn, respectively; and and are dV/dt and dn/dt, respectively. This system is singular along the fold-curves, ∂f / ∂V = 0. Hence, we desingularize the system by rescaling time with τ = t / (∂f / ∂V). The desingularized system is given by, Where V´ and n´ are V and n derived from the new timescale τ. The phase portrait of Eqs 15 and 16 are identical to that of Eqs 13 and 14, respectively; but the direction of flow in the unstable slow manifold (∂f / ∂V > 0) is reversed because of time rescaling. In general, the desingularized system exhibits two forms of singularities—ordinary and folded. An ordinary singularity is an equilibrium of the reduced system (Eqs 13 and 14) and is given by,

By contrast, folded singularities are categorized as node, saddle, or saddle-node based on the type of equilibrium of the desingularized system. Folded singularities are obtained by,

Although folded singularities are not the equilibria of the reduced flow, they make it possible for the reduced flow to cross the fold-curve in finite time. These solutions are called “singular canard”. When arriving at the fold-curve, the absence of folded singularities causes the reduced flow to jump along the fast fibers and move away from the subthreshold region without producing any STOs [62]. A particular type of three-dimensional canards (canard of folded node) is responsible for STOs in MMOs [60]. In the present model, one folded node singularity (with negative eigenvalues) on the lower fold-curve (L) and/or upper curve (L+) exists under control conditions (Fig 4). Under this condition, all parameters are represented by the physiological values introduced by Acker et al. [33].

Fig 4

In control condition, the singularities, V-nullcline (red) and n-nullcline (green) for I = 17.5 (μA/cm2).

In the Region of interest (ROI), the existence of folded node on L- (filled blue cycle, d) indicates the occurrence of type I MMOs (left-trace and right-lower-trace). In this model, there is also a folded node on L+ (filled blue cycle, b), which shows there is also the type II MMOs in this model (right-upper-trace). In control condition, there is also a stable folded focus on L+ (a), folded saddle on L-(c) and saddle point on NN (e).

Results and discussion

The range of parameters

In this section, we determine the parameter ranges in which the desingularized system maintains a folded node singularity. We specifically examine singularities for various values of g, g, and g. Fig 4 shows that the desingularized system exhibits a three-branched n-nullcline (green curves; L, L─, and NN) and a single-branch V-nullcline (red curve) and satisfies fg(V,m) + fh(V,n) = 0.

The curve NN satisfies h(V,n) = 0, and both curves L and L─ satisfy f = 0. Ordinary singularities are placed at the intersection of V-nullcline with NN, whereas folded singularities are placed at the intersections of V-nullcline with L and L─. Fig 5 shows that changing g affects the singularities of the desingularized system by moving V-nullcline up or down. At g < −9.388 (μS/cm2), n-nullcline (L, L−, and NN) does not change and V- and n-nullclines do not intersect; thus, no critical point exists within this range. V-nullcline moves down with inceasing g. At g = −9.388 (μS/cm2), V-nullcline and L+ collide, resulting in the emergence of a saddle-node bifurcation point (SN1), which is called a type-I folded saddle-node and is regarded as a standard saddle-node bifurcation of folded singularities [72, 73]. With further increase in g, two equilibria, namely folded focus (A1) and folded saddle (B1), originate from SN1 and diverge from each other on L+. In addition, a stable node on NN curve (C1) emerges; this node is an ordinary singularity and represents the equilibrium of the Eqs 4−6. As g is increased to 0.514 (μS/cm2), the folded saddle B1 moves to the right and the stable node C1 moves down and to the left. At g = 0.514 (μS/cm2), a transcritical bifurcation (TR1) occurs where B1 and C1 collide and coalesce. This folded singularity bifurcation is a type II folded saddle-node [72, 73].

Fig 5

The folded and ordinary singularities on nullclines for changing g (μS/cm2).

Filled cyan circles (E1 and B2-B8) on L+ and L- curves are stable folded node and unfilled cyan circles (B1 and D1-D5) are folded saddle. Filled blue triangles (SN1 and SN2) are type-I folded saddle-node bifurcation (standard saddle-node bifurcation). Filled blue diamond (HB) is Hopf bifurcation for g = 2.937 (μS/cm2). Filled blue rectangulars (TR1 and TR2) are type-II folded saddle-node bifurcation (transcritical bifurcation). Violet-color filled and unfilled cycles (A1-A8) are stable and unstable folded fuci respectively. Filled and unfilled cyan circles (C1-C6) on NN curve are stable and saddle point respectively, which are ordinary singhularities of system (15) and (16).

As a consequence of this bifurcation, the folded saddle point becomes the folded node. At g between 0.514 and 1.66 (i.e., g = 1 (μS/cm2)), the equilibria on L+ are the folded focus (A3) and the folded node (B2). Moreover, the equilibrium on NN (C3) becomes a saddle point. At g = 1.66 (μS/cm2), the second type I folded saddle-node bifurcation (SN2) occurs on L−. At g between 1.66 and 2.29 (i.e., g = 2 (μS/cm2)), the folded focus (A5) and the folded node (B4) occur on L+, and a saddle point (C3) occurs on NN. Similarly, a folded saddle (D1) and a folded node (E1) occur on L−, which originate from SN2 and subsequently diverge from each other.

The appearance of the folded node on L− corresponds to an MMO, which is a combination of STO and spiking. As g is increased to 2.29 (μS/cm2), the second type II folded saddle-node bifurcation occurs on L−, where E1 and C3 coalesce at a transcritical bifurcation (TR2). At g > 2.29 (μS/cm2), the system does not show type I MMOs. At g between 2.29 and 2.94 (μS/cm2), a stable point (C4) and a folded saddle (E2) occur on NN and L−, respectively. At g = 2.94 (μS/cm2), the system undergoes a supercritical HB on L+, whereas at g > 2.94 (μS/cm2), an unstable folded focus (A8) occurs on L+.

Figs 6 and 7 show the effects of g on singularities of the desingularized system. Varied g values influence V-nullcline and n-nullcline in the (n, V)-phase-plane. At g < 3.35 (μS/cm2), only a stable node (A) exists (Fig 6(A)). At g = 3.35 (μS/cm2), a fold in the critical manifold emerges and SN1 and SN2 occur on the fold-curve (Fig 6(B)).

Fig 6

The folded and ordinary singularities on nullclines for changing g (μS/cm2).

a) g = 1 (μS/cm2), b) g = 3.36 (μS/cm2), c) g = 20 (μS/cm2), and d) g = 40.86 (μS/cm2). The phase space of this system contains a single branch V-nullcline (red line) and three branches n-nullcline (green line) including L+ (upper horizontal line), L- (lower horizontal line) and NN (vertical curve). The color convention for equilibria is similar to Fig 5.

Fig 7

The folded and ordinary singularities on nullclines for changing g (μS/cm2).

a) g = 50 (μS/cm2), b) g = 64.57 (μS/cm2), c) g = 80 (μS/cm2), and d) g = 104.76 (μS/cm2). The phase space of this system contains a single branch V-nullcline (red line) and three branches n-nullcline (green line) including L+ (upper horizontal line), L- (lower horizontal line), and NN (vertical curve). The color convention for equilibria is similar to Fig 6.

With further increase in g, the fold-curve splits into two L+ and L− fold-curves. At g between 3.35 and 40.85 (μS/cm2), a folded saddle and a folded node occur on L+, two folded saddles occur on L−, and a stable point occurs on NN (Fig 6(C)). At g = 40.85 (μS/cm2), the system undergoes a supercritical HB on L+, as shown in Fig 6(D). In addition, the folded saddle on L− coalesces with the stable node on NN through a type II folded saddle-node bifurcation. With further increase in g (i.e., g = 50 (μS/cm2)), a folded focus occurs on L+, a folded node occurs on L−, and a saddle point occurs on NN (Fig 7(A)). Similar to g, the appearance of the folded node on L− corresponds to a type I MMO, which is a combination of STO and spiking.

At g = 64.57 (μS/cm2), the folded saddle and folded node coalesce at a type I folded saddle-node bifurcation (Fig 7(B)). As g is increased further (i.e., g = 80 (μS/cm2)), no folded singularity occurs on L− (Fig 7(C)). At g = 104.46 (μS/cm2), a transcritical bifurcation (type-II folded saddle-node bifurcation) occurs on L+, and two equilibria, namely saddle and unstable points, occur on NN (Fig 7(D)). At g > 104.46 (μS/cm2), the only equilibrium, which is an ordinary stable point, occurs on NN.

On the other hand, changing g does not affect the position of both nullclines, but affects the folded singularities. In this case, the type I folded saddle-node bifurcation at g = 5.92 (μS/cm2) and the type-II folded saddle-node bifurcation at g = 42.6 (μS/cm2) occur on L−. At 5.92 < g < 42.6 (μS/cm2), a folded node on L− and a folded node on L+ occur. This result indicates that both type I and type II MMOs occur in this range. Analysis of this phenomenon is beyond the scope of this paper.

Significance of MMOs

The MMOs are significant because STOs in both subthreshold (type I MMOs) and suprathreshold (type II MMOs) ranges influence the amplification of synaptic inputs, the sensitivity of the neuron to transient inputs, the regulation of neuronal firing rate, and the synchronization of neural network at a specific frequency [11, 12, 38, 74]. In this section, we investigate the effect of I on the STO phenomenon in both types of MMOs from electrophysiological point of view.

In the previous section, the numerical analysis of a desingularized system (Eqs 15 and 16) was presented. According to these results, pyramidal neurons display two types of MMOs in two different ranges of g. Specifically, for 0.514 < g < 1.66 (i.e., g = 1 (μS/cm2)) and for 1.66 < g < 2.29 (i.e., g = 2 (μS/cm2)), there exist one type I MMO and one type II MMO, respectively. These two ranges are within the physiological ranges defined in previous experimental studies [4, 33, 75]. Fig 8 shows the evolution of membrane potential (V), n, and m for two values of these ranges of parameters.

Fig 8

The evolution of variables in the reduced model of pyramidal neurons for I = 17.2 (μA/cm2).

a) For g = 2 (μS/cm2), the membrane potential displays the type I MMO (upper trace), activation variable of K+-channel, n (middle trace), and activation variable of KS-channel, m (lower trace). b) For g = 1 (μS/cm2), the membrane potential shows the type II MMO (upper trace), activation variable of K+-channel, n (middle trace), and activation variable of KS-channel, m (lower trace).

For g = 2 (μS/cm2), the membrane potential starts to depolarize in response to the injection of input current, I (Fig 8(A)). As shown in Fig 2, the time constants of n (τ) and m (τ) are larger than the time constant of membrane potential (τ). Therefore, both the transient potassium (Hodgkin-Huxley potassium channel or K+-channel) and the slow non-activating potassium channel (KS-channel) are activated with a considerable delay. Particularly, KS-channel is activated with more delay than K+-channel because of its very large value of time constant (τ = 90 (msec)). Therefore, for the first few seconds of membrane oscillation, these channels cannot oppose the changes in membrane potential. As a result, the membrane potential (V) sharply rises to the depolarization values. At V ~ −55 (mV), however, the K+-channel and, with delay, the KS-channel are activated (n and m start to increase). The activation of these channels moves the positive K+ ions from the intracellular to the extracellular region, which in turn hyperpolarizes the membrane potential.

This displacement forms a force against the membrane depolarization which is similar to the behavior of resonator ion channels discussed in pyloric dilator neurons [35]. At V ~ −45 (mV), K+- and KS- channels can neutralize the effect of I and decrease the membrane potential; this, in turn, partially inactivates both ion channels. At this point, the partial inactivation of K+- and KS- channels allows I to increase the membrane potential to depolarization values, which in turn activates both ion channels, but with considerable delays. The interaction between I and the K+- and KS- channels leads to the small-amplitude oscillation (STOs) in the membrane potential.

In the last convex curve of STOs (before the action potential), K+- and KS- channels hyperpolarize the membrane potential more than before (downwards convex curve), which in turn, decreases the values of n and m more than the previous period (further inactivates K+- and KS- channels). Consequently, I depolarizes the membrane potential more sharply than the previous periods and triggers the action potential. Note that the time constant values of K+- and KS- channels are larger than the time constant of V. Therefore, they cannot be activated fast enough to avoid the depolarization of membrane potential towards the threshold of action potential.

Fig 8 (B) shows that the membrane potential of pyramidal neurons can display type II MMOs for g = 1 (μS/cm2). Decreasing g decreases the effect of the terms g. m. (V − E) in Eq 4. Therefore, the potassium channels (K+- and KS- channels) cannot counterbalance the effect of I on rapid depolarization of the membrane potential. Therefore, the injection of I sharply increases the membrane potential and triggers the action potential. At the top of the action potential, the K+- and KS- channels are properly activated to sharply decrease the membrane potential to −23 (mV).

It should be noted that the decrease in membrane potential is mostly caused by the K+-channel (variable n). At this point, the KS-channel does not reach its maximum activation value. The interaction between I and the potassium channels (mostly the K+-channel) leads to the small-amplitude oscillation (STOs) in the membrane potential. After several STOs, the activation value of the KS-channel reaches its maximum value. At this point, the rate of K+ ions moving from the intracellular to the extracellular region is very high, resulting in a sharp decrease of the membrane potential.

Although canard theory has been previously used to analyze the behavior of neuronal models under normal conditions [76], our study is the first to apply canard theory to examine a neuronal model under antiepileptic drug conditions. Additionally, our mathematical analysis shows that MMOs of type 1 and type 2 could be seen in pyramidal neurons when I is inhibited but in different ranges of g. To our knowledge, this is the first time the two types of MMOs have been demonstrated in a single neuronal model. Importantly, there is mounting evidence in literature that supports such predictions. For instance, both types of MMOs have been observed in pyramidal neurons under antiepileptic drug condition when I was blocked using pharmacological agents [13, 74, 77]. This shows that our results are biologically valid. These results indicate that the type of MMO depends on the maximum conductance (g) of I. Accordingly, this explains why MMOs of type 1 or type 2 have been seen in these experimental studies–probably due to the magnitude of I in these studies.

Conclusion

Although several studies have applied canard theory to different neural models and predicted important mathematical characteristics of these models, the discussion of the electrophysiological mechanisms that underlie the mathematical predictions has attracted less attention [64, 73, 78, 79]. This study provides a more in-depth investigation of the relationship between the mathematical characteristics and the electrophysiological mechanisms of MMOs in pyramidal neurons under antiepileptic drug conditions. It also showed that the magnitude of slow potassium currents is critical for determining the type of the emerging MMO. For the questions posed in this study, slow outward currents, as opposed to fast inward currents, have a major impact on the properties of the MMO phenomena.

Supplementary information: The ion channel gating kinetics of pyramidal neurons.

(DOCX)

Click here for additional data file.

MATLAB code of Fig 2.

(M)

Click here for additional data file.

MATLAB code of Fig 3.

(M)

Click here for additional data file.

MATLAB code of Fig 4 and Fig 5.

(M)

Click here for additional data file.

MATLAB code of Fig 6 and Fig 7.

(M)

Click here for additional data file.

MATLAB code of Fig 8.

(M)

Click here for additional data file.