Hamiltonian energy and coexistence of hidden firing patterns from bidirectional coupling between two different neurons.

In this paper, bidirectional-coupled neurons through an asymmetric electrical synapse are investigated. These coupled neurons involve 2D Hindmarsh-Rose (HR) and 2D FitzHugh-Nagumo (FN) neurons. The equilibria of the coupled neurons model are investigated, and their stabilities have revealed that, for some values of the electrical synaptic weight, the model under consideration can display either self-excited or hidden firing patterns. In addition, the hidden coexistence of chaotic bursting with periodic spiking, chaotic spiking with period spiking, chaotic bursting with a resting pattern, and the coexistence of chaotic spiking with a resting pattern are also found for some sets of electrical synaptic coupling. For all the investigated phenomena, the Hamiltonian energy of the model is computed. It enables the estimation of the amount of energy released during the transition between the various electrical activities. Pspice simulations are carried out based on the analog circuit of the coupled neurons to support our numerical results. Finally, an STM32F407ZE microcontroller development board is exploited for the digital implementation of the proposed coupled neurons model.

Introduction

The human brain is an organ that can exhibit extremely complex nonlinear behavior (Natarajan et al. 2004). It is constituted of a huge number of interconnected neurons. Neurons can be viewed as the elementary, structural, and functional elements of the central nervous system (Tsou et al. 1998). Neurons play an important role in the processes of recording, selection, storage, and data transfer in the brain’s activities (Leuthardt et al. 2004). To ease the understanding of the dynamics of electrical activities in the brain, several mathematical models of neurons have been proposed and investigated in the literature, among which the Hodgkin–Huxley neuron (Hodgkin and Huxley 1952), Izhikevich neuron (Izhikevich 2003a), Morris–Lecar neuron (Tsumoto et al. 2006), 2-D and 3-D Hindmarsh–Rose (HR) model (Hindmarsh and Rose 1982, 1984), FitzHugh–Nagumo (FN) model (Izhikevich and FitzHugh 2006), Hopfield neural network model (Njitacke et al. 2020a, 2020b; Doubla Isaac et al. 2020; Tabekoueng Njitacke et al. 2020) and Chay model (Chay 1985) just to name few. These various neurons have been intensively studied separately as well as coupled. Concerning the investigation of the single neuron, (Zhang et al. 2021) introduced a novel free-equilibrium HR model using memristive electromagnetic induction. During their investigations, the authors found that the proposed model of the neuron is capable of unusual and rare phenomena of hidden homogeneous extreme multistability, defined by the coexistence of an infinite number of firing activities with the same shape but at different positions. It is good to mention that such types of nonlinear behaviors involving multistability, hidden dynamics previously found in neurons were also able to occur in circuits and systems, such as the extended Lu system(Lai et al. 2018a), a unified chaotic system (Lai 2021), a two-memristor-based chaotic system (Lai et al. 2020a) as well as a non-equilibrium chaotic system (Lai et al. 2020b). In (Ngouonkadi et al. 2016), the authors studied in an extended model of a 4D HR neuron model the phenomenon of coexistence of firing patterns. The investigations revealed that the model is able to exhibit several types of complex phenomena, including symmetry breaking, period doubling, reverse period-doubling, and crisis scenarios. Hou et al. (Hou et al. 2021) have estimated the electrical activity in a neuron under a depolarization field. During their investigation, they found that there was a peak in the firing interval of neurons with the increase of stimulation current intensity. In the presence of an electric field, the firing pattern of the neuron was transformed from single bursting to intermittent multimodal bursting with the increase in frequency and amplitude of the electric field. Kafraj et al. (Kafraj et al. 2020) have introduced and investigated a three-variable memristive Izhikevich model to describe the behavior of neurons under electromagnetic induction and noise. Their model represented the effect of internal and external magnetic fields on neurons. The improved model without the external magnetic field was able to exhibit firing patterns such as regular spiking, resonator, chattering, fast-spiking, chaotic spiking, and chaotic bursts. The presence of the external magnetic field was able to change the firing pattern of the neuron, for example, shifting from chaotic firings to periodic ones or vice versa. It was also observed that the external field effect stimulated the neuron to fire, while in the absence of the external field, it was at rest.

Besides the above-mentioned advances, investigation of the coupled neurons has also attracted a lot of attention. For example, Pisarchik et al. (Pisarchik et al. 2018) explored the dynamics of two identical 3D HR coupled models through an asymmetric electrical coupling. The result indicated that the 6D neuronal oscillator obtained was able to exhibit the phenomenon of multistability with up to six firing activities. In a real brain, it is very difficult to have a connection between two identical neurons with identical parameters. Henceforth, Njitacke et al. (Tabekoueng Njitacke et al. 2020) have considered the dynamics of two coupled 2D HR neurons with a small discrepancy of their parameters. From their investigations, the authors found that their proposed simple 4D neural oscillator was able to exhibit traditional brain firing activities including bursting, spiking oscillations as well as the phenomenon of bistability under the variation of the symmetric coupling weight and the external stimulus. In (Wang et al. 2016) the authors coupled the neurons through various chemical synapses, namely: rapid excitatory and inhibitory synapses. The spatiotemporal patterns were analyzed in the coupled model. The results indicated that the rapid excitatory synapse coupling type is easier to produce periodic spatiotemporal patterns than the rapid inhibitory synapse coupling type, and the process was analyzed using the bifurcation of a unique neuron model. They also introduced the permutation entropy, defined as a measure of network firing complexity to explore the process of formation and transition of spatiotemporal patterns. Karthikeyan et al. (Rajagopal et al. 2020) have introduced a modified Hindmarsh-Rose neuron model having a fractional-order threshold magnetic flux. Their investigations were conducted both in the presence and absence of external electromagnetic induction. Besides, the emergence of the spiral waves in the network of the proposed model was studied. To find the effects of different factors on the formation and destruction of spiral waves, the external current, the coupling strength, and the external stimulus amplitude were varied. It has been observed that all of these parameters have significant impacts on the spiral waves.

As we all know, the nervous system has a large number of neurons with different biological structures and functions. Most memristive neurons and neural network models consider only identical neurons and a single biological function (Hou et al. 2021; Tabekoueng Njitacke et al. 2020; Goetze and Lai 2021; Joshi 2021; Mersing et al. 2021; Wouapi et al. 2021). Then, the previous literature is devoted to investigating the electrical activities of coupled neurons with either identical neurons or non-identical neurons governed by the same mathematical model. However, the dynamical behaviors of coupled neurons with different mathematical models were rarely reported in previous publications (Li et al. 2021). Thus, more different memristive neuron and neural network models need to be developed based on the different biological neuronal systems as well as different mathematical models (Lin et al. 2021).

Very recently, Li et al. (Li et al. 2021) reported the dynamics of coupled neurons made of HR and FN models under the magnetic field effect. Studies of the coupled neurons showed the phenomenon of firing pattern coexistence as well as phase synchronization. Before performing this work on coupled neurons consisting of HR and FN models under the effect of the magnetic field, it would have been interesting to see what happened to the original 2D HR neuron coupled to the original 2D FN neuron by asymmetric electrical synapses. This later questioning justifies the interest of this work since in a normal brain it is difficult to have a homogenous coupling between neurons as presented in Refs. (Tabekoueng Njitacke et al. 2020; Bao et al. 2019a; Ren et al. 2017; Njitacke et al. 2020).

It is well known that during the metabolic process of the neuronal system or biological system, energy is consumed so that neurons can save normal and continuous electric activity (Harris et al. 2012). The liberation and storage of energy of neuronal models can be estimated. Then it is interesting to detect the energy transmission and release dependencies on the electric activity mode in those neuronal models. In several works focused on the dynamics of neurons, Hamilton’s statistical function is usually exploited to the energy consumption by exploiting Helmholtz’s theorem. The said energy can be also exploited to explain some biological phenomena such as synchronization regime, spiking state, bursting state. This discussion on the energy of the neurons enables us to see the main drawback of the previous works addressed on coupled neurons, and it is solved in this work.

Then, in the present contribution, we have investigated the complex dynamics of a simple heterogeneous neural network made of the coupling between a 2D HR neuron model with a 2D FN neuron model through asymmetric electrical synapses. In addition, the electrical activity and the Hamilton energy in the coupled neurons are analyzed by considering the simultaneous effect of the asymmetrical electrical coupling. Based on Helmholtz’s theorem, the Hamilton energy function H and its derivative in terms of membrane potentials, recovery variables, and parameters of the model are provided (Lu et al. 2019; Ma et al. 2017; Xin-Lin et al. 2015). This work is very distinct from the previous research works, which focus on the homogeneous coupling between neurons.

The plan of this work is presented as follows: In Sect. 2, the heterogeneous model of the coupled neurons is established. In Sect. 3, the dynamics of the coupled neurons model is studied and its Hamilton energy analysed. In Sect. 4, Pspice implementation of the coupled neurons is realized. In Sect. 5, the microcontroller implementation of the coupled neurons is addressed using the STM32F407ZE microcontroller development board. Finally, we summarize the paper in Sect. 6.

The model and its basic dynamics studies

Model description

It is well known that biological neurons can develop several types of firing activities to reproduce brain dynamics. To achieve that goal several mathematical models of functional neurons have been proposed in the literature. Among others, we have the Hodgkin–Huxley neuron (Hodgkin and Huxley 1952), Hindmarsh–Rose (HR) neurons (Hindmarsh and Rose 1982), Morris–Lecar neuron (Tsumoto et al. 2006), Fitzhugh-Nagumo neuron (Izhikevich and FitzHugh 2006), Izhikevich neuron (Izhikevich 2003b), Hopfield neural networks (Njitacke et al. 2020a, 2019; Njitacke and Kengne 2018, 2019), and so on. Let us recall that reduced from the Hodgkin-Huxley neuron model (Hodgkin and Huxley 1952), the 2D HR neuron model was first proposed in 1982 (Hindmarsh and Rose 1982) and, is described as:where , and respectively indicate the membrane potential, the spiking variable, and steady current. Besides another 2-D simplification of the Hodgkin-Huxley (HH) model leads to the FitzHugh-Nagumo (FN) model described by Izhikevich and FitzHugh (2006)where , and respectively indicate the membrane potential, the retrieval variable, and the magnitude of stimulus current. In the brain as well as in the neural system, coupling plays an essential role in the processes of generation of particular rhythms and formation of memory. However, in all the previous works, the authors have considered the coupling between identical neurons. Although in practice, the coupling of neurons can be done between different families of neurons having different mathematical models. Motivated by this observation, the bidirectional coupling between the HR neuron and the FN neuron through electrical synapses of Fig. 1 is considered in this work.

Fig. 1

Bidirectional coupling between the two different families of 2D neurons

Based on that figure, the mathematical model of the coupled neurons is obtained as:

Here and are the potentials of the membrane in HR neuron and FN neuron respectively, and are retrieval variables related to a fast current of either or , and represent exterior input currents and and indicate the strengths of electrical coupling utilized as command parameters. Let us stress that in the considerations of this model the parameters are all positive and defined as: , , , ,, , , , , and .

Steady states and their stabilities

The steady states of our coupled neurons are obtained by solving the equation thus we obtain

For the considered coupled neurons, the steady states are numerically determined aswhere and can be obtained graphically from the intersection between the following curves

Considering some selected values of the synaptic weight the steady states of the coupled neurons model are shaped in Fig. 2. From this figure, it is evident that the coupled neurons model possesses only one steady state since both curves display one intersection point. It is also observed that the stability of the steady-state depends on the control synaptic weight . To have an idea of the stability of the coupled neurons, let consider the following Jacobian matrix:

Fig. 2

Evolution of the functions defined In Eq. 6 and Eq. 7 showing the steady states of the coupled neurons at the intersection for some selected values of synaptic weight

The eigenvalues at the steady-state are computed from the following equation:

Using some selected values of the electrical synaptic weight , the steady states and their stability are summarized in Table 1. From Table 1 it is evident that the coupled neurons can display either hidden dynamics associated with a stable equilibrium point or self-excited dynamics associated with unstable equilibria (Njitacke et al. 2020a; Lai et al. 2018b; Tsafack et al. 2020).

Table 1

Steady states for some values of with their corresponding eigenvalues and stabilities

Electrical synapse value	Steady states	Eigenvalues	System stability
m_2=0.5	P_01=0.3311,0.4519,-0.8186, - 0.0607	-5.6884 + 0.0000i-0.1190±1.5769i-0.0758 + 0.0000i	Stable
m_2=0.5301	P_02=0.3417,0.4162,-0.7904, - 0.0255	-5.3093 + 0.0000i-0.0938±1.5925i-0.0770 + 0.0000i	Stable
m_2=0.75	P_03=0.4123,0.1500,-0.5861,0.2299	-2.9905 + 0.0000i0.0661±1.6547i-0.0913 + 0.0000i	Unstable
m_2=0.9	P_04=0.4525, - 0.0238,-0.4554,0.3932	0.1328±1.6399i-1.8641 + 0.0000i-0.1146 + 0.0000i	Unstable
m_2=1.0	P_05=0.4757, - 0.1315,-0.3753,0.4934	0.1477±1.6126i-0.1466 + 0.0000i-1.2761 + 0.0000i	Unstable
m_2=0.54	P_06=0.3452,0.4042,-0.7812, - 0.0140	-5.1889 + 0.0000i-0.0855±1.5975i-0.0774 + 0.0000i	Stable
m_2=0.523	P_07=0.3393,0.4244,-0.7971, - 0.0339	-5.3981 + 0.0000i-0.0995±1.5891i-0.0767 + 0.0000i	Stable

Dissipation property

The model of the coupled neurons is dissipative if its divergence is negative i.e. . The divergence of the coupled neurons model defined in Eq. (3) also known as the volume contraction rate is evaluated as indicated in Eq. (10) to evaluate the dissipation property.where stands for the phase volume. To be sure that coupled model of neurons is dissipative, the result of the divergence should verify which imply that . This enables the volume element to be contracted and the flow to be confined into a volume element over time . In other words when , each volume involving the trajectories of the coupled neurons model contract to zero with the exponential rate . Consequently, the trajectories are all confined in a particular set of zero volume and the behavior is stabilized on a stable state in the state space (Zheng et al. 2004).

Hamiltonian energy of the coupled neurons

The field of study based on energy utilization in the brain has been considered experimentally and not theoretically due to the complex connections in the brain (Chuankui 2011). Given that there is no accurate method to establish the energy consumption and supply, Hamilton energy is usually exploited to identify the appurtenance of state on energy considering a unique neuron model or a coupled neurons model (Nabi et al. 2012). To Determine the energy of our coupled neurons , let us express its dynamical equations by Eq. (11) (Lu et al. 2019; Xin-Lin et al. 2015)

here and represent the conservative component and the dissipative component, respectively. represents the gradient matrix of an energy function . Indicates a skew-symmetric matrix and is a symmetric matrix. The Hamilton energy function can be determined as:

Thus, for the coupled bidirectional neurons given in Eq. (3), we have

From Eqs. (11) and (12), the Hamilton energy function of the coupled neurons can be expressed as follows:

The general solution of Eq. (15) is computed as:

The derivative of the Hamilton energy function with respect to the time is given by

After some algebraic manipulation, we obtained

In addition,

Then

From (19) and (22), it is easy to conclude that which validates the choice of the energy function. Regarding Eq. (16), it is evident that the Hamilton energy of the coupled neurons depends on its state variables (membrane potentials and recovery variables) and its parameters. More interestingly, it is obvious that the energy function defined in Eq. (16) depends on the external currents and as it was mentioned in some works addressed on the neurons model (Lu et al. 2019; Xin-Lin et al. 2015). Consequently, energy will be sufficient to maintain continuous electrical activities in the considered coupled neurons.

Result of numerical simulations

Dynamic behavior of the bidirectional coupled neurons

The well-known fourth-order Runge–Kutta integration algorithm is exploited in this work for various numerical investigations with a fixed step of . Two parameters Lyapunov exponent graphs when increasing respectively decreasing two electrical synapses in both directions enable to quickly explore the various firing activities that can occur into the proposed bidirectional coupled model of neurons.

As depicted in Fig. 3a, b the coupled neurons considered in this work can display several types of firing patterns among which resting patterns are characterized by , periodic patterns when and chaotic patterns; characterized by . In addition, a window of hysteretic neuronal activities can be recorded in the area called which materializes the region where the right and left diagrams are different. To have a deep and detailed understanding of what happens in our functional coupled neurons, the local maxima of the first neuron are determined. The corresponding bifurcations are computed versus the electrical synapse, which connects the FN neuron with the HR neuron, and the result is presented in Fig. 4a. Figure 4b represents the maximal Lyapunov exponent related to the neuronal activity shown by the bifurcation diagrams. Finally, in Fig. 4c, we have the variation of the volume contraction rate of the model versus the bifurcation parameter . It is trivial that for any value of we have which, justifies the dissipative nature of the proposed coupled neurons model. From Fig. 4a, b two sets of data are clearly identified; the data in blue is obtained when increasing the control parameter and the red data is obtained when decreasing the control parameter. From the blue diagram, the model is able to display three types of neuronal activity including resting patterns characterized by , periodic patterns when and chaotic patterns characterized by . In contrast, the diagrams in red enable us to identify only two types of firing patterns including periodic patterns for and chaotic patterns for . Besides, these various firing activities, a huge window of hysteretic dynamics in the coupled neurons model characterized by the phenomenon of the coexistence of several types of bifurcations is found.

Fig. 3

Two parameters Lyapunov exponent of the proposed model when the coupling parameters and are simultaneously varied in both direction. Initial conditions are . The left panel a is computed by increasing the control parameters and the right panel b is computed by decreasing control parameters

Fig. 4

Bifurcations of the coupled neurons with respect to the variation of the synaptic coupling in a with the corresponding diagram of maximum Lyapunov exponent in b. The evolution of the volume contraction rate of the coupled neurons versus is depicted in c. Initial conditions are with . The diagrams in red are obtained decreasing the synaptic coupling while those in blue are obtained when increasing synaptic coupling

Hidden coexistence of multiple firing patterns

Multistable dynamics have been investigated in some functional neuron models (Ngouonkadi et al. 2016; Pisarchik et al. 2018; Bao et al. 2019a, b). This phenomenon corresponds to the coexistence of several firing patterns in the state space of a given model of single or coupled neurons. Recall that in practice the initial condition of the neuron can be affected by the environment such as the electric field, the magnetic field (Zhang et al. 2020), the temperature (Xu et al. 2020) as well as the light based on the photosensitivity of the neuron(Liu et al. 2020). The coexisting patterns exhibited by a neuron can be hidden (Bao et al. 2019b) when the neuronal activity comes from either a model with a stable rest point or a model without rest points. Self-excited (Bao et al. 2019a) when the neuronal activity comes from a model with unstable rest points. The phenomenon of coexisting bifurcations exhibited by the coupled neurons model is studied in this work through an enlargement of the window of hysteretic behavior identified on the graphs of Fig. 4. This enlargement is depicted in the graphs of Fig. 5.

Fig. 5

Enlargement of the bifurcation diagrams of Fig. 4a and the equivalent graph of the maximum Lyapunov exponent of Fig. 4b. Initial conditions are

The simultaneous study of the bifurcation diagrams of Fig. 5a and the corresponding graphs of the maximum Lyapunov exponents of Fig. 5b allows us to identify many types of activities that coexist. Among others, the coexistence of resting patterns () with periodic patterns (), the coexistence of resting patterns () with chaotic patterns () and the coexistence of periodic patterns () with chaotic patterns ().

For example, when , the coupled neurons display the coexistence of periodic spiking and chaotic bursting, as shown in Fig. 6 with the 2D projection of the HR neuron variables. While the coexistence of periodic and chaotic spiking behavior is illustrated in Fig. 7 using the 2D projection of the FN neuron variables as an argument. These coexisting patterns involving periodic and chaotic patterns are supported using both phase portraits in (a) and the time evolution of the membrane potential of each neuron in (b), both for Fig. 6 and Fig. 7. For the set of neurons parameters enabling the neuronal activities of Fig. 6 and Fig. 7, the fixed point of the model is given by , and the eigenvalues associated with this rest point are given by As the coupled neurons model presents eigenvalues with two complex conjugate roots showing negative real parts, and two negative real roots, we can conclude that the model is stable for the considered set of parameters consequently; the coexistence of firing patterns found is hidden. When , the functional model of the coupled neurons under the considerations displays the coexistence of the resting pattern with a chaotic bursting as presented in Fig. 8 using the 2D projection of the HR neuron variable. Furthermore, the coexistence of the resting pattern with a chaotic spiking of Fig. 9 is also obtained for the same set of parameters but using the 2D projection of the FN neuron. These coexisting patterns are further supported using phase portraits as well as the time evolution of the membrane potential of each considered neuron. On the phase space, the resting pattern is characterized by a dot; while on time series is characterized by a straight line.

Fig. 6

Coexistence of two different types of patterns from the HR neuron including, a periodic spiking pattern (blue) and a chaotic bursting (red) for using two different initials conditions as depicted on the graph of phase portraits of Fig. 6a and their corresponding time series of Fig. 6b

Fig. 7

Coexistence of two different types of patterns from the FN neuron including, a periodic spiking (blue) and a chaotic spiking (red) for using two different initials conditions as depicted on the graph of phase portraits of a and their corresponding time series of b

Fig. 8

Coexistence of two different types of patterns from the HR neuron including, a resting-state (blue) and a chaotic bursting (red) for using two different initials conditions as depicted on the graph of phase portraits of a and their corresponding time series of b

Fig. 9

Coexistence of two different types of patterns from the FN neuron including, a resting-state (blue) and a chaotic spiking (red) for using two different initials conditions as depicted on the graph of phase portraits of a and their respective time series of b

For these coexisting patterns, the fixed point of the model is given by and the eigenvalues associated given by . For this set of parameters, it is clear that the coupled neurons model presents eigenvalues with two complex conjugate roots showing negative real parts, and two negative real roots, and then it is evident that it is stable consequently; display the phenomenon of the coexistence of hidden firing patterns.

Attraction basins associated with other initials

The basin of attraction corresponds to the set of initial values leading to each of the various coexisting firing patterns discovered in this work (Kengne et al. 2016a, b; Kengne et al. 2015). When coexisting behaviors involving periodic and chaotic patterns are captured using both phase portraits and time evolutions as arguments. The basin of attraction of each neuron model (either HR or FN) has been determined as presented in Fig. 10a, b. Figure 10a represents the set of initial conditions; in the plane that enables to obtain; the coexistence of the hidden patterns from the coupled neurons using the HR neuron when the initial conditions of the FN neuron model are all set to zero. While Fig. 10b displays the set of the initial conditions, in the plane that enables the capture of the coexistence of the hidden firing patterns of the coupled neurons from the FN neuron when the initial conditions of the HR neuron are all fixed at zero. When the functional model of the coupled neurons under consideration displays the coexistence of a resting pattern with a chaotic pattern. The set of the initial conditions; that enables to have each coexisting hidden firing activity are shown in Fig. 10c for the HR neuron and Fig. 10d for the FN neuron. From these figures, the basin of periodic or resting activities and chaotic activities are respectively marked in blue and red while the yellow color stands for the unbounded activity. These attraction basins clearly demonstrate that the distribution of the multistability in each neuron of the coupled model is associated with the initial conditions of the other neuron. In addition, the absence of the riddled structure on each basin of attraction shows that each attraction basin has a specific domain associated with each coexisting pattern thus, the occurrence of extreme events is excluded in the introduced model of coupled neurons (Chaudhuri and Prasad 2014).

Fig. 10

Basin of attractions, given the set of initial conditions associated with each of the firing patterns depicted in Fig. 6–9 for for a and b and for c and d. The red color is associated with the chaotic pattern; the blue color is associated with periodic or resting patterns, while the yellow color corresponds to the unbounded motion of the coupled neurons

Circuit implementation of the coupled neurons

In this section, we want to confirm the hidden firing patterns exhibited by the coupled neurons considered in this work. To achieve this goal, an analog circuit based on operational amplifiers TL082, analog multipliers AD633JN, resistors, and capacitors is built as presented in Fig. 11. The power supply of the operational amplifiers is symmetric. The circuit build for the coupled neurons is composed of three main blocks: one block is an emulator of the HR neuron (in the red enclosure), another block serves as an emulator of the FN neuron (in the green enclosure) and the third block plays the role of coupling between the above-mentioned blocks (in the blue enclosure). Using the well-known Kirchhoff’s electrical circuit law, the circuit equations of the coupled neurons are given as:

Fig. 11

Analog circuit implementation of the bidirectional coupled neurons

In the set of the nonlinear differential equation obtained in Eq. (23) and represent the state evolution of the membrane potential and the recovery variable respectively. indicate the outer parameters and the sources exploited for the generation of external stimulus.

Considering , ,, .

The comparison between Eq. (23) and Eq. (3) yields to:

, , , , , , , , , , , , , , .

When performing PSPICE simulations, the coexisting behaviors found in the previous section have also been reported. Figure 12 shows the coexisting time evolutions of the membrane potentials of the HR neuron in (a) and the one of the FN neuron in (b). From these figures, several kinds of coexisting behaviors are exhibited by coupled neurons. For example, the HR neuron displays the coexistence of the resting pattern with the chaotic bursting behavior while the FN neuron displays the coexistence of the resting pattern with the chaotic spiking behavior. All these coexisting patterns are obtained for the same group of circuit values but using different initial values of the capacitors. Good accordance is observed between numerical analysis (on the right-hand side of Fig. 12) with the PSPICE results (on the left-hand side of Fig. 12).

Fig. 12

Simulation issues of the analog circuit under PSPICE environment (left-hand side) showing coexisting time evolutions of the membrane potential of the coupled neurons model with their MATLAB equivalents (right-hand side). Those in a are associated with the HR neuron while those in b are associated with the FN neuron. Initial values of voltages are for resting behaviors and for the chaotic bursting and spiking behaviors obtained for for

Microcontroller based implementation

In this section, the theoretical and numerical dynamics will be reproduced practically using a microcontroller development board (Nestor et al. 2020). A microcontroller is a digital environment system comprising encapsulated in a single case: a microprocessor, memories (Flash, EEPROM, RAM, etc.), and I / O peripherals. This digital computer operates sequentially and can reach operating frequencies of several hundred megahertz. Although it has some drawbacks, the digital computer has several advantages over the analog computer: it is robust, reprogrammable, compact, flexible in setting the parameters and the initial conditions of the oscillator to be implemented with precision, just to name a few. The STM32F407ZE development board was selected as our digital computer due to its important resources, as it is mentioned in Table 2.

Table 2

Some important resources of the STM32F407ZE microcontroller for the implementation

Core	ARM Cortex™-M4 32-bit RISC Core
	Floating point unit (FPU)
	Frequency up to 168 MHz
	DSP instructions
	Adaptive real-time accelerator (ART)
Peripheral	Flash memory 512 Kb
	SRAM memory 192 + 4 Kb
	2 × 12-bit D/A converters

The synoptic diagram of our microcontroller computer is given in Fig. 13. This configuration has advantages in terms of cost, size, and gain in energy consumption.

Fig. 13

Block diagram of the digital computer with microcontroller (Black Board STM32F407ZE)

The implementation of the system consists of coding the flowchart of Fig. 14 in an evolved language (here C/C + +) using an integrated development environment (Arduino 1.8.13 in our case). The 4-byte real numbers (float) were used as the data format during experiments. The experimental set-up is presented in Fig. 15 using parameters and initial seed of Fig. 6 and Fig. 7.

Fig. 14

Flowchart of the 4th order Runge–Kutta integration method implemented with the STM32F407ZE

Fig. 15

Experimental setup of a digital computer based on a microcontroller and visualization of the signals using a digital oscilloscope

Following the parameters and initial conditions of Fig. 6, experimental coexisting attractors in several 2D projections, planes are obtained and results are illustrated in figure Fig. 16. High similarity can be observed between the numerical and experimental results.

Fig. 16

Comparative presentation of the experimental curves (first column and third column) with numerical curves (second column and fourth column). The first two columns are periodic curves while the last two are chaotic curves. Initial conditions and parameters are those of Fig. 6

Conclusion

In this contribution, a model of bidirectionally coupled neurons consisting of HR and FN neurons was considered and analyzed in depth. The equilibrium point of the coupled neurons was investigated. Their stability has revealed that for some values of the electrical connection, the model under consideration can display either self-excited firing patterns or hidden patterns. In addition, hidden phenomena such as the coexistence of chaotic bursting with periodic spiking, the coexistence of chaotic spiking with period spiking, the coexistence of chaotic bursting with a resting pattern, and the coexistence of chaotic spiking with a resting pattern are also found for some sets of parameters. For all the phenomena reported, the Hamiltonian energy of the model has been determined. It enables us to estimate the amount of energy released during the transition between various electrical activities. A series of PSPICE simulations have been carried out based on the analog circuit of the coupled neurons model to support our numerical results. Finally, a microcontroller implementation based on the STM32F407ZE development board has been used to further support our numerical results of the multistability. Since in real neurons, it is difficult to have the same homogenous dimension in a coupling between several neurons, for future work we are going to explore the complex dynamics of two different families of neurons having different dimensions.