Recently, in the past decade, high-frequency oscillations (HFOs), very high-frequency oscillations (VHFOs), and ultra-fast oscillations (UFOs) were reported in epileptic patients with drug-resistant epilepsy. However, to this day, the physiological origin of these events has yet to be understood. Our study establishes a mathematical framework based on bifurcation theory for investigating the occurrence of VHFOs and UFOs in depth EEG signals of patients with focal epilepsy, focusing on the potential role of reduced connection strength between neurons in an epileptic focus. We demonstrate that synchronization of a weakly coupled network can generate very and ultra high-frequency signals detectable by nearby microelectrodes. In particular, we show that a bistability region enables the persistence of phase-shift synchronized clusters of neurons. This phenomenon is observed for different hippocampal neuron models, including Morris–Lecar, Destexhe–Paré, and an interneuron model. The mechanism seems to be robust for small coupling, and it also persists with random noise affecting the external current. Our findings suggest that weakened neuronal connections could contribute to the production of oscillations with frequencies above 1000 Hz, which could advance our understanding of epilepsy pathology and potentially improve treatment strategies. However, further exploration of various coupling types and complex network models is needed.

We have built a mathematical framework to examine how a reduced neuronal coupling within an epileptic focus could lead to very high-frequency (VHFOs) and ultra-fast oscillations (UFOs) in depth EEG signals. By analyzing weakly coupled neurons, we found a bistability synchronization region where in-phase and anti-phase synchrony persist. These dynamics can be detected as very high-frequency EEG signals. The principle of weak coupling aligns with the disturbances in neuronal connections often observed in epilepsy; moreover, VHFOs are important markers of epileptogenicity. Our findings point to the potential significance of weakened neuronal connections in producing VHFOs and UFOs related to focal epilepsy. This could enhance our understanding of brain disorders. We emphasize the need for further investigations of weakly coupled neurons.

It is known that hippocampal pyramidal neurons can switch from integrators to resonators, that is, from class I excitability to class II (Prescott, Ratté, De Koninck, & Sejnowski, 2008). This dual operating mode allows for individual spikes, tonic spiking, and bursting, as well as different firing rates (Ermentrout & Terman, 2010; Izhikevich, 2006). Specifically, there is a maximum firing frequency above which the potential on the neuron membrane cannot oscillate faster; oscillations occur at a certain natural frequency for a given external stimulus. Although different types of neurons have different rates of action potential firing, these are physiologically limited (Gabbiani & Cox, 2017; Purves et al., 2019). The fI curve (frequency given a fixed input I) typically saturates to a finite level not much above 350 Hz, even for the fastest spiking neurons; see Wang et al. (2016). Using mathematical conductance-based models, we can explain the onset and offset of the frequency–input curve. For the offset, one typically encounters a Hopf bifurcation or a limit point of cycles (Izhikevich, 2006; Kuznetsov, 2023). Biophysically plausible parameter settings then show we cannot expect higher frequencies. The dependence of the individual neuron model dynamics on the external current and other parameters has been studied extensively (Cessac & Samuelides, 2007; Duan & Lu, 2006; Shilnikov, 2012; Tsumoto, Kitajima, Yoshinaga, Aihara, & Kawakami, 2006; Xing, Song, Wang, Yang, & Chen, 2022). This biophysical basis of the mechanism is also confirmed by comparison with a real neuronal signal (Prescott, De Koninck, & Sejnowski, 2008).

Nevertheless, hippocampal electroencephalographic signals (EEG) show that high-frequency oscillations (HFOs) are reported in the seizure onset zone and are connected with epileptic seizure generation in humans with focal epilepsy (Cimbalnik et al., 2018, 2020; Jacobs et al., 2008, 2012; Jiruska et al., 2017; Pail et al., 2020; Řehulka et al., 2019; Staba & Bragin, 2011; Worrell & Gotman, 2011). Moreover, the extent of HFOs correlates with seizure frequency and disease severity (Zijlmans, Jacobs, Zelmann, Dubeau, & Gotman, 2009). HFOs are most commonly distinguished by frequency bands as interictal high gamma (65–100 Hz), ripples (100–250 Hz), and fast ripples (250–600 Hz). Recently, very high-frequency oscillations (VHFOs, 600–2000 Hz; Brázdil et al., 2017; Travnicek et al., 2021; Vasickova et al., 2023), very fast ripples (VFRs, 600–1000 Hz), and ultra-fast ripples (UFRs, 1000–2000 Hz), and most recently ultra-fast oscillations (UFOs, over 2000 Hz) were reported in depth EEG recordings of patients with epilepsy (Brazdil et al., 2023). Since VHFOs were more spatially restricted in the brain than HFOs with lower frequencies, they have been suggested as novel biomarkers of the epileptogenic zone. Moreover, VFRs present intriguing synchronization phenomena, demonstrating strong phase-locking with slow oscillations (Hao et al., 2021). This highlights the crucial role synchronization mechanisms may play in the pathological progression and underlying propagation of epileptiform discharges within the hippocampal network during status epilepticus in temporal lobe epilepsy. However, the neurons themselves are unable to produce such high oscillation frequencies. This brings us to the question of whether multiple coupled neurons are capable of producing a synchronized signal of such high frequency.

It is known that synchronization occurs in oscillators for sufficiently strong coupling (Boccaletti, Pisarchik, Del Genio, & Amann, 2018; Pikovsky, Rosenblum, & Kurths, 2001; Strogatz, 2018; Wiggins, 2003) and even that it can occur through external noise (Lang, Lu, & Kurths, 2010; Nagai & Kori, 2010; Zhou & Kurths, 2002). At the same time, the dynamics of coupled oscillators may be complex, chaotic, or consist of chimera states (Abrams & Strogatz, 2004; Aihara, Takabe, & Toyoda, 1990; Calim, Hövel, Ozer, & Uzuntarla, 2018; Hoppensteadt & Izhikevich, 1997; Majhi, Bera, Ghosh, & Perc, 2019). Therefore, it is important to determine how coupled neuron models behave dynamically for physiological parameters and whether and, if so, how they can be synchronized with higher than natural frequencies, especially within VHFO and UFO frequency bands.

Within this paper, we consider near-identical neurons weakly connected through gap junctions. In this setup, we show that local field potentials can exhibit VHFOs and UFOs through phase synchronization. Moreover, we claim that it is not an exceptional phenomenon but quite the opposite in the case of very weak neuronal couplings. It arises for various connectivity types and neuron models. The mechanism is based on theoretical works on the symmetry of dynamical networks that can also be spatiotemporal and can lead to a rigidphase-shift synchrony (Golubitsky, Messi, & Spardy, 2016; Golubitsky, Romano, & Wang, 2012; Golubitsky & Stewart, 2006, 2016; Golubitsky, Stewart, & Török, 2005). We show that the occurrence of these higher frequencies is possible both in the case of different values of model parameters and in the case of random noise affecting the magnitude of the external current input.

In the Results section, we first discuss the Morris–Lecar neuron model, which we examine thoroughly. This includes a discussion of its fundamental dynamics, followed by an exploration of Morris–Lecar neuronal networks. We then analyze the behavior of two coupled Morris–Lecar neuron models, explain the phase shift, and reveal the concept of anti-phase collective synchrony in a model of a neuronal network. We show that we can achieve transient anti-phase synchrony by stimulating a subnetwork of neurons due to the existence of the bistability region. This concept can lead to an apparently higher frequency in the EEG signal. We want to demonstrate that the observed phenomena of phase-shift synchrony and bistability are generic for weak coupling and still hold for more physiological models of neurons. The analysis of the Morris–Lecar model is straightforward as the single-neuron model is two-dimensional. The hippocampal neuronal network, however, is considerably more complex than the simple networks composed of Morris–Lecar neuron models, which exhibit a specific course of spike dynamics. We show that the bistability and anti-phase synchronization appear robustly in network models with more complex neuronal dynamics. One model is an interneuron network that exhibits VHFOs and UFOs with anti-phase locking, and the other is a Destexhe–Paré neuron network that shows VHFOs transiently. Similarly to the Morris–Lecar model, investigating the dynamics of two coupled interneurons, we find regions with synchronization and bistability and present simulations of the collective anti-phase behavior. We discuss transient ultra-fast ripples and ultra-fast oscillations in a two-level neuronal network. Lastly, simulations of the collective anti-phase behavior are presented for the Destexhe–Paré neuron model. In the Discussion and Conclusions section, we detail our primary objective to investigate whether weakened connections between neurons in an epileptic focus could result in high, very high-frequency, and ultra-fast oscillations; we examine how our findings support this theory and how they provide a mathematical framework for interpreting EEG signals from patients with focal epilepsy, while also highlighting the need for future research on the effects of different coupling types, strengths, and parameters on the dynamics of coupled neurons.

Model Dynamics Explaining the Tonic Spiking

Let us recall basic neuron spiking mechanisms that lead to tonic spiking clearly described in Prescott, De Koninck, and Sejnowski (2008). Class I neurons continuously increase their firing rate in response to a steady increase in input current without any abrupt change in their behavior. This type of behavior arises from a saddle-node on an invariant curve (SNIC) bifurcation, characterized by a break of a limit cycle on a saddle-node equilibrium as the input current is varied. On the contrary, class II neurons exhibit a discontinuous increase in their firing rate in response to a change in the input current, with a threshold for tonic spiking emergence. This behavior occurs near a subcritical Hopf bifurcation when a stable spiking orbit already exists due to a saddle-node of limit cycles with respect to the input current. The bifurcation diagram presented in Figure 1 demonstrates that the Morris–Lecar system Equation 11 shows both bifurcations at onset and offset as we vary the external input external current Iext, as mentioned. Note that class III neurons do not possess tonic spiking dynamics with a given natural frequency; hence, we will not focus on their type of mechanism. From now on, we will assume the neuron exhibits tonic spiking dynamics with a certain given natural frequency under a suitable external stimulus current.

Figure 1.

Morris–Lecar model dynamics with respect to the applied external current Iext. Left: The bifurcation diagram shows a limit point (LP) for Iext ≈ 40 where a saddle-node on an invariant curve (SNIC) leads to a stable cycle, which exists for a wide range of external current values (tonic spiking region) until the limit point of cycles near Iext = 140. The branch of cycles turns and becomes unstable, and ends at a subcritical Hopf bifurcation (H). Beyond the LPC, neuronal excitation is impossible. The other parameter values are listed in Table 1, and the membrane capacitance is C = 5 for simplicity. Right: The natural firing frequency as a function of external current Iext for various membrane capacitances C. As the capacitance increases, only the frequency of the oscillation decreases but does not affect the existence of the limit cycle. Other parameter values are listed in Table 1.

Figure 1.

Morris–Lecar model dynamics with respect to the applied external current Iext. Left: The bifurcation diagram shows a limit point (LP) for Iext ≈ 40 where a saddle-node on an invariant curve (SNIC) leads to a stable cycle, which exists for a wide range of external current values (tonic spiking region) until the limit point of cycles near Iext = 140. The branch of cycles turns and becomes unstable, and ends at a subcritical Hopf bifurcation (H). Beyond the LPC, neuronal excitation is impossible. The other parameter values are listed in Table 1, and the membrane capacitance is C = 5 for simplicity. Right: The natural firing frequency as a function of external current Iext for various membrane capacitances C. As the capacitance increases, only the frequency of the oscillation decreases but does not affect the existence of the limit cycle. Other parameter values are listed in Table 1.

Close modal

External current and parameter changes influence the neuron natural frequency; see Figure 1. Membrane capacitance depends on the cell’s and membrane’s biophysical properties, including fluidity, thickness, permeability, and so forth. (Bakhtiari, Manshadi, Candas, & Beskok, 2023). Moreover, recently published studies (Patel, Tewari, Chaunsali, & Sontheimer, 2019; Tewari et al., 2018) propose pathological changes in the effective capacitance of neurons that may contribute to epilepsy. Therefore, the capacitance appears to be an appropriate parameter for the bifurcation analysis in relation to the emergence of HFOs, VHFOs, and UFOs in epilepsy patients. Similarly, one can investigate the dependence on other model parameters influencing the neuron frequency, for example, φ and Iext (e.g., see Liu, Liu, & Liu, 2014; Tsumoto et al., 2006; Xing et al., 2022) or gCa, gK, and gL that can also be influenced by drugs (see Macdonald, 1989; Nicholson, Blanche, Mansfield, & Tran, 2002).

Two Coupled Morris–Lecar Neuron Models

To get insight into the behavior of a large Morris–Lecar neuronal network, let us focus on the smallest possible network first. Consider a pair of bidirectionally coupled Morris–Lecar neurons (Equation 15) for N = 2. It is well known that if the ratio of the frequencies of two oscillators is irrational, the orbit lives on a two-dimensional invariant torus and is called asynchronous or quasi-periodic (Wiggins, 2003). For a strong enough coupling ε between them, the dynamics of both oscillators can become mutually synchronized. If the phases of two synchronized oscillators are similar or opposite, we refer to this as in-phase (IPS, Figure 2B) and anti-phase synchrony (APS, Figure 2C), respectively; see Pikovsky et al. (2001). For very different frequencies, even strong coupling might not be capable of synchronizing the dynamics of the coupled systems, see Figure 2A, while for close frequencies, the dynamics are quasi-periodic only for very weak coupling.

Figure 2.

Dynamics in the system (Equation 15) for ε = 0.05, various values of C1, and fixed C2 = 1 (see Figure 4). Left: projections of orbits (blue dotted line), and stable (blue solid line) and unstable (red dashed line) cycles onto the state space V1 × V2 found by continuation. Right: time course of V1 (blue) and V2 (red) of the blue orbits in the left panel. Other parameter values are listed in Table 1. Choosing C1 affects the dynamics: (A) C1 = 2.75, quasi-periodicity. (B) C1 = 1.60, a stable in-phase periodic solution (IPS; the bold blue cycle on the left state space V1 × V2). (C) C1 = 0.90, bistability, the system possesses stable in-phase and anti-phase periodic solutions (IPS and APS; the bold blue cycles on the left state space V1 × V2). (D) C1 = 0.581, bistability, anti-phase cycle lost stability through a period-doubling bifurcation (PD).

Figure 2.

Dynamics in the system (Equation 15) for ε = 0.05, various values of C1, and fixed C2 = 1 (see Figure 4). Left: projections of orbits (blue dotted line), and stable (blue solid line) and unstable (red dashed line) cycles onto the state space V1 × V2 found by continuation. Right: time course of V1 (blue) and V2 (red) of the blue orbits in the left panel. Other parameter values are listed in Table 1. Choosing C1 affects the dynamics: (A) C1 = 2.75, quasi-periodicity. (B) C1 = 1.60, a stable in-phase periodic solution (IPS; the bold blue cycle on the left state space V1 × V2). (C) C1 = 0.90, bistability, the system possesses stable in-phase and anti-phase periodic solutions (IPS and APS; the bold blue cycles on the left state space V1 × V2). (D) C1 = 0.581, bistability, anti-phase cycle lost stability through a period-doubling bifurcation (PD).

Close modal

Phase Shift for Two Weakly Coupled Morris–Lecar Neuron Models

Two identical oscillators form a symmetrical dynamical system with a symmetrical periodic solution. The symmetry of the system is roughly preserved even for non-identical oscillators due to the continuous dependence of the solutions on parameters (see Perko, 2001), that is, coupled neurons with similar phase dynamics form a nearly symmetrical system with a stable limit cycle such that the phases of both tonically spiking neurons are mutually shifted, although having a common frequency roughly the magnitude of the natural frequencies of the individual neurons. This phenomenon, so-called phase-locking, and the theoretical background for this mechanism is based on Golubitsky and Stewart’s results of equivariant bifurcation theory allowing groupoid formalism for networks of systems of differential equations (Golubitsky et al., 2012, 2016; Golubitsky & Stewart, 2006, 2016; Golubitsky et al., 2005; Nijholt, Rink, & Sanders, 2019). The symmetry of the network with a T-periodic solution implies the existence and persistence of a rigid phase-shift synchrony. In our case of two nearly identical weakly coupled oscillators (Equation 15), there may also exist solutions synchronized with time shift T/2 for appropriate system parameters: the attracting stable cycle, corresponding to the presence of an anti-phase synchronous state (see Figure 2C). Simultaneously, the system possesses a solution reflecting in-phase synchrony where the neurons oscillate almost identically; see Figures 2B, 2C (middle), and 2D (middle).

Both branches of stable periodic solutions, that is, in-phase and anti-phase synchrony, persist inside relatively wide regions of the parameter space; see Figure 3. This fact is also depicted in Figure 4 where the bifurcation diagrams with respect to C1 and C2, and ε and C1, respectively, are shown. More accurately, Figure 3 precisely corresponds to the dashed sections at C2 = 1 and ε = 0.05 in Figure 4. The regions, so-called Arnold tongues or resonance tongues (Kuznetsov, 2023; Wiggins, 2003), are delineated by borders belonging to the limit point of cycle (LPC) manifolds. Outside the Arnold tongues, the difference in capacitances C1 and C2 implies a synchronous state with heterogeneous frequencies (see Figure 1), quasi-periodic oscillations, or more complex dynamics. The LPC manifolds were computed for parameters listed in Table 1. Since interchanging C1 and C2 mirrors the synchronization regions, the black and blue solid curves in Figure 4A are symmetrical. Cusps in the intersections of two branches of LPC manifolds are co-dimension two bifurcation points (Kuznetsov, 2023). Although the LPC curves in Figure 4 were computed for given parameters from Table 1, the generality of bifurcation theory guarantees that these LPC manifolds persist for nearby parameter values.

Figure 3.

Frequency of the limit cycles in the system (Equation 15) of two coupled Morris–Lecar neuron models with respect to C1. The remaining parameter values are identical (see Table 1), C2 = 1, and ε = 0.05. The upper and bottom blue branches correspond to the stable in-phase and anti-phase synchrony, respectively. The dashed sections at the given parameter values in Figure 4 correspond to these depicted branches. Notice the bistability, present for a relatively wide range of C1. For completeness, we remark that there exist more and more complex periodic solutions of the system (Equation 15) for C1 between the LPC3 and PD1 points, and between PD2 and PD3 points that are not depicted in this figure. The dashed line at C1 = 1.2 denotes the frequencies at IPS near 30 Hz and APS near 26 Hz that reappear in the simulations depicted in Figures S1 and S2 in the Supporting Information.

Figure 3.

Frequency of the limit cycles in the system (Equation 15) of two coupled Morris–Lecar neuron models with respect to C1. The remaining parameter values are identical (see Table 1), C2 = 1, and ε = 0.05. The upper and bottom blue branches correspond to the stable in-phase and anti-phase synchrony, respectively. The dashed sections at the given parameter values in Figure 4 correspond to these depicted branches. Notice the bistability, present for a relatively wide range of C1. For completeness, we remark that there exist more and more complex periodic solutions of the system (Equation 15) for C1 between the LPC3 and PD1 points, and between PD2 and PD3 points that are not depicted in this figure. The dashed line at C1 = 1.2 denotes the frequencies at IPS near 30 Hz and APS near 26 Hz that reappear in the simulations depicted in Figures S1 and S2 in the Supporting Information.

Close modal
Figure 4.

Arnold tongues in parameter spaces (A) C1 × C2 and (B) C1 × ε for two coupled Morris–Lecar neuron models (Equation 15). The regions are delineated by blue and black solid LPC curves corresponding to the IPS and APS, respectively. The blue dash-dotted curves refer to PD of the anti-phase cycle. The background color indicates the dominant frequency in a simulated composed signal V = V1 + V2. In the region bounded by the PD curves, the system may exhibit more complex, even chaotic dynamics. The red region corresponds to the stable anti-phase cycle, while the other solutions (yellow) have a frequency close to the natural frequency. Compare with Figure 2: (A) quasi-periodicity outside Arnold tongues; (B) IPS inside the outer Arnold tongue; (C) bistability in the red region delineated by LPC curves and PD curves; (D) APS lost stability due to PD, frequency halves.

Figure 4.

Arnold tongues in parameter spaces (A) C1 × C2 and (B) C1 × ε for two coupled Morris–Lecar neuron models (Equation 15). The regions are delineated by blue and black solid LPC curves corresponding to the IPS and APS, respectively. The blue dash-dotted curves refer to PD of the anti-phase cycle. The background color indicates the dominant frequency in a simulated composed signal V = V1 + V2. In the region bounded by the PD curves, the system may exhibit more complex, even chaotic dynamics. The red region corresponds to the stable anti-phase cycle, while the other solutions (yellow) have a frequency close to the natural frequency. Compare with Figure 2: (A) quasi-periodicity outside Arnold tongues; (B) IPS inside the outer Arnold tongue; (C) bistability in the red region delineated by LPC curves and PD curves; (D) APS lost stability due to PD, frequency halves.

Close modal
Table 1.

Typical physiological parameter values for the Morris–Lecar neuron model of a hippocampal pyramidal cell (Gutkin & Ermentrout, 1998)

 SettingUnitMeaning
C 0.1–2.6 μF/cm2 membrane capacitance 
gL mS/cm2 maximum leak conductance 
gCa mS/cm2 maximum Ca2+ conductance 
gK mS/cm2 maximum K+ conductance 
VL −60 mV equilibrium potential of leak channel 
VCa 120 mV equilibrium potential of Ca2+ channel 
VK −80 mV equilibrium potential of K+ channel 
β1 −1.2 mV tuning parameters for Ca2+ activation function 
β2 18 mV 
β3 10 mV tuning parameters for K+ activation function and time function τw 
β4 17.4 mV 
φ 1/15 s−1 reference frequency 
Iext 43 μA/cm2 externally applied direct current 
 SettingUnitMeaning
C 0.1–2.6 μF/cm2 membrane capacitance 
gL mS/cm2 maximum leak conductance 
gCa mS/cm2 maximum Ca2+ conductance 
gK mS/cm2 maximum K+ conductance 
VL −60 mV equilibrium potential of leak channel 
VCa 120 mV equilibrium potential of Ca2+ channel 
VK −80 mV equilibrium potential of K+ channel 
β1 −1.2 mV tuning parameters for Ca2+ activation function 
β2 18 mV 
β3 10 mV tuning parameters for K+ activation function and time function τw 
β4 17.4 mV 
φ 1/15 s−1 reference frequency 
Iext 43 μA/cm2 externally applied direct current 

Similarly, the bifurcation diagram in Figure 4B shows the Arnold tongues with respect to the coupling strength ε and membrane capacitance C1. In this case, they emanate from the point C1 = C2 (other parameters are identical) since the limit case ε = 0 belongs to the uncoupled pair of neurons with identical frequencies. Generally, nonidentical similar neurons will have similar intrinsic frequencies, and the Arnold tongue would emanate from a cusp point with the same frequency, but the membrane capacitances C1, C2, as well as other parameters, would differ slightly. Moreover, results presented in Kobelevskiy (2008), which deals with Morris–Lecar models that incorporate time-delayed gap-junctional coupling, suggest that similar dynamics will emerge in such models. Furthermore, an increase in delay may even play a stabilizing role.

The background color in the bifurcation diagrams in Figure 4A and 4B indicates the dominant frequency detected in a composed signal V = V1 + V2 simulated on a time interval [200, 1,200], discarding 200 ms transients, for parameter values from a 401 × 401 grid and appropriate initial conditions. These regions may not match precisely since the two-cycle created by crossing the PD curves has a similar spike shape to the cycle corresponding to the APS (see Figure 2D). Furthermore, around the point [1, 0] in Figure 4B, that is, for very weak coupling and similar membrane capacitances, the observed dominant frequency of ∼60 Hz in the composition is caused by long-lasting transient dynamics.

Anti-Phase Synchrony in a Model of Morris–Lecar Neuronal Network Leading to Higher Collective Frequency

A similar approach can be used for more complex networks of (nearly) identical neurons with symmetry. The phenomenon leads to complete collective synchrony or phase-shift synchrony in the network of neurons.

Experiments, simulations, and theory show that anti-phase collective synchronization is possible (Chowdhury, Rakshit, Buldu, Ghosh, & Hens, 2021; Kawamura, Nakao, Arai, Kori, & Kuramoto, 2010; Sebek, Kawamura, Nott, & Kiss, 2019). Ermentrout and Terman (2010) previously reported on the anti-phase dynamics of two Morris–Lecar neurons with synaptic coupling; however, our findings demonstrate that this phenomenon is not exceptional and has the potential to lead to a rise of higher frequency in the summed signal. Furthermore, we have observed similar collective dynamics in a network, highlighting the robustness of this phenomenon across different scales of neuronal systems.

Consider a neuronal network model (Equation 15) of N = 50 weakly coupled Morris–Lecar neurons. Based on the results presented in the previous subsection, one can expect the presence of multistability in this system, that is, the coexistence of multiple stable solutions corresponding to various synchronous states depending on the initial conditions and varied model parameters ε and Ci, as is depicted in Figures S1 and S2 in the Supporting Information.

Specifically, for the demonstration of the collective anti-phase behavior, let us suppose identical coupling strength ε = 0.1/50 = 0.002. Further, we will consider that the heterogeneity of neurons results from differences in their membrane capacitances (one may assume similar heterogeneity in other parameters, as mentioned before). Let Ci be independent and follow a truncated normal (TN) distribution defined on the interval [0.75, 1.05] to avoid unrealistic or extreme values, that is, Ci ∼ TN(0.90, 0.052, 0.75, 1.05), while the remaining parameters stay identical (see Table 1). We base our choice for the mean and standard deviation of neuron membrane capacitances Ci on direct measurements from Gentet, Stuart, and Clements (2000). To account for possible pathological changes near the epileptic focus, we also explored high values of the standard deviation. Finally, let the population consist of two subpopulations of roughly the same size, 𝒫1 = {1, …, 25}, 𝒫2 = {26, …, 50}, that have random initial conditions that partially overlap
Vi0TN35525020,i𝒫1𝒫2,wi0TN0.010.025200.085,i𝒫1,TN0.060.025200.135,i𝒫2.
(1)
Also, we set Iext ∼ N(43, 12).

Figure 5 shows that this setting leads to anti-phase synchronous oscillations of the two clusters, resulting in a double frequency in the composed signal, as one can see in the periodograms of the clusters and the composition, respectively. Namely, the observed frequency increased to ∼56 Hz. Notice that individual neurons can jump from one cluster to the other due to the noise, while the global behavior remains unchanged.

Figure 5.

(A) Collective APS of two subpopulations in a neuronal network model (Equation 15) composed of 50 all-to-all coupled Morris–Lecar neurons with heterogeneous membrane capacitances Ci ∼ TN(0.90, 0.052, 0.75, 1.05) and random initial conditions in Equation 1. This setting leads to a doubled dominant frequency in (B) the composed signal V = i=150Vi. Zooming in on the dynamics of (C) all neurons, (D) composed signal, and the corresponding periodograms of (E) individual neurons and (F) the composed signal. The robustness of the anti-phase behavior is demonstrated by a noisy applied current Iext ∼ N(43, 12) at each simulation step Δt = 0.01. The remaining parameter values are identical (see Table 1), ε = 0.002.

Figure 5.

(A) Collective APS of two subpopulations in a neuronal network model (Equation 15) composed of 50 all-to-all coupled Morris–Lecar neurons with heterogeneous membrane capacitances Ci ∼ TN(0.90, 0.052, 0.75, 1.05) and random initial conditions in Equation 1. This setting leads to a doubled dominant frequency in (B) the composed signal V = i=150Vi. Zooming in on the dynamics of (C) all neurons, (D) composed signal, and the corresponding periodograms of (E) individual neurons and (F) the composed signal. The robustness of the anti-phase behavior is demonstrated by a noisy applied current Iext ∼ N(43, 12) at each simulation step Δt = 0.01. The remaining parameter values are identical (see Table 1), ε = 0.002.

Close modal

We acknowledge that the very high-frequency events are only transient phenomena in the experimental observations. We have just shown that with small levels of noise, the anti-phase solution persists. However, higher levels typically lead to complete synchronization (Neiman, Schimansky-Geier, Cornell-Bell, & Moss, 1999; Pikovsky et al., 2001; Zhou & Kurths, 2003). Nevertheless, Figure 6 demonstrates that higher frequencies may temporally appear through external inputs. Let us change the applied current to Iext = 43 + pi(t) + σ dW with σ = 2.0 and pi = 40 if t ∈ [716, 718] ms in Equation 15 for i ∈ 𝒫1. This stimulus is applied only to half of the population, that is, pi = 0 for i ∈ 𝒫2. We initialize the population into the anti-phase behavior and subsequently observe that the population synchronizes within 100 ms. Looking at the same time interval as before, the external stimulus arrives to halve the population evoking anti-phase behavior that lasts for ∼200 ms; then the population synchronizes again. As the periodograms and time series show, it is evident that the higher frequencies are temporally present. This transient phase solution also appears if we stimulate fewer neurons but then lasts only a few cycles (data not shown). This transient phenomenon depends on the timing of the stimulus as a stimulus during different phases of the period does not induce the anti-phase solution.

Figure 6.

(A) Transient collective APS in a neuronal network model (Equation 15) composed of 50 all-to-all coupled Morris–Lecar neurons with heterogeneous membrane capacitances Ci ∼ TN(1.04, 0.032, 0.95, 1.13) and random initial conditions in Equation 1. The parameter values remain unchanged from the previous setting, with ε = 0.002. The APS is evoked by 2 ms lasting additional stimulus at t = 716 (arrow). Zooming in on the dynamics of (B) all neurons, (D) composed signal, and the corresponding periodograms of (C) individual neurons and (E) the composed signal.

Figure 6.

(A) Transient collective APS in a neuronal network model (Equation 15) composed of 50 all-to-all coupled Morris–Lecar neurons with heterogeneous membrane capacitances Ci ∼ TN(1.04, 0.032, 0.95, 1.13) and random initial conditions in Equation 1. The parameter values remain unchanged from the previous setting, with ε = 0.002. The APS is evoked by 2 ms lasting additional stimulus at t = 716 (arrow). Zooming in on the dynamics of (B) all neurons, (D) composed signal, and the corresponding periodograms of (C) individual neurons and (E) the composed signal.

Close modal

Model Dynamics for Two Coupled Interneurons

The dynamics of two interacting interneurons coupled through gap-junctional connections has already been studied in Chow and Kopell (2000) and Skinner, Zhang, Velazquez, and Carlen (1999). At low coupling strengths and very high firing rates, the synchronous state is unstable, and a pair of cells fires in anti-phase synchrony, while for a lower range of frequencies, the in-phase and anti-phase synchronies may be bistable (Chow & Kopell, 2000). Moreover, Skinner et al. (1999) found that blocking the inhibition in a network of interneurons coupled by gap junctions and inhibitory synapses can lead to anti-phase bursting with higher frequency.

Similarly, as in the case of Morris–Lecar neuron, let us focus briefly on the dynamics of the interneuron model with respect to parameters C1, C2, and ε. One can proceed analogously with respect to other parameters affecting the frequency of the membrane potential. Figure 7 shows Arnold tongues in parameter spaces C1 × C2 and C1 × ε (fI curves for a single interneuron for various values of C are provided in the Supporting Information (see Figure S4); the default parameter choice of Iext = 24 leads to a frequency of 335 Hz). The blue LPC curves delineate a pair of Arnold tongues corresponding to two symmetrical in-phase solutions; the black tongue refers to the anti-phase synchrony. For the sake of completeness, the red dashed line denotes the LPC curves of two unstable periodic solutions related to the APS. The background color illustrates the dominant frequency detected in a composed signal V = V1 + V2 simulated on a time interval [500, 1,500], discarding 500 ms transients, for parameter values from a 401 × 401 grid and appropriate initial conditions. These regions do not match precisely, for example, around the point [1, 0] in Figure 7B, that is, for very weak coupling and similar membrane capacitances, where the observed dominant frequency of ∼670 Hz in the composition is caused by transient dynamics (similar to the behavior depicted in Figure 8).

Figure 7.

Arnold tongues in parameter spaces (A) C1 × C2 and (B) C1 × ε for two coupled interneuron models (Equation 17 and 19) with Iext = 24. The regions are delineated by blue and black LPC curves corresponding to the IPS and APS , respectively. The red dashed curves refer to the LPC bifurcation of two unstable periodic solutions. The background color indicates the dominant frequency in a simulated composed signal V = V1 + V2, reaching up to ∼750 Hz. The yellow region inside the APS tongue delimited by black curves belongs to the signals where the multiple frequency persists, but is not dominant.

Figure 7.

Arnold tongues in parameter spaces (A) C1 × C2 and (B) C1 × ε for two coupled interneuron models (Equation 17 and 19) with Iext = 24. The regions are delineated by blue and black LPC curves corresponding to the IPS and APS , respectively. The red dashed curves refer to the LPC bifurcation of two unstable periodic solutions. The background color indicates the dominant frequency in a simulated composed signal V = V1 + V2, reaching up to ∼750 Hz. The yellow region inside the APS tongue delimited by black curves belongs to the signals where the multiple frequency persists, but is not dominant.

Close modal
Figure 8.

(A) Transient anti-phase dynamics of two subpopulations in a neuronal network model (Equation 17 and 19) composed of 50 all-to-all coupled interneurons with heterogeneous membrane capacitances Ci ∼ TN(1, 0.032, 0.91, 1.09), coupling (Equation 2), and initial conditions (Equation 3). This configuration leads to the rise of VHFOs in (B) the composed signal V = i=150Vi. Dynamics of (C) all neurons are zoomed, (D) composed, and depicted with periodograms of (E) individual neurons and (F) the composed signal. The robustness of the anti-phase behavior is demonstrated by a noisy applied direct current Iext ∼ N(20, 12) at each simulation step Δt = 0.01. The remaining parameter values are identical.

Figure 8.

(A) Transient anti-phase dynamics of two subpopulations in a neuronal network model (Equation 17 and 19) composed of 50 all-to-all coupled interneurons with heterogeneous membrane capacitances Ci ∼ TN(1, 0.032, 0.91, 1.09), coupling (Equation 2), and initial conditions (Equation 3). This configuration leads to the rise of VHFOs in (B) the composed signal V = i=150Vi. Dynamics of (C) all neurons are zoomed, (D) composed, and depicted with periodograms of (E) individual neurons and (F) the composed signal. The robustness of the anti-phase behavior is demonstrated by a noisy applied direct current Iext ∼ N(20, 12) at each simulation step Δt = 0.01. The remaining parameter values are identical.

Close modal

Simulation of the Transient Collective Anti-Phase Behavior in a Neuronal Network of Coupled Interneurons

Consider the neuronal network (Equation 17 and 19) composed of 50 interacting interneurons formed into two clusters, 𝒫1 = {1, …, 25} and 𝒫2 = {26, …, 50}, respectively. Further, for the demonstration of the transient anti-phase collective behavior in this system, let the coupling strength take the form
εij=0.001,ij𝒫1×𝒫1𝒫2×𝒫2,0.0002,else,
(2)
where the first row refers to the coupling within each subpopulation, while the second one indicates a weaker connection between them.
Finally, assume that the clusters start in roughly opposite phases, specifically,
Vi0=40,i𝒫1,40,i𝒫2,hi0=0.25,ni0=0.5,i𝒫1𝒫2.
(3)

Figure 8 reveals the transient anti-phase behavior of these clusters. The periodograms suggest that this dynamics may result in the emergence of VHFOs in the simulation of the composed EEG signal. Namely, summing the signals of individual interneurons resulted in an oscillation with a frequency of ∼610 Hz. Moreover, although transient, the VHFOs may sustain for a relatively long time interval, significantly exceeding the duration of VHFOs observed in real EEG signals (Brázdil et al., 2017; Cimbalnik et al., 2020; Travnicek et al., 2021).

Transient UFRs and UFOs in a Two-Level Neuronal Network of Coupled Interneurons

Once more, consider the neuronal network (Equation 17 and 19) composed of four interacting interneurons formed into two clusters, 𝒫1 = {1, 2} and 𝒫2 = {3, 4}, respectively. Suppose that the network connectivity is given by the coupling
εij=0.0125,ij𝒫1×𝒫1𝒫2×𝒫2,0.0025,else;
(4)
see Figure 9. And finally, assume the following initial conditions:
Vi0=40,i13,40,i24,hi0=0.25,i𝒫1,0.27,i𝒫2,ni0=0.50,i𝒫1,0.52,i𝒫2.
(5)
Figure 9.

Network topology leading to apparent UFRs; see Figure 10.

Figure 9.

Network topology leading to apparent UFRs; see Figure 10.

Close modal

As one can see in Figure 10, this configuration leads to the emergence of transient UFRs in the composed signal V = i=14Vi. Specifically, in spite of the noisy applied direct current, the signal simulation contains a segment reporting a frequency of ∼1335 Hz. Despite being transient, the simulated UFRs exhibit a prolonged duration compared to the UFRs and UFOs observed in actual EEG records (Brazdil et al., 2023; Travnicek et al., 2021).

Figure 10.

(A) Transient anti-phase dynamics of two subpopulations in a neuronal network model (Equation 17 and 19) composed of four all-to-all coupled interneurons with heterogeneous membrane capacitances C1 = 0.998, C2 = 0.999, C3 = 1, C4 = 1.001, coupling (Equation 4), and initial conditions (Equation 5). This configuration results in the emergence of UFRs in (B) the composed signal V = i=14Vi. Dynamics of (C) all neurons are zoomed, (D) composed, and depicted with periodograms of (E) individual neurons and (F) the composed signal. The robustness of the anti-phase behavior is demonstrated by a noisy applied direct current Iext ∼ N(24, 12) at each simulation step Δt = 0.01. The remaining parameter values are identical.

Figure 10.

(A) Transient anti-phase dynamics of two subpopulations in a neuronal network model (Equation 17 and 19) composed of four all-to-all coupled interneurons with heterogeneous membrane capacitances C1 = 0.998, C2 = 0.999, C3 = 1, C4 = 1.001, coupling (Equation 4), and initial conditions (Equation 5). This configuration results in the emergence of UFRs in (B) the composed signal V = i=14Vi. Dynamics of (C) all neurons are zoomed, (D) composed, and depicted with periodograms of (E) individual neurons and (F) the composed signal. The robustness of the anti-phase behavior is demonstrated by a noisy applied direct current Iext ∼ N(24, 12) at each simulation step Δt = 0.01. The remaining parameter values are identical.

Close modal
Moreover, let us extend the previous neuronal network model by a third pair of interneurons. Specifically, consider six interacting interneurons organized into three clusters, 𝒫1 = {1, 2}, 𝒫2 = {3, 4}, and 𝒫3 = {5, 6}, respectively; assume that the network connectivity is determined by the coupling
εij=0.01,ij𝒫1×𝒫1𝒫2×𝒫2𝒫3×𝒫3,0.001,else;
(6)
see Figure 11. Finally, let the initial state of the system satisfy
Vi0=40,i135,40,i246,hi0=0.25,i𝒫1,0.26,i𝒫2,0.27,i𝒫3,ni0=0.50,i𝒫1,0.51,i𝒫2,0.52,i𝒫3,
(7)
characterizing anti-phase synchrony within each cluster and slightly shifted phases between clusters.
Figure 11.

Network topology leading to apparent UFOs; see Figure 12.

Figure 11.

Network topology leading to apparent UFOs; see Figure 12.

Close modal

As Figure 12 shows, this setting results in the emergence of UFOs in the composed signal V = i=16Vi. Namely, despite the noisy applied direct current, after a transient period, the simulated signal composition oscillates with a frequency of ∼2 kHz.

Figure 12.

Stable dynamics of three subpopulations in a neuronal network model (Equation 17 and 19) comprising six all-to-all coupled interneurons with heterogeneous membrane capacitances Ci ∼ TN(1, 0.0012, 0.997, 1.003), coupling (Equation 6), and initial conditions (Equation 7). This configuration results in the birth of UFOs in (A), (B) the composed signal V = i=16Vi. Dynamics of (C) all neurons are shown in detail, (D) composed, and depicted with periodograms of (E) individual neurons and (F) the composed signal. The robustness of the regular shift-phase behavior is demonstrated by a noisy applied direct current Iext ∼ N(24, 0.032) at each simulation step Δt = 0.01. The remaining parameter values are identical.

Figure 12.

Stable dynamics of three subpopulations in a neuronal network model (Equation 17 and 19) comprising six all-to-all coupled interneurons with heterogeneous membrane capacitances Ci ∼ TN(1, 0.0012, 0.997, 1.003), coupling (Equation 6), and initial conditions (Equation 7). This configuration results in the birth of UFOs in (A), (B) the composed signal V = i=16Vi. Dynamics of (C) all neurons are shown in detail, (D) composed, and depicted with periodograms of (E) individual neurons and (F) the composed signal. The robustness of the regular shift-phase behavior is demonstrated by a noisy applied direct current Iext ∼ N(24, 0.032) at each simulation step Δt = 0.01. The remaining parameter values are identical.

Close modal

Simulation of the Collective Anti-Phase Behavior in a Model of Destexhe–Paré Neuronal Network

Let us take into account the neuronal network (Equation 20 and 21) comprised of 50 coupled Destexhe–Paré neurons formed into two clusters, 𝒫1 = {1, …, 25} and 𝒫2 = {26, …, 50}, respectively. Unlike previous cases, here we assume that the distribution of capacitances differs in the mean value between these two groups. Namely, let
CiTN0.950.0320.861.04,i𝒫1,andCiTN10.0320.911.09,i𝒫2.
(8)
Let the network connectivity take the form
εij=0.005,ij𝒫1×𝒫1𝒫2×𝒫2,0.001,else,
(9)
characterizing the coupling strength within each subpopulation and between them, respectively.
And finally, for the demonstration of transient anti-phase behavior in this network, suppose that the clusters start in roughly opposite phases; specifically, let
Vi0=75,i𝒫1,35,i𝒫2,mi0=0.5,hi0=0.2,ni0=0.4,mM,i0=0.24,i𝒫1𝒫2.
(10)

The transient anti-phase behavior of these clusters is presented in Figure 13. At the bottom, one can see signal segments with the examined dynamics. The periodograms corresponding to these segments demonstrate that the VHFOs are present in the composed EEG signal simulation. To be more specific, the combination of signals from each neuron leads to an oscillation with a frequency of ∼740 Hz (fI curves for a single neuron for various values of C are provided in the Supporting Information (see Figure S5); the default parameter choice of Iext = 40 leads to a frequency of 360 Hz). Although transient, the VHFOs can persist for a relatively long period, significantly surpassing the length of VHFOs detected in actual EEG signals (Brázdil et al., 2017; Brazdil et al., 2023; Cimbalnik et al., 2020; Travnicek et al., 2021).

Figure 13.

(A) Transient anti-phase dynamics of subpopulations 𝒫1, 𝒫2 in a neuronal network model Equation 20 and 21 composed of 50 all-to-all coupled Destexhe–Paré neurons with heterogeneous membrane capacitances (Equation 8), coupling (Equation 9), and initial conditions (Equation 10). Dynamics of (B) composed signal, zoomed dynamics of (C) individual neurons, and (D) composed signal are depicted with periodograms of (E) individual neurons and (F) the composed signal. The robustness of the anti-phase behavior is demonstrated by a noisy applied direct current Iext ∼ N(40, 12) at each simulation step Δt = 0.01. The remaining parameter values are identical.

Figure 13.

(A) Transient anti-phase dynamics of subpopulations 𝒫1, 𝒫2 in a neuronal network model Equation 20 and 21 composed of 50 all-to-all coupled Destexhe–Paré neurons with heterogeneous membrane capacitances (Equation 8), coupling (Equation 9), and initial conditions (Equation 10). Dynamics of (B) composed signal, zoomed dynamics of (C) individual neurons, and (D) composed signal are depicted with periodograms of (E) individual neurons and (F) the composed signal. The robustness of the anti-phase behavior is demonstrated by a noisy applied direct current Iext ∼ N(40, 12) at each simulation step Δt = 0.01. The remaining parameter values are identical.

Close modal

The primary objective of our paper was to establish a rigorous basis for the possibility that a reduction in the strength of the connections between neurons in an epileptic focus could produce high, very high-frequency, and ultra-fast oscillations measured at microelectrodes in its vicinity. Although such oscillations are biomarkers used for diagnostic purposes in presurgical evaluation, there is still no definitive way to differentiate between physiological (most commonly associated with the normal physiological function of cognitive or sensory processing (see Pail et al., 2020) and pathological HFOs in EEG signals (Frauscher et al., 2018). There are recent studies (Cimbalnik et al., 2018, 2020; Nejedly et al., 2019) that introduced algorithms for the detection of pathological HFOs using pathological events visually marked by expert reviewers inside the seizure onset zone of patients with focal epilepsy (detection training and machine learning), but finding a way to distinguish between physiological and pathological HFOs remains challenging. In contrast, VHFOs and UFOs are likely to be important markers of epileptogenicity (Brázdil et al., 2017), and a plausible mechanistic explanation of such high frequencies in the LFP signal adds to the understanding of the pathophysiology and improving the treatment of drug-resistant focal epilepsy or the development of new less invasive tools for its treatment.

We established that synchronization of a weakly coupled neuronal network occurs in a large parameter region. The phenomenon of multistability allows for the synchronization of in-phase, anti-phase, or other spatiotemporally symmetric phase shifts. This enables the network to generate tonic spiking at its own frequency and simultaneously produce spikes at a regularly shifted phase. This phenomenon potentially creates an event that can be recorded as a very high-frequency or ultra-fast oscillation by nearby electrodes. One assumes the recorded signal is a weighted sum of the nearby activity (see Katzner et al., 2009; Mitzdorf, 1985), but here, for simplicity, all neurons contribute equally as this choice does not affect the dominant frequency. We simulated the network for different parameters in physiological ranges and subjected the system to an external current with random noise around a reference value. We found that synchronized VHFOs and UFOs are possible for very weakly coupled nearly identical neurons, persisting even with small random noise affecting the magnitude of the external current input. Our findings resonate with prior work highlighting already fast ripples as a network phenomenon in epileptic rodents, driven by out-of-phase synchronized firing across neuron pools (see Jefferys et al., 2012).

In this study, we implemented gap-junctional coupling into our model because the reactions of the dynamics would manifest more slowly if chemical synapses were utilized, particularly in the context of modeling very high-frequency oscillations. Gap junctions, as electrical synapses, are known for their ability to support rapid and bidirectional communication between neurons, which makes them suitable for modeling high-frequency oscillatory dynamics. However, we acknowledge the intricate, multifaceted nature of neuronal networks and therefore concede that both gap junctions and chemical synapses might play significant roles in the overall dynamic behavior. Nevertheless, we have shown that the phenomenon of frequency multiplication is generic and model-independent. There are other, more complex models with gap-junctional coupling but random (Erdös–Renyi graph) connectivity that could show the same (see Ghosh et al., 2020). Initial simulations showed emergent network oscillations with higher frequencies (up to 70 Hz instead of the individual 10 Hz) as the coupling strength increased. We observed that the network fired in clusters, but neurons did not always participate in the active cluster. So, while our mechanism seems relevant, a complete characterization requires a separate study. Future studies might benefit from models that concurrently integrate both types of synapses to capture a broader range of neuronal dynamics and interactions.

We used gap-junctional coupling of near-identical neurons of the same type for our analysis, but other parameters can also be varied with similar results. Specifically, we have already examined rigorously that slight changes in the sodium conductance gNa and reversal potential VNa in the interneuron model do not affect the obtained results. This will be reported elsewhere. Moreover, we found that excitation with a noisy external current gives analogous results. The approach is applicable to a multilevel network comprising subgroups of neurons that already exhibit VHFOs through anti-phase spiking within each subgroup. As a result, the presence of weak coupling between these subgroups gives rise to transient UFRs and UFOs in the simulated EEG signal. However, there is still much to be explored, including whether VHFOs and UFOs can be simulated in the EEG signal for a more complex network composed of multiple different types of neurons and how the network connectivity matrix affects the dynamics. Moreover, the persisting stability of anti-phase oscillations in a multilevel network is currently undetermined, as their stability and resilience have been established in single-level and two-level networks with very weak coupling. However, it is likely that the stability region will diminish in the multilevel scenario.

Our results support the hypothesis of a spatially limited pathological phenomenon, which is also manifested in VHFOs and UFOs. In addition, the coupling type and its strength can affect the dynamics of the coupled neurons, with strong enough coupling leading to in-phase synchrony or small phase-shift synchrony and very weak coupling allowing stable anti-phase regime, that may, in that case, temporally appear in the EEG signal. This result is consistent with the observation of altered connectivity in epileptic brain networks (see DeSalvo, Douw, Tanaka, Reinsberger, & Stufflebeam, 2014).

In conclusion, our study provides a mathematical framework for investigating the occurrence of HFOs, VHFOs, and UFOs in the EEG signals of patients with focal epilepsy. Our findings suggest that a reduction in the strength of the connections between neurons in an epileptic focus can produce very high-frequency signals, which could aid in understanding the pathology and improving the treatment of focal epilepsy. However, there is still much to be explored regarding the coupling types, their strength, and other parameters that could affect the dynamics of the coupled neurons.

Morris–Lecar Neuron Model

The Morris–Lecar model was developed by Morris and Lecar (1981) to reproduce the variety of oscillatory behavior in relation to Ca2+ and K+ conductance (van Putten, 2020). Although it was originally created to model the potential in the muscle fiber of the giant barnacle, the Morris–Lecar model is today used for modeling brain pyramidal neurons due to its simplicity and ability to model various dynamics (Lecar, 2007; Prescott, De Koninck, & Sejnowski, 2008; Prescott, Ratté, et al., 2008). The model is described by the following two-dimensional system of differential equations
CddtV=IextgLVVLgCamVVVCagKwVVK,ddtw=wVwτwV,
(11)
where V is the membrane potential [mV] and w represents the activation variable for K+, that is, the probability that the K+ channel is conducting. The parameter C corresponds to the cell membrane capacitance.
In comparison to the K+ current, the Ca2+ current changes rapidly, and its activation is assumed to be instantaneous. The activation function is modeled as
mV=121+tanhVβ1β2,
(12)
where β1 and β2 denote the potential at which the Ca2+ current is half-activated and the slope of the activation voltage dependence, respectively.
Similarly, for the voltage-dependent steady state w, we use the form
wV=121+tanhVβ3β4,
(13)
where β3 and β4 have analogous meanings to β1 and β2, respectively. The timescale of the w variable is given by
τwV=12φcoshVβ32β41,
(14)
where the parameter φ is a temperature factor adjusting the relative timescale of V and w.

The system parameters gL, gCa, and gK represent the maximum leak, Ca2+, and K+ electrical conductances through membrane channels, respectively, and VL, VCa, and VK are reversal potentials of the specific ion channels. Finally, Iext denotes the externally applied direct current, whichdirect can be used as a single parameter determining a class of the studied cell, that is, whether we deal with class I or class II neuron according to the classification proposed by Hodgkin (1948).

Unless otherwise stated, we make the parameter choices as listed in Table 1. These physiological parameter values were taken from Gutkin and Ermentrout (1998). Characteristic values for the hippocampal cell membrane capacitance C are found to be around 1 μF/cm2 (Destexhe & Paré, 1999; Ermentrout & Terman, 2010; Hartveit, Veruki, & Zandt, 2019; Keener & Sneyd, 2009; Tewari et al., 2018). Therefore, within the framework of this article, we consider these values.

Model of Morris–Lecar Neuronal Network

Timescales of synapses of neurons and gap junctions differ. Specifically, synaptic coupling refers to the mechanism whereby neurotransmitters are released from a presynaptic neuron and subsequently bind to receptors located on a postsynaptic neuron, inducing changes in the latter’s membrane potential. This process is relatively slow, taking several milliseconds to occur. In contrast, gap junctions provide direct electrical communication between neurons, allowing for very rapid transmission of electrical signals with virtually no delay in the communication between the coupled neurons. Given our objective to model high frequencies, we opted for a network model that incorporates gap junctions, which are known to be present in the mammalian central nervous system, including the hippocampus (Draguhn, Traub, Schmitz, & Jefferys, 1998; Fukuda & Kosaka, 2000). Morris–Lecar neuronal networks with gap junctions have recently been studied with regard to epilepsy (Volman, Perc, & Bazhenov, 2011; Naze, Bernard, & Jirsa, 2015). HFOs have also been reported in a computational study with Hodgkin–Huxley type models with gap junctions by Helling, Koppert, Visser, and Kalitzin (2015). Here, we examine a small population of neurons with gap junction coupling as it has been found that the epileptic brain tissue may be functionally isolated from surrounding brain regions (see Klimes et al., 2016; Warren et al., 2010). Moreover, the review of Jin and Zhong (2011) focuses on gap junctions in the pathophysiology of epilepsy, describing their role in the generation, synchronization, and maintenance of seizures, and points to gap junctions as promising targets for the development of antiseizure medication, emphasizing the need for further research in this field to refine our understanding and treatment of epilepsy.

Let us examine a neuronal network model of N electrically coupled Morris–Lecar neurons (Equation 11) described by
CiddtVi=IextgLViVLgCamViViVCagKwiViVKj=1NεijViVj,ddtwi=wViwiτwVi,
(15)
where i ∈ {1, …, N} and K = εiji,j=1N represents a coupling matrix determining the network topology, that is, its connectivity. Within this article, concerning the Morris–Lecar model, we focus on the case of an all-to-all coupled neuronal network with identical coupling strength
εij=ε=0.1N.
(16)
The functions m(V), w(V), and τw(V) are defined in Equation 12, 13, and 14, respectively. Initial conditions for gating variable: w1(0) = w2(0) = 0.04 unless specified otherwise.

Interneuron Model

Consider a neuronal network composed of N electrically coupled interneurons (White, Chow, Rit, Soto-Treviño, & Kopell, 1998). Let the membrane potential Vi [mV] of neuron i satisfy the equation
CiddtVi=IextILINaIKj=1NεijViVj,
(17)
where Ci represents the membrane capacitance, Iext is the externally applied current, and IL, INa, and IK stand for the Hodgkin–Huxley type leak, Na+, and K+ currents, respectively. The last term represents the gap-junctional coupling between neuron i and all other neurons in the network, where εij is the gap junction conductance. The Na+ and K+ currents are given by
INa=gNami3hiViVNaandIK=gKni4ViVK,
(18)
where gNa and gK are the maximum conductances, mi and hi represent gating variables for Na+ channels, and ni is a gating variable for K+ channels. The leak current is given by IL = gL(ViVL), where gL stands for the leak conductance and VL is the leak reversal potential.
The activation variable mi is assumed to tend rapidly to its steady state. Hence one can substitute it by its asymptotic value mi = m(Vi) = (1 + exp [−0.08(Vi + 26)])−1. Further, the gating variables hi and ni obey
ddthi=hVihiτhViandddtni=nViniτnVi
(19)
with
hV=11+exp0.13V+38,τhV=0.61+exp0.12V+67,
nV=11+exp0.045V+10,τnV=0.5+21+exp0.045V50,
respectively.

In this article, we consider parameter values C ≈ 1 μS/cm2, gL = 0.1 mS/cm2, gNa = 30 mS/cm2, gK = 20 mS/cm2, VL = −60 mV, VNa = 45 mV, and VK = −80 mV. These values are within physiological ranges and give the high-frequency firing rates typical of hippocampal interneurons (White et al., 1998).

Destexhe–Paré Neuron Model

The last model of a hippocampal pyramidal cell we consider within this paper comes from Destexhe and Paré (1999). Currents dynamics are described by Hodgkin–Huxley type models with kinetics based on a model of hippocampal pyramidal cells from Traub and Miles (1991), adjusted to match voltage-clamp data of cortical pyramidal cells (Huguenard, Hamill, & Prince, 1988), and a noninactivating K+ current was described in Mainen, Joerges, Huguenard, and Sejnowski (1995).

Once again, we examine a neuronal network comprising N coupled neurons. Let the membrane potential Vi [mV] of neuron i be described by equation
CiddtVi=IextILINaIKdrIMj=1NεijViVj,
(20)
where Ci denotes the membrane capacitance; Iext is the externally applied direct current, and IL, INa, IKdr, and IM represent the leak, Na+, and K+ currents. The sum term stands for the gap junction coupling between neuron i and all other neurons in the network with the gap junction conductance εij. As usual, the leak current is given by IL = gL(ViVL), where gL is the leak conductance, and VL denotes the leak reversal potential. As mentioned, all gating variables p ∈ {m, h, n, mM} obey first-order kinetics
ddtpi=αpVi1piβpVipi
(21)
with functions αp(Vi) and βp(Vi) taking the form specific for each variable.
The voltage-dependent Na+ current was described by Traub and Miles (1991) and is given by
INa=gNami3hiViVNa,
(22)
where gNa stands for its maximum conductance, VNa is its reversal potential, and mi and hi represent gating variables satisfying Equation 21 with
αmV=0.32VVT13expVVT13/41,βmV=0.28VVT40expVVT40/51,
αhV=0.128expVVTVS17/18,βhV=41+expVVTVS40/5
for VT = −58 mV and VS = −10 mV fitting the voltage-clamp data.
The current corresponding to the delayed rectifier K+ channel was studied by Traub and Miles (1991) and is given by
IKdr=gKdrni4ViVK,
(23)
where gKdr is the maximum conductance, VK represents its reversal potential, and ni obeys Equation 21 with
αnV=0.032VVT15expVVT15/51,βnV=0.5expVVT10/40.
And finally, the noninactivating current was described in Mainen et al. (1995) and takes the form
IM=gMmM,iViVK,
(24)
where gM denotes the maximum conductance and mM obeys Equation 21 with
αmMV=0.0001V+301expV+30/9,βmMV=0.0001V+301expV+30/9.

Within this paper, we consider the following physiological parameter values for a hippocampal pyramidal cell taken from Destexhe and Paré (1999). Specifically, C ≈ 1 μS/cm2, gL = 0.019 mS/cm2, gNa = 120 mS/cm2, gKdr = 100 mS/cm2, gM = 2 mS/cm2, VL = −65 mV, VNa = 55 mV, and VK = − 85 mV.

We used the standard ode45 solver in matlab for deterministic settings, while for simulations with noisy input, we add noise to the external current at each simulation step with Δt = 0.01 using the Euler–Maruyama method. Bifurcation curves have been computed using the numerical bifurcation package MatCont (Dhooge, Govaerts, & Kuznetsov, 2003; Dhooge, Govaerts, Kuznetsov, Meijer, & Sautois, 2008).

The authors thank Pavel Jurák and Jan Cimbálník for their insightful comments and critical reading of the work before submission.

Supporting information for this article is available at https://doi.org/10.1162/netn_a_00351. Code is available at https://gitlab.ics.muni.cz/9607/UFOs.

Lenka Přibylová: Conceptualization; Formal analysis; Funding acquisition; Methodology; Project administration; Supervision; Validation; Writing – original draft. Jan Ševčík: Data curation; Formal analysis; Methodology; Software; Visualization; Writing – original draft. Veronika Eclerová: Conceptualization; Methodology; Validation; Writing – review & editing. Petr Klimeš: Validation; Writing – review & editing. Milan Brázdil: Supervision; Validation; Writing – review & editing. Hil Meijer: Methodology; Supervision; Validation; Writing – review & editing.

Lenka Přibylová, Masarykova Univerzita (https://dx.doi.org/10.13039/501100010653), Award ID: MUNI/G/1213/2022. Jan Ševčík, Masarykova Univerzita (https://dx.doi.org/10.13039/501100010653), Award ID: MUNI/G/1213/2022. Veronika Eclerová, Masarykova Univerzita (https://dx.doi.org/10.13039/501100010653), Award ID: MUNI/G/1213/2022. Milan Brázdil, Masarykova Univerzita (https://dx.doi.org/10.13039/501100010653), Award ID: MUNI/G/1213/2022, Grantová Agentura České Republiky (https://dx.doi.org/10.13039/501100001824), Award ID: GA22-28784S, and European Union – Next Generation EU, Award ID: LX22NPO5107 (MEYS). Petr Klimeš, Grantová Agentura České Republiky (https://dx.doi.org/10.13039/501100001824), Award ID: GA22-28784S, and European Union – Next Generation EU, Award ID: LX22NPO5107 (MEYS).

Hopf bifurcation (H):

A change in a dynamical system that is related to cycle birth from equilibrium.

Limit point of cycles (LPC):

A point in a dynamical system representing the fold of a sequence of periodic orbits or cycles as a parameter changes.

Phase-locking:

A state in a dynamical system where oscillators synchronize their frequencies, maintaining a constant relative phase difference over time.

Local field potentials:

Electrophysiological signals reflecting the summed electrical activity of a group of neurons in a local area.

Rigid:

In the context of dynamical systems, refers to behaviors that do not vary or change under a set of conditions.

Phase-shift synchrony:

A synchronous state where oscillating systems align their cycles but with a constant time delay or phase difference.

Saddle-node on an invariant curve (SNIC) bifurcation:

A critical event in a dynamical system where a stable limit cycle is created or destroyed alongside an equilibrium point.

Saddle-node equilibrium (LP):

A state in a dynamical system where two equilibrium points, a stable and an unstable one, merge and annihilate each other.

Subcritical Hopf bifurcation:

A system transition where an unstable limit cycle exists before the bifurcation point and a stable equilibrium point loses stability after it.

Bifurcation:

A critical point in a dynamical system where a small change in system parameter values results in a qualitative change in its behavior.

Two-dimensional invariant torus:

A set of solutions in a dynamical system that form a toroidal shape, wherein trajectories neither converge nor diverge but remain on the torus.

Quasi-periodic:

A behavior in a dynamical system with oscillations at two or more incommensurable (irrational ratio) frequencies, creating a nonrepeating pattern.

Period-doubling bifurcation (PD):

A transition in a dynamical system where a small change in parameters causes the period of oscillations to double, often preceding chaos.

Arnold tongues:

Regions in parameter space indicating synchronization or frequency-locking ranges in dynamical systems, typically with a characteristic ’tongue’ shape.

Synaptic coupling:

The process where activity in one neuron influences another through synaptic transmission, facilitating communication and synchronization between neurons.

Gap junction:

A specialized intercellular connection that allows direct electrical and chemical communication between neighboring cells, facilitating synchronization of activity.

Abrams
,
D. M.
, &
Strogatz
,
S. H.
(
2004
).
Chimera states for coupled oscillators
.
Physical Review Letters
,
93
(
17
),
174102
. ,
[PubMed]
Aihara
,
K.
,
Takabe
,
T.
, &
Toyoda
,
M.
(
1990
).
Chaotic neural networks
.
Physics Letters A
,
144
(
6–7
),
333
340
.
Bakhtiari
,
S.
,
Manshadi
,
M. K.
,
Candas
,
M.
, &
Beskok
,
A.
(
2023
).
Changes in electrical capacitance of cell membrane reflect drug partitioning-induced alterations in lipid bilayer
.
Micromachines
,
14
(
2
),
316
. ,
[PubMed]
Boccaletti
,
S.
,
Pisarchik
,
A. N.
,
Del Genio
,
C. I.
, &
Amann
,
A.
(
2018
).
Synchronization: From coupled systems to complex networks
.
Cambridge University Press
.
Brázdil
,
M.
,
Pail
,
M.
,
Halámek
,
J.
,
Plešinger
,
F.
,
Cimbálník
,
J.
,
Roman
,
R.
, …
Jurák
,
P.
(
2017
).
Very high-frequency oscillations: Novel biomarkers of the epileptogenic zone
.
Annals of Neurology
,
82
(
2
),
299
310
. ,
[PubMed]
Brazdil
,
M.
,
Worrell
,
G. A.
,
Travnicek
,
V.
,
Pail
,
M.
,
Roman
,
R.
,
Plesinger
,
F.
, …
Jurák
,
P.
(
2023
).
Ultra fast oscillations in the human brain and their functional significance
.
medRxiv
.
Calim
,
A.
,
Hövel
,
P.
,
Ozer
,
M.
, &
Uzuntarla
,
M.
(
2018
).
Chimera states in networks of type-I Morris–Lecar neurons
.
Physical Review E
,
98
(
6
),
062217
.
Cessac
,
B.
, &
Samuelides
,
M.
(
2007
).
From neuron to neural networks dynamics
.
European Physical Journal Special Topics
,
142
(
1
),
7
88
.
Chow
,
C. C.
, &
Kopell
,
N.
(
2000
).
Dynamics of spiking neurons with electrical coupling
.
Neural Computation
,
12
(
7
),
1643
1678
. ,
[PubMed]
Chowdhury
,
S. N.
,
Rakshit
,
S.
,
Buldu
,
J. M.
,
Ghosh
,
D.
, &
Hens
,
C.
(
2021
).
Antiphase synchronization in multiplex networks with attractive and repulsive interactions
.
Physical Review E
,
103
(
3
),
032310
. ,
[PubMed]
Cimbalnik
,
J.
,
Brinkmann
,
B.
,
Kremen
,
V.
,
Jurak
,
P.
,
Berry
,
B.
,
Gompel
,
J. V.
, …
Worrell
,
G.
(
2018
).
Physiological and pathological high frequency oscillations in focal epilepsy
.
Annals of Clinical and Translational Neurology
,
5
(
9
),
1062
1076
. ,
[PubMed]
Cimbalnik
,
J.
,
Pail
,
M.
,
Klimes
,
P.
,
Travnicek
,
V.
,
Roman
,
R.
,
Vajcner
,
A.
, &
Brazdil
,
M.
(
2020
).
Cognitive processing impacts high frequency intracranial EEG activity of human hippocampus in patients with pharmacoresistant focal epilepsy
.
Frontiers in Neurology
,
11
,
578571
. ,
[PubMed]
DeSalvo
,
M. N.
,
Douw
,
L.
,
Tanaka
,
N.
,
Reinsberger
,
C.
, &
Stufflebeam
,
S. M.
(
2014
).
Altered structural connectome in temporal lobe epilepsy
.
Radiology
,
270
(
3
),
842
848
. ,
[PubMed]
Destexhe
,
A.
, &
Paré
,
D.
(
1999
).
Impact of network activity on the integrative properties of neocortical pyramidal neurons in vivo
.
Journal of Neurophysiology
,
81
(
4
),
1531
1547
. ,
[PubMed]
Dhooge
,
A.
,
Govaerts
,
W.
, &
Kuznetsov
,
Y. A.
(
2003
).
MatCont: A MATLAB package for numerical bifurcation analysis of ODEs
.
ACM Transactions on Mathematical Software
,
29
(
2
),
141
164
.
Dhooge
,
A.
,
Govaerts
,
W.
,
Kuznetsov
,
Y. A.
,
Meijer
,
H. G.
, &
Sautois
,
B.
(
2008
).
New features of the software MatCont for bifurcation analysis of dynamical systems
.
Mathematical and Computer Modelling of Dynamical Systems
,
14
(
2
),
147
175
.
Draguhn
,
A.
,
Traub
,
R.
,
Schmitz
,
D.
, &
Jefferys
,
J.
(
1998
).
Electrical coupling underlies high-frequency oscillations in the hippocampus in vitro
.
Nature
,
394
(
6689
),
189
192
. ,
[PubMed]
Duan
,
L.
, &
Lu
,
Q.
(
2006
).
Codimension-two bifurcation analysis on firing activities in Chay neuron model
.
Chaos, Solitons & Fractals
,
30
(
5
),
1172
1179
.
Ermentrout
,
B.
, &
Terman
,
D. H.
(
2010
).
Mathematical foundations of neuroscience
(
Vol. 35
).
Springer
.
Frauscher
,
B.
,
von Ellenrieder
,
N.
,
Zelmann
,
R.
,
Rogers
,
C.
,
Nguyen
,
D. K.
,
Kahane
,
P.
, …
Gotman
,
J.
(
2018
).
High-frequency oscillations in the normal human brain
.
Annals of Neurology
,
84
(
3
),
374
385
. ,
[PubMed]
Fukuda
,
T.
, &
Kosaka
,
T.
(
2000
).
Gap junctions linking the dendritic network of gabaergic interneurons in the hippocampus
.
Journal of Neuroscience
,
20
(
4
),
1519
1528
. ,
[PubMed]
Gabbiani
,
F.
, &
Cox
,
S. J.
(
2017
).
Mathematics for neuroscientists
.
Academic Press
.
Gentet
,
L. J.
,
Stuart
,
G. J.
, &
Clements
,
J. D.
(
2000
).
Direct measurement of specific membrane capacitance in neurons
.
Biophysical Journal
,
79
(
1
),
314
320
. ,
[PubMed]
Ghosh
,
S.
,
Mondal
,
A.
,
Ji
,
P.
,
Mishra
,
A.
,
Dana
,
S. K.
,
Antonopoulos
,
C. G.
, &
Hens
,
C.
(
2020
).
Emergence of mixed mode oscillations in random networks of diverse excitable neurons: The role of neighbors and electrical coupling
.
Frontiers in Computational Neuroscience
,
14
,
49
. ,
[PubMed]
Golubitsky
,
M.
,
Messi
,
L. M.
, &
Spardy
,
L. E.
(
2016
).
Symmetry types and phase-shift synchrony in networks
.
Physica D: Nonlinear Phenomena
,
320
,
9
18
.
Golubitsky
,
M.
,
Romano
,
D.
, &
Wang
,
Y.
(
2012
).
Network periodic solutions: Patterns of phase-shift synchrony
.
Nonlinearity
,
25
(
4
),
1045
.
Golubitsky
,
M.
, &
Stewart
,
I.
(
2006
).
Nonlinear dynamics of networks: The groupoid formalism
.
Bulletin of the American Mathematical Society
,
43
(
3
),
305
364
.
Golubitsky
,
M.
, &
Stewart
,
I.
(
2016
).
Rigid patterns of synchrony for equilibria and periodic cycles in network dynamics
.
Chaos: An Interdisciplinary Journal of Nonlinear Science
,
26
(
9
),
094803
. ,
[PubMed]
Golubitsky
,
M.
,
Stewart
,
I.
, &
Török
,
A.
(
2005
).
Patterns of synchrony in coupled cell networks with multiple arrows
.
SIAM Journal on Applied Dynamical Systems
,
4
(
1
),
78
100
.
Gutkin
,
B. S.
, &
Ermentrout
,
G. B.
(
1998
).
Dynamics of membrane excitability determine interspike interval variability: A link between spike generation mechanisms and cortical spike train statistics
.
Neural Computation
,
10
(
5
),
1047
1065
. ,
[PubMed]
Hao
,
J.
,
Cui
,
Y.
,
Niu
,
B.
,
Yu
,
L.
,
Lin
,
Y.
,
Xia
,
Y.
, …
Guo
,
D.
(
2021
).
Roles of very fast ripple (500–1000 Hz) in the hippocampal network during status epilepticus
.
International Journal of Neural Systems
,
31
(
4
),
2150002
. ,
[PubMed]
Hartveit
,
E.
,
Veruki
,
M. L.
, &
Zandt
,
B.-J.
(
2019
).
Capacitance measurement of dendritic exocytosis in an electrically coupled inhibitory retinal interneuron: An experimental and computational study
.
Physiological Reports
,
7
(
15
),
e14186
. ,
[PubMed]
Helling
,
R.
,
Koppert
,
M.
,
Visser
,
G.
, &
Kalitzin
,
S.
(
2015
).
Gap junctions as common cause of high-frequency oscillations and epileptic seizures in a computational cascade of neuronal mass and compartmental modeling
.
International Journal of Neural Systems
,
25
(
6
),
1550021
. ,
[PubMed]
Hodgkin
,
A. L.
(
1948
).
The local electric changes associated with repetitive action in a non-medullated axon
.
Journal of Physiology
,
107
(
2
),
165
181
. ,
[PubMed]
Hoppensteadt
,
F. C.
, &
Izhikevich
,
E. M.
(
1997
).
Weakly connected neural networks
(
Vol. 126
).
Springer Science & Business Media
.
Huguenard
,
J.
,
Hamill
,
O.
, &
Prince
,
D.
(
1988
).
Developmental changes in Na+ conductances in rat neocortical neurons: Appearance of a slowly inactivating component
.
Journal of Neurophysiology
,
59
(
3
),
778
795
. ,
[PubMed]
Izhikevich
,
E.
(
2006
).
Dynamical systems in neuroscience
.
MIT Press
.
Jacobs
,
J.
,
LeVan
,
P.
,
Chander
,
R.
,
Hall
,
J.
,
Dubeau
,
F.
, &
Gotman
,
J.
(
2008
).
Interictal high-frequency oscillations (80–500 Hz) are an indicator of seizure onset areas independent of spikes in the human epileptic brain
.
Epilepsia
,
49
(
11
),
1893
1907
. ,
[PubMed]
Jacobs
,
J.
,
Staba
,
R.
,
Asano
,
E.
,
Otsubo
,
H.
,
Wu
,
J.
,
Zijlmans
,
M.
, …
Gotman
,
J.
(
2012
).
High-frequency oscillations (HFOs) in clinical epilepsy
.
Progress in Neurobiology
,
98
(
3
),
302
315
. ,
[PubMed]
Jefferys
,
J. G.
,
De La Prida
,
L. M.
,
Wendling
,
F.
,
Bragin
,
A.
,
Avoli
,
M.
,
Timofeev
,
I.
, &
da Silva
,
F. H. L.
(
2012
).
Mechanisms of physiological and epileptic HFO generation
.
Progress in Neurobiology
,
98
(
3
),
250
264
. ,
[PubMed]
Jin
,
M.-M.
, &
Zhong
,
C.
(
2011
).
Role of gap junctions in epilepsy
.
Neuroscience Bulletin
,
27
(
6
),
389
406
. ,
[PubMed]
Jiruska
,
P.
,
Alvarado-Rojas
,
C.
,
Schevon
,
C. A.
,
Staba
,
R.
,
Stacey
,
W.
,
Wendling
,
F.
, &
Avoli
,
M.
(
2017
).
Update on the mechanisms and roles of high-frequency oscillations in seizures and epileptic disorders
.
Epilepsia
,
58
(
8
),
1330
1339
. ,
[PubMed]
Katzner
,
S.
,
Nauhaus
,
I.
,
Benucci
,
A.
,
Bonin
,
V.
,
Ringach
,
D. L.
, &
Carandini
,
M.
(
2009
).
Local origin of field potentials in visual cortex
.
Neuron
,
61
(
1
),
35
41
. ,
[PubMed]
Kawamura
,
Y.
,
Nakao
,
H.
,
Arai
,
K.
,
Kori
,
H.
, &
Kuramoto
,
Y.
(
2010
).
Phase synchronization between collective rhythms of globally coupled oscillator groups: Noiseless nonidentical case
.
Chaos: An Interdisciplinary Journal of Nonlinear Science
,
20
(
4
),
043110
. ,
[PubMed]
Keener
,
J.
, &
Sneyd
,
J.
(
2009
).
Mathematical physiology: I: Cellular physiology
.
Springer
.
Klimes
,
P.
,
Duque
,
J.
,
Brinkmann
,
B.
,
Van Gompel
,
J.
,
Stead
,
M.
,
St. Louis
,
E.
, …
Worrell
,
G.
(
2016
).
The functional organization of human epileptic hippocampus
.
Journal of Neurophysiology
,
115
(
6
),
3140
3145
. ,
[PubMed]
Kobelevskiy
,
I.
(
2008
).
Bifurcation analysis of a system of Morris-Lecar neurons with time delayed gap junctional coupling
.
(Unpublished master’s thesis)
.
University of Waterloo
.
Kuznetsov
,
Y. A.
(
2023
).
Elements of applied bifurcation theory
, (4th ed.,
Vol. 112
).
New York
:
Springer
.
Lang
,
X.
,
Lu
,
Q.
, &
Kurths
,
J.
(
2010
).
Phase synchronization in noise-driven bursting neurons
.
Physical Review E
,
82
(
2
),
021909
. ,
[PubMed]
Lecar
,
H.
(
2007
).
Morris–Lecar model
.
Scholarpedia
,
2
(
10
),
1333
.
Liu
,
C.
,
Liu
,
X.
, &
Liu
,
S.
(
2014
).
Bifurcation analysis of a Morris-Lecar neuron model
.
Biological Cybernetics
,
108
(
1
),
75
84
. ,
[PubMed]
Macdonald
,
R. L.
(
1989
).
Antiepileptic drug actions
.
Epilepsia
,
30
,
S19
S28
. ,
[PubMed]
Mainen
,
Z. F.
,
Joerges
,
J.
,
Huguenard
,
J. R.
, &
Sejnowski
,
T. J.
(
1995
).
A model of spike initiation in neocortical pyramidal neurons
.
Neuron
,
15
(
6
),
1427
1439
. ,
[PubMed]
Majhi
,
S.
,
Bera
,
B. K.
,
Ghosh
,
D.
, &
Perc
,
M.
(
2019
).
Chimera states in neuronal networks: A review
.
Physics of Life Reviews
,
28
,
100
121
. ,
[PubMed]
Mitzdorf
,
U.
(
1985
).
Current source-density method and application in cat cerebral cortex: Investigation of evoked potentials and EEG phenomena
.
Physiological Reviews
,
65
(
1
),
37
100
. ,
[PubMed]
Morris
,
C.
, &
Lecar
,
H.
(
1981
).
Voltage oscillations in the barnacle giant muscle fiber
.
Biophysical Journal
,
35
(
1
),
193
213
. ,
[PubMed]
Nagai
,
K. H.
, &
Kori
,
H.
(
2010
).
Noise-induced synchronization of a large population of globally coupled nonidentical oscillators
.
Physical Review E
,
81
(
6
),
065202
. ,
[PubMed]
Naze
,
S.
,
Bernard
,
C.
, &
Jirsa
,
V.
(
2015
).
Computational modeling of seizure dynamics using coupled neuronal networks: Factors shaping epileptiform activity
.
PLOS Computational Biology
,
11
(
5
),
e1004209
. ,
[PubMed]
Neiman
,
A.
,
Schimansky-Geier
,
L.
,
Cornell-Bell
,
A.
, &
Moss
,
F.
(
1999
).
Noise-enhanced phase synchronization in excitable media
.
Physical Review Letters
,
83
(
23
),
4896
.
Nejedly
,
P.
,
Cimbalnik
,
J.
,
Klimes
,
P.
,
Plesinger
,
F.
,
Halamek
,
J.
,
Kremen
,
V.
, …
Jurák
,
P.
(
2019
).
Intracerebral EEG artifact identification using convolutional neural networks
.
Neuroinformatics
,
17
(
2
),
225
234
. ,
[PubMed]
Nicholson
,
G. M.
,
Blanche
,
T.
,
Mansfield
,
K.
, &
Tran
,
Y.
(
2002
).
Differential blockade of neuronal voltage-gated Na+ and K+ channels by antidepressant drugs
.
European Journal of Pharmacology
,
452
(
1
),
35
48
. ,
[PubMed]
Nijholt
,
E.
,
Rink
,
B.
, &
Sanders
,
J.
(
2019
).
Center manifolds of coupled cell networks
.
SIAM Review
,
61
(
1
),
121
155
.
Pail
,
M.
,
Cimbálník
,
J.
,
Roman
,
R.
,
Daniel
,
P.
,
Shaw
,
D. J.
,
Chrastina
,
J.
, &
Brázdil
,
M.
(
2020
).
High frequency oscillations in epileptic and non-epileptic human hippocampus during a cognitive task
.
Scientific Reports
,
10
(
1
),
18147
. ,
[PubMed]
Patel
,
D. C.
,
Tewari
,
B. P.
,
Chaunsali
,
L.
, &
Sontheimer
,
H.
(
2019
).
Neuron-glia interactions in the pathophysiology of epilepsy
.
Nature Reviews Neuroscience
,
20
(
5
),
282
297
. ,
[PubMed]
Perko
,
L.
(
2001
).
Differential equations and dynamical systems
(
Vol. 7
).
Springer Science & Business Media
.
Pikovsky
,
A.
,
Rosenblum
,
M.
, &
Kurths
,
J.
(
2001
).
Synchronization: A universal concept in nonlinear sciences
.
Cambridge University Press
.
Prescott
,
S. A.
,
De Koninck
,
Y.
, &
Sejnowski
,
T. J.
(
2008
).
Biophysical basis for three distinct dynamical mechanisms of action potential initiation
.
PLoS Computational Biology
,
4
(
10
),
e1000198
. ,
[PubMed]
Prescott
,
S. A.
,
Ratté
,
S.
,
De Koninck
,
Y.
, &
Sejnowski
,
T. J.
(
2008
).
Pyramidal neurons switch from integrators in vitro to resonators under in vivo-like conditions
.
Journal of Neurophysiology
,
100
(
6
),
3030
3042
. ,
[PubMed]
Purves
,
D.
,
Augustine
,
G. J.
,
Fitzpatrick
,
D.
,
Hall
,
W.
,
LaMantia
,
A.-S.
, &
White
,
L.
(
2019
).
Neuroscience
(6th ed.).
Oxford University Press
.
Řehulka
,
P.
,
Cimbalnik
,
J.
,
Pail
,
M.
,
Chrastina
,
J.
,
Hermanová
,
M.
, &
Brázdil
,
M.
(
2019
).
Hippocampal high frequency oscillations in unilateral and bilateral mesial temporal lobe epilepsy
.
Clinical Neurophysiology
,
130
(
7
),
1151
1159
. ,
[PubMed]
Sebek
,
M.
,
Kawamura
,
Y.
,
Nott
,
A. M.
, &
Kiss
,
I. Z.
(
2019
).
Anti-phase collective synchronization with intrinsic in-phase coupling of two groups of electrochemical oscillators
.
Philosophical Transactions of the Royal Society A
,
377
(
2160
),
20190095
. ,
[PubMed]
Shilnikov
,
A.
(
2012
).
Complete dynamical analysis of a neuron model
.
Nonlinear Dynamics
,
68
(
3
),
305
328
.
Skinner
,
F.
,
Zhang
,
L.
,
Velazquez
,
J. P.
, &
Carlen
,
P.
(
1999
).
Bursting in inhibitory interneuronal networks: A role for gap-junctional coupling
.
Journal of Neurophysiology
,
81
(
3
),
1274
1283
. ,
[PubMed]
Staba
,
R. J.
, &
Bragin
,
A.
(
2011
).
High-frequency oscillations and other electrophysiological biomarkers of epilepsy: Underlying mechanisms
.
Biomarkers in Medicine
,
5
(
5
),
545
556
. ,
[PubMed]
Strogatz
,
S. H.
(
2018
).
Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering
.
CRC Press
.
Tewari
,
B. P.
,
Chaunsali
,
L.
,
Campbell
,
S. L.
,
Patel
,
D. C.
,
Goode
,
A. E.
, &
Sontheimer
,
H.
(
2018
).
Perineuronal nets decrease membrane capacitance of peritumoral fast spiking interneurons in a model of epilepsy
.
Nature Communications
,
9
(
1
),
4724
. ,
[PubMed]
Traub
,
R. D.
, &
Miles
,
R.
(
1991
).
Neuronal networks of the hippocampus
(
Vol. 777
).
Cambridge University Press
.
Travnicek
,
V.
,
Jurak
,
P.
,
Cimbalnik
,
J.
,
Klimes
,
P.
,
Daniel
,
P.
, &
Brazdil
,
M.
(
2021
).
Ultra-fast oscillation detection in EEG signal from deep-brain microelectrodes
. In
2021 43rd annual international conference of the IEEE Engineering in Medicine & Biology Society (EMBC)
(pp.
265
268
). ,
[PubMed]
Tsumoto
,
K.
,
Kitajima
,
H.
,
Yoshinaga
,
T.
,
Aihara
,
K.
, &
Kawakami
,
H.
(
2006
).
Bifurcations in Morris–Lecar neuron model
.
Neurocomputing
,
69
(
4–6
),
293
316
.
van Putten
,
M.
(
2020
).
Dynamics of neural networks
.
Springer
.
Vasickova
,
Z.
,
Klimes
,
P.
,
Cimbalnik
,
J.
,
Travnicek
,
V.
,
Pail
,
M.
,
Halamek
,
J.
, …
Brazdil
,
M.
(
2023
).
Shadows of very high-frequency oscillations can be detected in lower frequency bands of routine stereoelectroencephalography
.
Scientific Reports
,
13
(
1
),
1065
. ,
[PubMed]
Volman
,
V.
,
Perc
,
M.
, &
Bazhenov
,
M.
(
2011
).
Gap junctions and epileptic seizures—Two sides of the same coin?
PLoS One
,
6
(
5
),
e20572
. ,
[PubMed]
Wang
,
B.
,
Ke
,
W.
,
Guang
,
J.
,
Chen
,
G.
,
Yin
,
L.
,
Deng
,
S.
, …
Shu
,
Y.
(
2016
).
Firing frequency maxima of fast-spiking neurons in human, monkey, and mouse neocortex
.
Frontiers in Cellular Neuroscience
,
10
,
239
. ,
[PubMed]
Warren
,
C.
,
Hu
,
S.
,
Stead
,
M.
,
Brinkmann
,
B.
,
Bower
,
M.
, &
Worrell
,
G.
(
2010
).
Synchrony in normal and focal epileptic brain: The seizure onset zone is functionally disconnected
.
Journal of Neurophysiology
,
104
(
6
),
3530
3539
. ,
[PubMed]
White
,
J. A.
,
Chow
,
C. C.
,
Rit
,
J.
,
Soto-Treviño
,
C.
, &
Kopell
,
N.
(
1998
).
Synchronization and oscillatory dynamics in heterogeneous, mutually inhibited neurons
.
Journal of Computational Neuroscience
,
5
(
1
),
5
16
. ,
[PubMed]
Wiggins
,
S.
(
2003
).
Introduction to applied nonlinear dynamical systems and chaos
(2nd ed.).
New York
:
Springer
.
Worrell
,
G.
, &
Gotman
,
J.
(
2011
).
High-frequency oscillations and other electrophysiological biomarkers of epilepsy: Clinical studies
.
Biomarkers in Medicine
,
5
(
5
),
557
566
. ,
[PubMed]
Xing
,
M.
,
Song
,
X.
,
Wang
,
H.
,
Yang
,
Z.
, &
Chen
,
Y.
(
2022
).
Frequency synchronization and excitabilities of two coupled heterogeneous Morris-Lecar neurons
.
Chaos, Solitons & Fractals
,
157
,
111959
.
Zhou
,
C.
, &
Kurths
,
J.
(
2002
).
Noise-induced phase synchronization and synchronization transitions in chaotic oscillators
.
Physical Review Letters
,
88
(
23
),
230602
. ,
[PubMed]
Zhou
,
C.
, &
Kurths
,
J.
(
2003
).
Noise-induced synchronization and coherence resonance of a Hodgkin-Huxley model of thermally sensitive neurons
.
Chaos: An Interdisciplinary Journal of Nonlinear Science
,
13
(
1
),
401
409
. ,
[PubMed]
Zijlmans
,
M.
,
Jacobs
,
J.
,
Zelmann
,
R.
,
Dubeau
,
F.
, &
Gotman
,
J.
(
2009
).
High frequency oscillations and seizure frequency in patients with focal epilepsy
.
Epilepsy Research
,
85
(
2–3
),
287
292
. ,
[PubMed]

Competing Interests

Competing Interests: The authors have declared that no competing interests exist.

Author notes

Senior supervisor.

Handling Editor: Michael Breakspear

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. For a full description of the license, please visit https://creativecommons.org/licenses/by/4.0/legalcode.

Supplementary data