Weak coupling of neurons enables very high-frequency and ultra-fast oscillations through the interplay of synchronized phase shifts

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 f–I 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.

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 I_ext, 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.

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 I_ext (e.g., see Liu, Liu, & Liu, 2014; Tsumoto et al., 2006; Xing et al., 2022) or g_Ca, g_K, and g_L that can also be influenced by drugs (see Macdonald, 1989; Nicholson, Blanche, Mansfield, & Tran, 2002).

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.

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 C₁ and C₂, and ε and C₁, respectively, are shown. More accurately, Figure 3 precisely corresponds to the dashed sections at C₂ = 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 C₁ and C₂ 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 C₁ and C₂ 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.

Similarly, the bifurcation diagram in Figure 4B shows the Arnold tongues with respect to the coupling strength ε and membrane capacitance C₁. In this case, they emanate from the point C₁ = C₂ (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 C₁, C₂, 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 = V₁ + V₂ 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.

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 C_i, as is depicted in Figures S1 and S2 in the Supporting Information.

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.

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 I_ext = 43 + p_i(t) + σ dW with σ = 2.0 and p_i = 40 if t ∈ [716, 718] ms in Equation 15 for i ∈ 𝒫₁. This stimulus is applied only to half of the population, that is, p_i = 0 for i ∈ 𝒫₂. 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.

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 C₁, C₂, 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 C₁ × C₂ and C₁ × ε (f–I curves for a single interneuron for various values of C are provided in the Supporting Information (see Figure S4); the default parameter choice of I_ext = 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 = V₁ + V₂ 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 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).

As one can see in Figure 10, this configuration leads to the emergence of transient UFRs in the composed signal V = $∑i=14$V_i. 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).

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

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 (f–I curves for a single neuron for various values of C are provided in the Supporting Information (see Figure S5); the default parameter choice of I_ext = 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).

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 g_Na and reversal potential V_Na 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.

The system parameters g_L, g_Ca, and g_K represent the maximum leak, Ca²⁺, and K⁺ electrical conductances through membrane channels, respectively, and V_L, V_Ca, and V_K are reversal potentials of the specific ion channels. Finally, I_ext 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/cm² (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.

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.

In this article, we consider parameter values C ≈ 1 μS/cm², g_L = 0.1 mS/cm², g_Na = 30 mS/cm², g_K = 20 mS/cm², V_L = −60 mV, V_Na = 45 mV, and V_K = −80 mV. These values are within physiological ranges and give the high-frequency firing rates typical of hippocampal interneurons (White et al., 1998).

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).

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/cm², g_L = 0.019 mS/cm², g_Na = 120 mS/cm², g_Kdr = 100 mS/cm², g_M = 2 mS/cm², V_L = −65 mV, V_Na = 55 mV, and V_K = − 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.