The study of neuron models and their networks is a riveting topic for many researchers worldwide because it allows to glimpse the fundamental processes using accessible methodology. The paper considers dynamics of small networks of Hodkin-Huxley neurons, namely a chain of three neurons and a small-world-like network of seven neurons. The ensembles of neurons are represented by systems of ordinary differential equations, so the research has been conducted numerically. It has been found that complex quasi-periodic and chaotic regimes may arise in the systems, and the existense of such regimes is caused by the inner parameters of the systems, such as individual currents of the neurons and the coupling between them. This research contributes to the fundamental understanding of signal propagation in networks of neuron models and may provide insight into the physiology of real neuronal systems.
Hodgkin-Huxley neuronsynchronizationneural oscillationsneural networksnetwork physiologyThe authors declare that financial support was received for the research and/or publication of this article. This work is supported by the Brain Program of the IDEAS Research Center.section-in-acceptanceNetworks of Dynamical SystemsIntroduction
Specialists from a wide range of scientific fields are interested in the processes occurring in the brain. Attempts to understand the human brain’s functioning increasingly incorporate concepts from physics, mathematics, computer science, mathematical biology, and related disciplines. So, the development of synergetic approaches (Haken, 1977) has provided a new perspective which can advance the investigation of such processes.
A significant number of fundamental interdisciplinary studies have been conducted by Hermann Haken, one of the founders of synergetics. A general description of the brain’s operational principles from a synergetics standpoint is detailed in (Haken, 2013) and the works cited therein, while there is a number of works that inquire about specific phenomena. For instance, some works cover the issue of spontaneous synchronization of neuronal spiking (Haken, 1977; Osipov et al., 2007), and the work (Haken, 1983) is devoted to the competition between oscillatory modes in neural networks and the emergence of a dominant mode.
Synergetics’ methodology and approaches have evolved significantly over the past decades, and modern science has come closer to understanding the functioning principles of the brain. However, the brain is a very complex object, and often we can only model a small part of it or are forced to be restricted with general small-world models (Watts and Strogatz, 1998).
Many research groups suggest phenomenological brain models of different complexity in order to understand its general behavior (Dimulescu et al., 2025; Cakan and Obermayer, 2020; Rosenblum, 2024). However, because these models are not directly grounded in the macroscopic characteristics of real neurons, key dynamical differences may exist between the numerical models and actual biological neurons.
Many studies exploring the dynamics of complex brain networks tend to utilize abstract neuron models, for example FitzHugh–Nagumo or integrate-and-fire models, as their partial elements (Rybalova et al., 2023; Wang et al., 2021; Anesiadis and Provata, 2022; Tyloo, 2024; Augustin and Obermayer, 2017). This preference is generally motivated by the relative simplicity of such systems: they often incorporate no more than three dynamic variables, are dimensionless, and are therefore easier to solve numerically.
In contrast, the Hodgkin–Huxley neuron (Hodgkin and Huxley, 1952) is a physiologically plausible model of spike generation, as it is founded on a macroscopic description of the neuronal membrane’s dynamics, which allows the researcher to draw parallels with parameters from real systems. This leads us to believe that the dynamics of networks constructed with Hodgkin–Huxley neurons are worthy of their own study.
Many works about the Hodgkin–Huxley networks consider random coupling structures (Majhi et al., 2025), which clearly have a high degree of similarity to the structures in a real brain. However, studying such structures poses a number of challenges, both technical and interpretative: such topology requires a large number of calculations and it is unclear which characteristics of the resulting system are suitable for drawing parallels with the real brain. In the presented paper we consider small, elementary structures consisting of only a few Hodgkin-Huxley neuron models with diffusive coupling (Sun et al., 2008). This topology is similar to the so-called ”small world network” concept with short distanse between the nodes, high clustering coefficient and connections through hub.
According to the principles of synergetics, the foundation of self-organization is the emergence of a new order and the increase in complexity of systems through random deviations in the states of their elements and subsystems. Such fluctuations are usually neutralized through negative feedback loops, which ensure the preservation of the structure and the system’s near-equilibrium state. However, in more complex open systems, due to the influx of energy from outside, deviations increase over time causing the effect of collective behavior of elements and subsystems. Ultimately, this process leads either to the destruction of the previous structure or to the emergence of a new order. In this work, we set the goal of understanding how simple neuron-like structures of Hodgkin-Huxley models behave, in order to further generalize the acquired knowledge to more complex nonequilibrium systems and move on to nonequilibrium ensembles composed of such elementary subnetworks. So, we dedicate this article to the memory of Hermann Haken, the founder of the Synergetics, who was actively involved in the study of the most complex processes occurring in the brain.
The paper is organised as follows. Section 2 introduces the model under consideration, explains its characteristics and describes the design of the numerical experiments. Section 3 presents the results of the experiments in full detail. Its Subsection 3.1 describes the results of the experiments for a chain of three coupled Hodgkin-Huxley neurons, while the Subsection 3.2 provides insight into the dynamics of a small network of seven coupled Hodgkin-Huxley neurons. Section 4 discusses the weaknesses and the prospects of the research, and Section 5 summarizes the findings.
Model and methods
This paper focuses on the regime formation and signal propagation in small ensembles of Hodgkin-Huxley neurons (Hodgkin and Huxley, 1952). The ensembles are defined by the following set of equations:dxidt=1CmḡKn4xi−xK+ḡNam3hxi−Na++ḡlxi−xl+Iexti+∑j=0,j≠iNwijxj−xi,dnidt=αnixi1−ni−βnixini.dmidt=αmixi1−mi−βmiximi,dhidt=αhixi1−hi−βhixihi.
Here, the first equation describes the dynamics of the neuron membrane potential x, dependent on the ion current flowing through ion channels in the membrane. The next three equations are responsible for the amount of open ion channels, regulating the ion currents. The α(x) and β(x) functions along with all the other parameters (Cm=1μF/cm2, xK=+12 mV, xNa=−115 mV, xl=−10.63 mV, ḡK=36 mmho/cm2, ḡNa=120 mmho/cm2 and ḡl=0.3 mmho/cm2) are taken directly from the original paper (Hodgkin and Huxley, 1952).
The sum ∑j≠iNwij(xj−xi) in (Equation 1) introduces linear electric coupling with wij determining coupling strength. The terms wij are portrayed as elements of coupling matrix W, which allows one to determine any kind of network topology with ease.
It is well-known that the value of the external current density Iext affects the change of the regime in a solitary Hodgkin-Huxley system (Hodgkin and Huxley, 1952; Rinzel and Miller, 1980), and the neuron can switch to one of three regimes. In the original Hodgkin-Huxley system, at values of Iext≤0μA/cm2 oscillations are hindered and the neuron is silent; at 0<Iext<8μA/cm2 the system shows the excitatory regime and there is a stable focus on the projection of the phase portrait in the x(n) plane (Figure 1a). But at Iext=8μA/cm2, a supercritical Andronov-Hopf bifurcation occurs and the neuron demonstrates self-oscillations, so a stable limit cycle arises on the projection of the phase space in the x(n) plane (Figure 1b).
Dynamical regimes in a solitary Hodgkin-Huxley neuron, projections of the phase space on the (x,n) plane (left column) and x(t) time realisations (right column). Excitatory regime for Iext3 = 3 μA/cm2(a) and self-oscillating regime for Iext3 = 12 μA/cm2(b).
Two panels labeled (a) and (b). Panel (a) shows two purple graphs: Left is a plot of n versus x in millivolts, depicting a loop; right is x versus time in milliseconds, showing a sharp decline and stabilization. Panel (b) displays two green graphs: Left is a plot of n versus x, forming a closed loop; right is x versus time, showing periodic spikes.
In this paper we aim to establish the influence of external current densities Iext and values of coupling strength w on the formation of various dynamical regimes and their propagation within a chain of three and a small ensemble of seven Hodgkin-Huxley neurons. Thereby, the initial conditions had to be fixed at certain values described in the following sections in order to exclude their influence on the dynamics of the considered ensembles. For an ensemble of two Hodgkin-Huxley neurons it has been shown that the regime can depend on the x0 initial conditions, while the other initial conditions n0, m0 and h0 have no influence on the regime formation in the system (Bogatenko et al., 2025).
In order to estimate the synchrony in a chain of three elements, Pearson correlation coefficient (Equation 2) (Pearson, 1896; Dunn and Clark, 1986; Rodgers and Nicewander, 1988) has been used:ρ=∑x1i−x1̄x2i−x2̄∑x1i−x1̄2∑x2i−x2̄2where x1i and x2i are the elements of the compared series x1 and x2, and x1̄ and x2̄ are the mean values of the series x1 and x2. The Pearson correlation coefficient takes values in the range [−1;1], where 1 means perfect positive or in-phase correlation, −1 means perfect out-of-phase correlation, and 0 means no correlation between the variables. The computation of the correlation coefficient was carried out for complete time realizations. Transient time periods were included in the calculations because they may contain critical information about the regime, for instance, single spikes in excitatory regime.
The research is conducted numerically utilising Runge-Kutta 4th order method (Runge, 1895; Kutta, 1901) over a time interval T=40000 with a step of 0.01 for each time realisation. A set of programs in C was used to carry out numerical integration, and the graphs were plotted with Gnuplot.
ResultsDynamical effects in a chain of three coupled Hodgkin-Huxley neurons
In this section a chain of three coupled Hodgkin-Huxley neurons of the topology shown in Figure 2 is under consideration. Here, the first neuron x1 receives Iext1=12μA/cm2 and self-oscillates, the second one x2 shows excitatory mode for Iext2=3μA/cm2, and Iext3 for the third neuron x3 is varied. Also, the second neuron x2 always influences the neurons x1 and x3 with a deliberately low value of coupling strength wout=0.1, while w21=w23=w are varied. Thus, we analyze the influence of external current density Iext3 and coupling strength w on the dynamical regimes and their synchrony in this chain. Iext3 were changed in the interval [-15; 15] μA/cm2, and w took values in [0; 4]. Initial conditions are the same for all the neurons: x01=x02=x03=10, n01=n02=n03=0.1, m01=m02=m03=0.01, h01=h02=h03=0.01, which allows us to exclude their influence on the system’s dynamics.
Schematic image of an ensemble of three Hodgkin-Huxley neurons under consideration.
Diagram showing three interconnected nodes labeled \(x_1\), \(x_2\), and \(x_3\). Arrows between nodes are labeled with weights \(w_{21}\), \(w_{23}\), and \(w_{out}\), with \(w_{out} = 0.1\). Graphs within nodes depict signal patterns.
Figure 3 shows the Pearson correlation coefficient maps for each of the three pairs of neurons in the chain, while Figures 4–6 show some typical time realizations of x(t) and projections of the phase portraits onto the (x,n) plane. In general, one can note that the absolute value of correlation coefficient ρ is equal or close to one for x1 and x2 on the most of the considered parameter plane, while it is generally closer to 0 for the other pairs. It means that the neurons x1 and x2 often show correlated behaviour, however, the negative value of the coefficient indicates antiphase oscillations. In all, the presence of areas with negative correlation coefficient values on the three maps is caused by a delay occurring in the neurons (Figure 4).
Pearson correlation coefficient maps on the plane (Iext3,w) of the three pairs of neurons in a chain: x1 and x2(a), x2 and x3(b), x1 and x3(c). Other parameters and initial conditions: Iext1=12μA/cm2, Iext2=3μA/cm2, x01=x02=x03=10, n01=n02=n03=0.1, m01=m02=m03=0.01, h01=h02=h03=0.01.
Three side-by-side heatmaps labeled (a), (b), and (c) illustrate data variations. The x-axis ranges from -15 to 15, and the y-axis labeled as \( w \) from 0 to 4. The color gradient from yellow to purple indicates different values of \( \rho \) from 1 to -1. Distinct vertical bands are visible, with subtle changes in color across the maps.
An example of delay in a chain of three coupled Hodgkin-Huxley neurons for Iext3=1.0μA/cm2 and w=0.2: x(t) time realisations and projections of the phase space on the (x,n) plane of the neurons x1 (blue), x2 (red) and x3 (green). Other parameters and initial conditions: Iext1=12μA/cm2, Iext2=3μA/cm2, x01=x02=x03=10, n01=n02=n03=0.1, m01=m02=m03=0.01, h01=h02=h03=0.01.
Four panels display different graph analyses. The first panel shows voltage over time with three color-coded curves in blue, red, and green. The next three panels illustrate phase space plots with each plot focusing on a different variable: blue for \(x_1\), red for \(x_2\), and green for \(x_3\), each demonstrating loop patterns against voltage.
Besides, on each of the maps in Figure 3a vertical region of complete positive correlation (ρ=1) is noticeable for Iext3 = 12 μA/cm2, in which the neurons are completely synchronous. The presence of this region precisely at the value Iext3 = 12 μA/cm2 is due to the fact that the first neuron receives an external current of exactly this magnitude: Iext1 = 12 μA/cm2, and the self-oscillations arising in it suppress the excitable mode of the second neuron (Iext2 = 3 μA/cm2). Note that the suppression of the excitable mode occurs smoothly in terms of the coupling strength w: in the absence of coupling, as well as at small coupling strength values w≪1, the correlation coefficient is close to zero and the neurons are not synchronous. Then, with an increase in the coupling strength, the correlation coefficient begins to increase and at w≈1.5 it becomes equal to 1 (Figure 3). This effect has been shown earlier for two coupled Hodgkin-Huxley neurons (Bogatenko et al., 2025). On the correlation map of neurons x1 and x3 (Figure 3c), the coefficient ρ equals to 1 for all the values of coupling strength in the considered range at Iext3 = 12 μA/cm2, since in this case these two neurons show the same regime throughout the experiment.
What the correlation maps in Figure 3 do not show is boundaries between oscillation modes. However, as shown in some representative examples, the neurons tend to exhibit oscillations of varying complexity throughout the entire parameter plane under consideration. Self-oscillations of classical spikes can be synchronized with high accuracy (ρ=1) (Figure 5b), but there also are oscillations that resemble spikes in shape but have a much smaller amplitude of up to 5 mV (Figure 5a). Also, quasi-periodic oscillation regimes of varying degrees of complexity can be realized in the system (Figure 6).
An example of self-oscillating regimes in a chain of three coupled Hodgkin-Huxley neurons for Iext3=−14.0μA/cm2, w=0.2(a) and Iext3=12.0μA/cm2, w=0.7(b). x(t) time realisations and projections of the phase space on the (x,n) plane of the neurons x1 (blue), x2 (red) and x3 (green). Other parameters and initial conditions: Iext1=12μA/cm2, Iext2=3μA/cm2, x01=x02=x03=10, n01=n02=n03=0.1, m01=m02=m03=0.01, h01=h02=h03=0.01.
Two panels, (a) and (b), each containing four graphs. Panel (a) shows a plot of voltage \(x, \text{mV}\) versus time \(t, \text{msec}\) with three curves and three phase space plots \(x, \text{mV}\) vs \(n\) for \(x_1\), \(x_2\), and \(x_3\). Panel (b) has a similar structure, with the voltage versus time showing three overlapping curves and repeated phase space plots for \(x_1\), \(x_2\), and \(x_3\), colored blue, red, and green, respectively.
An example of quasiperiodic regimes in a chain of three coupled Hodgkin-Huxley neurons for Iext3=1.0μA/cm2, w=0.4(a) and Iext3=1.0μA/cm2, w=0.9(b). x(t) time realisations and projections of the phase space on the (x,n) plane of the neurons x1 (blue), x2 (red) and x3 (green). Other parameters and initial conditions: Iext1=12μA/cm2, Iext2=3μA/cm2, x01=x02=x03=10, n01=n02=n03=0.1, m01=m02=m03=0.01, h01=h02=h03=0.01.
Two panels labeled (a) and (b) feature time series and phase portraits of voltage dynamics. Left graphics show voltage vs. time with overlapping curves in blue, red, and green. Right graphics display phase portraits for three variables \( X_1 \), \( X_2 \), and \( X_3 \) with corresponding colors. Panel (a) presents smoother dynamics compared to (b), highlighting differing system responses.Signal propagation in a network of seven coupled Hodgkin-Huxley neurons
Now let us construct a larger network of Hodgkin-Huxley neurons, consisting of two smaller chains connected unidirectionally via a hub (Figure 7). Within the chains, the neurons are connected bidirectionally, and the connection values between them remain fixed in all numerical experiments.
Schematic image of an ensemble of seven Hodgkin-Huxley neurons under consideration.
Diagram showing a network of nodes connected by arrows labeled with the letter "w". Central node labeled "x_hub" connects to nodes "x2" and "x5". Nodes "x1", "x3", "x4", and "x6" surround "x2" and "x5", each displaying a unique waveform.
At the beginning of each experiment, in the first chain one of the outer neurons is in the self-oscillation mode (Iext1 = 12 μA/cm2), and the other two show the excitatory mode, but with a small detuning in the value of the external current density: Iext2 = 3 μA/cm2, Iext3 = 2 μA/cm2. Also, the outer neurons intentionally exert a weak influence on the central neuron (w21=w23=0.1), while in the opposite direction, the coupling strength is deliberately set high – w12=w32=2.5. Such a combination of external currents and coupling strength values forces a certain regime to establish in this layer: all the three neurons are nearly synchronous in the excitatory mode (Figures 8a–c).
Structures in the first (green) and the second (blue) layer of the considered network of seven Hodgkin-Huxley neurons in each experiment: x(t) time realisations of x1(a), x2(b), x3(c), x4(d), x5(e), and x6(f). Parameters: Iext1 = Iext4 = 12 μA/cm2, Iext2 = Iext5 = 3 μA/cm2, Iext3 = 2 μA/cm2, Iext6 = 7 μA/cm2; w21=w23=0.1, w12=w32=2.5; w54=w56=0.1, w45=w65=0.3.
Six graphs display voltage over time, labeled (a) to (f). The y-axis represents voltage in millivolts, and the x-axis represents time in milliseconds. Graphs (a), (b), and (c) show initial voltage spikes stabilizing to a steady state, while graphs (d), (e), and (f) display repetitive voltage peaks.
In the second chain, there is a similar balance of connection strengths as in the first one, but here the central neuron acts on the outer neurons with a weaker coupling strength compared with the first layer: w54=w56=0.1, w45=w65=0.3. Here, similar regimes are realized in the neurons: Iext4 = 12 μA/cm2, Iext5 = 3 μA/cm2, Iext6 = 7 μA/cm2, however, we note that the current density Iext6 received by neuron x6 is of a prethreshold value for the Andronov-Hopf bifurcation. Now, self-oscillating regime establishes in this chain at the beginning of every experiment (Figures 8d–f).
Here we aspire to establish the influence of the coupling strength between the layers through the hub and the hub’s regime on the dynamics that is translated to the second chain. This goal motivates the choice of constant values of the coupling strength and the current density values within the chains. For the initial conditions in the system in all the experiments described in this section, a single set of values with a uniform distribution over the intervals [9.5;10.5] was chosen for the set of seven variables x, [0.095;0.105] for the set of seven variables n, and [0.0095;0.0105] for the sets of variables m and h (seven values each). The initial conditions were chosen in this way in order to introduce a degree of dissimilarity and bring the numerical experiment closer to real systems.
Figure 9 shows the regime map for the hub neuron on the parameter plane (Ihub, w). One can see that the hub neuron can show one of three regimes: excitatory, single burst, and period-1 self-oscillations. Self-oscillations are observed in a limited parameter range with low coupling strength w<0.5 and current density values greater than 7 μA/cm2 (Figure 10a). As the coupling strength increases, the region of self-oscillations decreases—here one can see how the regime established in the first chain suppresses the self-oscillation regime in the hub neuron. The excitatory regime predominates over most of the plane under consideration (Figure 10c). The boundary between these regimes is expressed in the regime of single burst generation (Figure 10b).
Regime map of the hub neuron xhub of the considered network of seven Hodgkin-Huxley neurons on the parameter plane (Ihub,w). Regime scale: EXC–excitatory regime, BUR–single burst with subsequent silence, P1 – periond-1 oscillations. Parameters are stated in Section 3.2.
Heatmap showing regime transitions as a function of \( I_{\text{hub}} \) in microamperes per square centimeter and \( w \). The map has predominant red, indicating the BUR regime, with a small yellow section at higher \( I_{\text{hub}} \) values indicating the P1 regime. The y-axis represents \( w \) from 0 to 3.
An example of typical regimes in the hub neuron xhub of the considered network of seven Hodgkin-Huxley neurons: x(t) time realisations and projections of the phase space on the (x,n) plane. Period-1 oscillations (a) (Ihub = 14 μA/cm2, w=0.2), single burst (b) (Ihub = 14 μA/cm2, w=0.6), and excitatory regime (c) (Ihub = 14 μA/cm2, w=1.4) are shown. Other parameters are stated in Section 3.2.
Three panels labeled (a), (b), and (c) with pairs of graphs in each. Panel (a) shows voltage over time and a trajectory on a phase plane. Panel (b) shows damped oscillations on a time series and a loop-shaped trajectory. Panel (c) shows rapid stabilization over time and a spiral-shaped phase plane trajectory. Each graph is distinctly labeled with axes in millivolts, milliseconds, and a variable eta.
Regime maps for the neurons in the second layer on the same parameter plane are shown in Figure 11. Within the considered parameter ranges, neurons can exhibit a variety of simpler and more complex regimes: excitatory regime (EXC), single burst with subsequent silence (BUR), single burst of low-amplitude oscillations (up to 5 mV) (LA BUR), period-1 self-oscillations (P1), low-amplitude self-oscillations (up to 5 mV) (LA P1), period-2 self-oscillations (P2), quasi-periodic (QUA), and chaotic (CH) oscillations.
Regime map of the neurons x5(a) and x6(b) of the second layer of the considered network of seven Hodgkin-Huxley neurons on the parameter plane (Ihub,w). Regime scale: BUR–single burst with subsequent silence, LA BUR–single burst of low-amplitude oscillations (up to 5 mV), P1 – periond-1 oscillations, LA P1 – low-amplitude self-oscillations (up to 5 mV), P2 – period-2 self-oscillations, QUA–quasi-periodic oscillations, CH–chaotic oscillations. Parameters are stated in Section 3.2.
Two side-by-side contour plots labeled (a) and (b) show regime distributions over a grid with axes labeled \(I_{hub}, \, \mu A/cm^2\) and \(w\). Various colored regions represent different regimes: orange, blue, green, purple, red, and others, with a legend on the right mapping colors to regime names, including CH, QUA, P2, LA P1, P1, LA BUR, and BUR.
The regime maps for neurons of the second chain x5 and x6 have a more complex structure than the map for the hub neuron. Here, for small values of the coupling strength with the hub neuron w<0.3, the self-oscillation mode of period 1 prevails in the system - for these parameter values neuron x4 transmits its self-oscillation mode to neurons x5 and x6 (Figure 12). In the range of values w∈[0.3,0.5], there is a region of transient modes, when the value of the coupling strength with the hub neuron w becomes greater than the value of the coupling strength within the chain w45=w65=0.3. Here we can observe various complex regimes: period-2 oscillations, quasi-periodic and chaotic oscillations (Figure 13).
An example of a self-oscillating regime in the second layer of the considered network of seven Hodgkin-Huxley neurons for Ihub = 6 μA/cm2 and w=0.1: x(t) time realisations and projections of the phase space on the (x,n) plane of the neurons x4 (blue), x5 (red) and x6 (green). Other parameters are stated in Section 3.2.
Graphs displaying voltage data over time and phase diagrams. The first two graphs show oscillating voltage patterns over time in milliseconds with voltages ranging from negative one hundred to twenty millivolts. The last three graphs depict phase diagrams with x and n axes, labeled \( x_4 \), \( x_5 \), and \( x_6 \), exhibiting looped trajectories in blue, red, and green, respectively. The phase diagrams display dynamic voltage relationships.
Examples of complex transient regimes in the second layer of the considered network of seven Hodgkin-Huxley neurons: x(t) time realisations and projections of the phase space on the (x,n) plane of the neurons x4 (blue), x5 (red) and x6 (green). Cases for Ihub = 4 μA/cm2, w=0.4(a), Ihub = 14 μA/cm2, w=0.5(b), and Ihub = 7 μA/cm2, w=0.5(c) are shown. Other parameters are stated in Section 3.2.
Graphs showing time series and phase plane plots of dynamical systems in three panels labeled (a), (b), and (c). Each panel has two time series plots on the left with oscillations in blue, red, and green. The right side has three phase plane plots with loops in blue, red, and green, marked as \(x4\), \(x5\), and \(x6\). Axes are labeled with time in milliseconds and voltage in millivolts or a variable \(n\).
With further increase in the coupling strength with the hub neuron w, the regime map becomes more regular. Here, one can see the predominance of two modes: period-1 self-oscillations of low amplitude (Figure 14) and single burst with subsequent silence (Figure 14). Low-amplitude period-1 self-oscillations occur with an amplitude smaller than in the classical self-oscillatory mode of the system and are approximately 5 mV. Interestingly, the boundary between these two modes has the form of a hyperbolic or exponential function and is expressed by the mode of generating a low-amplitude burst (Figure 14). It is also interesting to note that for neuron x6, this mode of generating a small-amplitude spike train is not only transient, as for neuron x5, but also a stable mode, which is realized at negative values of the current density Ihub (Figure 9).
Examples of stable regimes in the second layer of the considered network of seven Hodgkin-Huxley neurons: x(t) time realisations and projections of the phase space on the (x,n) plane of the neurons x4 (blue), x5 (red) and x6 (green). Cases for Ihub = 7 μA/cm2, w=0.1(a), Ihub = 7 μA/cm2, w=1.1(b), and Ihub = 7 μA/cm2, w=2.3(c) are shown. Other parameters are stated in Section 3.2.
Graphs depicting varying neuronal action potentials and phase plots across three sections labeled (a), (b), and (c). Each section presents time versus potential graphs and phase diagrams for variables \(x_4\), \(x_5\), and \(x_6\), in blue, red, and green, respectively. Panels illustrate differences in neuronal activity with distinct oscillatory and non-oscillatory behaviors over time intervals, showing diverse trajectories in the phase spaces. The \(x\) axis represents membrane potential in millivolts, and the \(t\) axis represents time in milliseconds.Discussion
The obtained results demonstrate that elementary Hodgkin-Huxley neuronal networks, renowned for their physiological fidelity, are capable of generating a rich repertoire of complex dynamics based solely on their internal structure and connectivity, without the need for complex external driving forces. Our study of two different topologies reveals several fundamental principles that govern the emergence of this complexity.
A key limitation of the current study is its deliberate focus on a set of specific, predetermined network topologies, which, while providing a controlled framework for analysis, may not fully capture the architectural diversity found in biological neural systems. Future research must therefore prioritize scaling these investigations to larger, more complex networks to determine if the observed dynamics–such as the interplay between sodium-channel persistence and delayed-rectifier potassium currents in shaping burst termination–are preserved or fundamentally altered. It will be crucial to examine how these dynamics are influenced by other overarching network characteristics. For instance, introducing a hierarchical organization, where clusters of neurons with specific firing properties are nested within larger functional modules, could reveal novel emergent computational states. Furthermore, the inclusion of adaptive, plastic synapses based on spike-timing-dependent plasticity (STDP) rules would transform the model from a static circuit into a dynamic, learning system. This would allow us to investigate whether the intrinsic neuronal properties governed by the Hodgkin-Huxley formalism interact synergistically or competitively with experience-dependent synaptic changes to guide network development and stability.
The choice of coupling in the model, presented in Equation 1, is also in line with maintaining the scalability of our research to larger and complex networks and is due to two main reasons. Firstly, its mathematical simplicity and well-defined properties provide a clear and interpretable framework for isolating the effects of network topology on the emergent dynamics of intrinsically bursting neurons. This allows us to distinguish effects stemming from the network architecture from those arising from the complexities of synaptic transmission. Secondly, for the purpose of a qualitative analysis of mode generation in elementary neuromorphic structures, diffusion coupling offers a computationally efficient yet biologically grounded foundation. While models incorporating detailed chemical synapses are essential for simulating specific cortical circuits, their complexity would obscure the primary focus of this study: the interplay between topology and intrinsic neuronal dynamics. Certainly, there are numerous ways to introduce coupling between Hodgkin-Huxley neurons [see (Rossoni et al., 2005; Ma and Tang, 2017; Petousakis et al., 2023; Catterall et al., 2012) and many others], including special software (Hines and Carnevale, 1997; Carnevale and Hines, 2006; Bologna et al., 2022; Willms, 2002), all of which are biologically justified to varying degrees. More complex models of coupled Hodgkin-Huxley neurons are typical primarily for modeling specialized cortical regions and parts of the nervous system and, therefore, are beyond the scope of the problem being solved here.
In assessing the biological plausibility of the model, it is important to consider the key output parameters it generates. Characteristic current values range from a few to tens of nanoamperes and fall well within the recognized physiological range for mammalian neurons (Kole and Stuart, 2012; Stuart and Sakmann, 1994). Moreover, the temporal dynamics with spike durations of 1–2 ms accurately reflect those observed in experimental recordings from various cortical and hippocampal cell types (McCormick et al., 1985). The firing rates obtained in our simulations, often in the gamma range (20–80 Hz) under moderate inputs, are possibly related to oscillations associated with cognitive processing (Buzsaki and Draguhn, 2004). Although the model may exhibit higher rates under extreme forcing forces–a known property of the underlying Hodgkin–Huxley formalism, this does not detract from its usefulness for studying the fundamental network dynamics of spikes initiation, propagation, and synchronization, which was the primary goal of this work. Therefore, we are confident that the model provides a physiologically sound platform for investigating the interaction between the intrinsic properties of neurons and network topology.
Conclusion
In a small chain of coupled Hodgkin-Huxley neurons, we observe that the classical self-oscillatory mode, although dominant, is not exclusive. The emergence of more complex quasi-periodic oscillations and the observed delay in signal propagation along the network underscore a key conclusion: the intrinsic properties of the network itself are the primary source of complexity. Specifically, it is the precise interplay of three factors that regulates these dynamic modes: network topology, neurons’ coupling strength, and the combination of individual neuron current densities. It is important because it suggests that even relatively simple, localized neural circuits harbor a latent capacity to form complex temporal patterns, which may be fundamental to processes such as central pattern generation or sensory processing.
This principle is further reinforced in the small-world ensemble with a hub. Here, the hub neuron acts as a nonlinear processor and integrator, enabling the system to generate not only quasi-periodic, but also chaotic oscillations. Parameter induced transitions between dynamical regimes (silence, periodic, quasi-periodic, and chaotic) have well-defined boundaries in the regime maps for the parameters Ihub and w. But this predictability is non-trivial–it indicates that, despite the potential for chaos, the system’s behavior is determined by an underlying order that can be controlled by specific biological analogs–the excitability of the key neuron and the coupling strength.
More broadly, our results are consistent with the synergetic view of the brain as a self-organizing system. The spontaneous emergence and predictable transitions between complex oscillatory modes, driven by intrinsic network parameters rather than external stimuli, provide a plausible model for the dynamic switching of neural ensembles between functional states. The hub-based model, in particular, offers a mechanistic explanation for how local changes (e.g., neuromodulation affecting Ihub and w) can lead to global shifts in network dynamics, which may be important for understanding the role of highly connected nodes in neural networks.
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Author contributions
TB: Formal Analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft. KS: Conceptualization, Methodology, Software, Validation, Writing – original draft. GS: Conceptualization, Methodology, Supervision, Writing – review and editing.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be interpreted as a potential conflict of interest.
The handling editor ES declared a past co-authorship with the author GS.
The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.
Generative AI statement
The authors declare that no Generative AI was used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Edited by:Eckehard Schöll, Technical University of Berlin, Germany
Reviewed by:Semen Kurkin, Immanuel Kant Baltic Federal University, Russia
Anatoly Karavaev, Saratov Branch of the Institute of RadioEngineering and Electronics of Russian Academy of Sciences, Russia
ReferencesAnesiadisK.ProvataA. (2022). Synchronization in multiplex leaky integrate-and-fire networks with nonlocal interactions. Front. Netw. Physiology2, 910862. 10.3389/fnetp.2022.91086236926067AugustinL. J. B. F. M.ObermayerK.BaumannF. (2017). Low-dimensional spike rate models derived from networks of adaptive integrate-and-fire neurons: comparison and implementation. PLoS Comput. Biol.13, e1005545. 10.1371/journal.pcbi.100554528644841BogatenkoT.SergeevK.StrelkovaG. (2025). The role of coupling and external current in two coupled hodgkin–huxley neurons. Chaos An Interdiscip. J. Nonlinear Sci.35, 023149. 10.1063/5.024343339977306BolognaL.SmirigliaR.LupascuC.AppukuttanS.DavisonA.IvaskaG. (2022). The ebrains hodgkin-huxley neuron builder: an online resource for building data-driven neuron models. Front. Neuroinformatics16, 991609. 10.3389/fninf.2022.99160936225653BuzsakiG.DraguhnA. (2004). Neuronal oscillations in cortical networks. Science304, 1926–1929. 10.1126/science.109974515218136CakanC.ObermayerK. (2020). Biophysically grounded mean-field models of neural populations under electrical stimulation. PLoS Comput. Biol.16, e1007822. 10.1371/journal.pcbi.100782232324734CarnevaleN.HinesM. (2006). The NEURON book. Cambridge University Press.CatterallW.RamanI.RobinsonH.SejnowskiT.PaulsenO. (2012). The hodgkin-huxley heritage: from channels to circuits. J. Neurosci.32, 14064–14073. 10.1523/JNEUROSCI.3403-12.201223055474DimulescuC.StrömsdörferR.FlöelA.ObermayerK. (2025). On the robustness of the emergent spatiotemporal dynamics in biophysically realistic and phenomenological whole-brain models at multiple network resolutions. Front. Netw. Physiology5, 1589566. 10.3389/fnetp.2025.158956640861379DunnO. J.ClarkV. A. (1986). Applied statistics: analysis of variance and regression. John Wiley and Sons, Inc.HakenH. (1977). Brain dynamics: synchronization and activity patterns in pulse-coupled neural nets with delays and noise. Berlin: Springer.HakenH. (1983). Synergetics: an introduction, vol. 1 of springer series in synergetics. 3rd rev. and enl. Berlin: Springer.HakenH. (2013). Principles of brain functioning: a synergetic approach to brain activity. Berlin: Springer.HinesM.CarnevaleN. (1997). The neuron simulation environment. Neural. Comput.9, 1179–1209. 10.1162/neco.1997.9.6.11799248061HodgkinA.HuxleyA. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol.117, 500–544. 10.1113/jphysiol.1952.sp00476412991237KoleM.StuartG. (2012). Signal processing in the axon initial segment. Neuron73, 235–247. 10.1016/j.neuron.2012.01.00722284179KuttaW. (1901). Beitrag zur näherungsweisen integration totaler differentialgleichungen. Z. für Math. Phys.46, 435–453.MaJ.TangJ. (2017). A review for dynamics in neuron and neuronal network. Nonlinear Dyn.89, 1569–1578. 10.1007/s11071-017-3565-3MajhiS.GhoshS.PalP. K.PalS.PalT. K.GhoshD. (2025). Patterns of neuronal synchrony in higher-order networks. Phys. Life Rev.52, 144–170. 10.1016/j.plrev.2024.12.01339753012McCormickD.ConnorsB.LighthallJ.PrinceD. (1985). Comparative electrophysiology of pyramidal and sparsely spiny stellate neurons of the neocortex. J. Neurophysiology54, 782–806. 10.1152/jn.1985.54.4.7822999347OsipovG.KurthsJ.ZhouC. (2007). Synchronization in oscillatory networks. Berlin: Springer.PearsonK. (1896). Vii. mathematical contributions to the theory of evolution.—iii. regression, heredity, and panmixia. Philos. Trans. R. Soc. Lond Ser. A187, 253–318.PetousakisK.-E.ApostolopoulouA.PoiraziP. (2023). The impact of hodgkin–huxley models on dendritic research. J. Physiol.601, 3091–3102. 10.1113/JP28275636218068RinzelJ.MillerR. (1980). Numerical calculation of stable and unstable periodic solutions to the Hodgkin-Huxley equations. Math. Biosci.49, 27–59. 10.1016/0025-5564(80)90109-1RodgersJ. L.NicewanderW. A. (1988). Thirteen ways to look at the correlation coefficient. Amer Stat.42, 59–66. 10.2307/2685263RosenblumM. (2024). Feedback control of collective dynamics in an oscillator population with time-dependent connectivity. Front. Netw. Physiology4, 1358146. 10.3389/fnetp.2024.135814638371453RossoniE.ChenY.DingM.FengJ. (2005). Stability of synchronous oscillations in a system of hodgkin-huxley neurons with delayed diffusive and pulsed coupling. Phys. Rev. E71, 061904. 10.1103/PhysRevE.71.06190416089762RungeC. D. T. (1895). Über die numerische auflösung von differentialgleichungen. Math. Ann.46, 167–178. 10.1007/BF01446807RybalovaE.BogatenkoT.BukhA.VadivasovaT. (2023). The role of coupling, noise and harmonic impact in oscillatory activity of an excitable fitzhugh-nagumo oscillator network. Izv. Saratov Univ. Phys.4, 294–306. 10.18500/1817-3020-2023-23-4-294-306StuartG.SakmannB. (1994). Active propagation of somatic action potentials into neocortical pyramidal cell dendrites. Nature367, 69–72. 10.1038/367069a08107777SunX.PercM.LuQ.KurthsJ. (2008). Spatial coherence resonance on diffusive and small-world networks of hodgkin-huxley neurons. Chaos18, 023102. 10.1063/1.290040218601469TylooM. (2024). Resilience of the slow component in timescale-separated synchronized oscillators. Front. Netw. Physiology2, 1399352. 10.3389/fnetp.2024.139935238962160WangZ.XuY.LiY.KapitaniakT.KurthsJ. (2021). Chimera states in coupled hindmarsh-rose neurons with stable noise. Chaos, Solit. Fractals148, 110976. 10.1016/j.chaos.2021.110976WattsD.StrogatzS. (1998). Collective dynamics of ’small-world’ networks. Nature393, 440–442. 10.1038/309189623998WillmsA. (2002). Neurofit: software for fitting hodgkin-huxley models to voltage-clamp data. J. Neurosci. Methods121, 139–150. 10.1016/s0165-0270(02)00227-312468004