Hamiltonian energy computation and complex behavior of a small heterogeneous network of three neurons: circuit implementation

Brain functions are sometimes emulated using some analog integrated circuits based on the organizational principle of natural neural networks. Neuromorphic engineering is the research branch devoted to the study and realization of such circuits with striking features. In this contribution, a novel small network of three neurons is introduced and investigated. The model is built from the coupling between two 2D Hindmarsh–Rose neurons through a 2D FitzHugh–Nagumo neuron. Thus, a heterogeneous coupled network is obtained. The biophysical energy released by the network during each electrical activity is evaluated. In addition, nonlinear analysis tools such as two-parameter Lyapunov exponent, bifurcation diagrams, the graph of the largest Lyapunov exponent, phase portraits, time series, as well as the basin of attractions are used to numerically investigate the network. It is found that the model can experience hysteresis justified by the simultaneous existence of three distinct electrical activities using the same set of parameters. Finally, the circuit implementation of the network is addressed in PSPICE to further support the obtained results.


Introduction
Brain functions are sometimes emulated using some analog integrated circuits based on the organizational principle of natural neural networks (NNN). Neuromorphic engineering is the research branch devoted to the study and realization of such circuits with striking features [1]. Some of these circuits are obtained from some mathematical models of the neurons that have been developed. Hodgkin-Huxley neuron [2], Izhikevich neuron [3], Morris-Lecar neuron [4], 2-D Hindmarsh-Rose (HR), 3-D-HR model [5,6], FitzHugh-Nagumo (FN) model [7], Hopfield neural network model [8][9][10][11], Chay model [12], and the Rulkov model [13], just to name a few, are among those neurons. Recall that in brain dynamics, neurons have a significant role in the brain's recording, selection, storage, learning, thinking, and data transfer processes. From the work addressed to a single neuron, it has been found that it can exhibit resting dynamics, spiking, and bursting dynamics [14]. Numerous studies on neuron dynamics have looked at single neurons, two-coupled neurons, and neural networks. Considering the case of single neurons, several studies have been done. Among others, [15] introduced a novel equilibrium-free HR model using memristive electromagnetic induction. During their investigations, the authors found that the proposed model of the neuron was to exhibit the phenomena of hidden homogeneous extreme multistability. That phenomenon is characterized by the simultaneous existence of an infinite number of firing activities of the same shape but at different positions. In [16], 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 could exhibit several types of complex phenomena, including symmetry breaking, period doubling, reverse period-doubling, and crisis scenarios. Hou et al. [17] estimated the neuron's electrical activity in a depolarization field. In their investigation, they detected a spike in the firing interval of the neuron when the stimulus strength increased. In addition, the neuron-firing pattern was able to switch from a single burst to an intermittent multimodal burst in the presence of an electric field as the amplitude and frequency of the electric field increased.
Kafraj et al. [18] proposed and studied a modified Izhikevich model to characterize the activity of neurons in the presence of electromagnetic induction and noise. Their study revealed several complex phenomena in the proposed model. They also found that the firing activities of the neuron were caused by an outer field, whereas it was at rest without the effect of that outer field. Some intriguing works on neurons have also designed models that take into account natural phenomena such as the effect of light [19], temperature [20], and sound [21]. Very recently, Xu et al. introduced and investigated a memristor-based novel model of non-autonomous Rulkov neuron. Afterward, the results of their numerical simulations, as well as experimental investigations, revealed that the proposed memristive model was able to exhibit extreme multistable behavior [13].
Concerning coupled neurons or networks, Pisarchik et al. [22] have explored the dynamics of two identical coupled 3D HR neurons. The model is built by exploiting a bidirectional asymmetric electrical synapse. Njitacke et al. investigated successively the bidirectional coupling between two non-identical 2D HR neurons in the absence and the presence of electromagnetic radiation [23,24]. Karthikeyan et al. [25] have introduced a modified Hindmarsh-Rose neuron model having a fractional-order threshold magnetic flux. Their investigations were conducted considering the effect of a magnetic field. Besides, the development of the spiral waves in the suggested network was explored. The variation of the parameters of the proposed network enabled them to see how they affected the birth and death of spiral waves.
Baran et al. [26] investigated the phenomenon of spike-timing-dependent-plasticity (STDP) learning using chemical and memristive synapses on two coupled HR neurons. Their results were obtained based on numerical simulation, and they showed that in neuromorphic studies of nonlinear processes like chaos and hyperchaos, the memristive device might be employed as an alternative synapse. In [27], using a sigmoid function as the recovery variable, the authors investigated the complex dynamics of a set of coupled FitzHugh-Nagumo neurons. The coupling connection was obtained using electromagnetic flux coupling. Their analyses revealed that for higher magnetic coupling strengths, spiral waves were destroyed, giving birth to the periodic regimes in the individual nodes. The synchronization process and chimera state were analyzed and addressed in a set of FitzHugh-Nagumo neurons under the photosensitive effect [28]. These chimera states were also demonstrated in a thermosensitive FitzHugh-Nagumo neuronal network [29]. In [30], considering the effects of electrical synapses and autapses, the authors have explored the behaviors of the considered network of Hindmarsh-Rose neurons. Using numerical simulations, they were able to show the occurrence of the chimera states in the considered network. More interestingly, traveling chimera patterns have been demonstrated in a 2-D HR neuronal network by Simo et al. [31]. A set of 3D HR neurons under an electric field has been investigated in [32]. As a result, traveling chimera states and multicluster oscillating breathers could be found in the absence of the field, whereas chimera states, multichimera states, alternating chimera states, and multicluster traveling chimeras could be found in the presence of the field.
In all of these works that address neurons, twocoupled neurons, or networks of neurons, energy is required for each firing activity as well as the transition between electrical activities. Indeed, during the metabolic process of the neuronal system or biological system, energy is consumed so that neurons can save normal and continuous electric activity [33]. The liberation and storage of energy of neuron models can be estimated. Then, it is interesting to determine the amount of energy released during transitions between the electric activity of neurons and neuron networks [31,34]. In several works focused on the dynamics of neurons, using Helmholtz's theorem, a Hamilton's function is derived and used to calculate energy consumption. The previously mentioned energy can also be used to explain biological processes like bursting regimes, spiking states, and synchronization states [31,34]. This discussion on energy is generally focused on a single neuron or homogeneous neural network. Recall that this Hamiltonian method has been widely used by Ma and collaborators to estimate the energy released by neurons during electrical activity [34][35][36].
So, from all the above-mentioned work, it is obvious that several scientific contributions have been devoted to the research and study of a large number of models of neurons or neural networks, taking into account environmental effects. Among them, we have the effect of external current [22,23], the magnetic field effect [37][38][39], the electric field effect [32,40], photosensitive effect [19], piezoelectric ceramics or auditory effect [21,35], and thermosensitive effect [20], just to name a few.
It is well known that the nervous system consists of a multitude of neurons. Some of them have specific biological functions and structures. Most works on coupled neurons or neural network models only considers identical neurons with a single biological function [17,23,39,[41][42][43]. Then, from the previous literature, it is obvious that many research efforts have been devoted to the study of the behaviors of coupled neurons. These coupled neurons are built with either identical or nonidentical neurons described by the same physical equation. In contrast, previous research has rarely studied the dynamic behavior of coupled neurons, as described by different mathematical models [44]. Therefore, more diverse neural network models, as well as memristive neural networks, must be constructed based on neurons with different mathematical models [45]. This is why, in this contribution, we propose and investigate the complex behavior of a small heterogeneous network of three neurons. That small network is constructed by coupling two traditional 2D HR neurons through a 2D FN neuron. Using the well-known Helmholtz theorem, the energy of the network is estimated to justify the nonlinear behaviors found during the analysis of the proposed model. So the remaining part of this work is structured as follows: In Sect. 2, the mathematical model of the proposed small network is established and its Hamilton energy is calculated. In Sect. 3, numerical simulations are used to investigate the nonlinear behavior of the coupled network. In Sect. 4, a PSPICE circuit of the network of three neurons is realized. Lastly, in Sect. 5, we summarize the paper.
2 Design and analysis of the proposed network

Presentation of the model
The brain as a nonlinear system is made up of an interconnection of several types of neurons with several functions. These various types of neurons, as well as neuron networks, are explored based on their mathematical models. The Hodgkin-Huxley neuron [30], the Hindmarsh-Rose (HR) neuron [31], the Morris-Lecar neuron [17], the FitzHugh-Nagumo neuron [32], the Izhikevich neuron [33], and the Hopfield neural networks [34][35][36][37] are among the mathematical models. In addition to the preceding research, the complicated functioning of the brain has been examined using homogeneous coupling (the pairing of similar neurons) between several of these quoted neuron models. However, a practical brain model is made of the interconnection of neurons with different functions and different mathematical models. Those neurons can also be subjected to perturbations induced by fields (magnetic or electric), temperature, light, and sound. This is why, in this contribution, the small network in Fig. 1 is proposed. It is made up of the bidirectional coupling between two HR neurons connected through a FN neuron.
Based on that coupling scheme in Fig. 1, the mathematical model of that small network is given in Eq. (1).
Here x i are the potentials of the membrane in HR neuron and FN neuron, respectively, y i are retrieval variables related to a fast current of either Na þ or K þ , i i represent exterior input currents and m ij indicate the strengths of electrical coupling utilized to connect Let us stress that in the considerations of this model the parameters are all positive and defined as:

Rest points of the network
The rest points are very important for the analysis of the stability of the considered network. As a result, they allow justification if the dynamics of the considered network are self-excited or hidden. We recall that self-excited dynamics is associated with a system having unstable rest points, whereas hidden dynamics is associated with a system without rest point or with stable rests point generally. To determine the rest points of the considered network, let us consider _ x i ¼ _ y i ¼ 0 and then, Eq. (2) is obtained.
After computing the equilibrium of each state variable of the coupled neurons with the software Maple 18, we discover that they are in the form a þ jb, which is a complex number (two cases are shown in the appendix), implying that this small network has a hidden dynamics.

2.3
The energy needed for firing activities.
The energy produced during the transition of electrical activity is very important and has been estimated in neural dynamics [46]. Given the lack of an accurate method for determining energy consumption and supply, Hamilton energy based on Helmholtz's theorem is commonly used to identify the appurtenance of state on energy in the context of a neuron or networks of neurons [47]. To determine the energy of our small network made of three neurons F x ð Þ, let us express its dynamical equations by Eq. (3) [48,49].
For the small network of Eq. (1), the conservative component and the dissipative component can be expressed as:  Fig. 1 The topological configuration of the small network of three neurons used for this study 3   2   6  6  6  6  6  6  6  4   3   7  7  7  7  7  7  7  5 ð6Þ Using Eqs. (5) and (6), the energy H x 1 ; y 1 ; x 2 ; y 2 ; x 3 ; y 3 ð Þof the small neuron of three neurons will verify the following equation: The general solution of Eq. (7) can be expressed as: The Hamilton energy function's derivative with respect to time is denoted by It is simple to demonstrate that after some algebraic manipulation, Remark that Eq. (10) is also equal to rH T F d since with From (10) to (12), it can be concluded that dH dt ¼ rH T F d x ð Þ which reinforces the selection of the energy function. Regarding Eq. (8), it is evident that the Hamilton energy of the small network of neurons is a function of all state variables and all parameters. Even more fascinating, it is clear that the energy function defined in Eq. (8) is clearly dependent on all the external currents i i as suggested in [48,49]. As a result, there will be enough energy to sustain continuous electrical activity in the small network of neurons under consideration.

Numerical analysis
In this section, our main objective is to investigate the global dynamics of the proposed network of three neurons using the fourth-order Runge-Kutta algorithm with a fixed step time of 5 Â 10 À3 . The parameters as well as the electrical coupling weight of the network are selected in the extended precision mode for better accuracy of the calculations. Nonlinear analysis tools such as two-parameter Lyapunov exponent graphs, the graph of the Hamiltonian energies, bifurcation diagrams, phase portraits; time series as well as the basin of attractions are used to characterize the complex dynamics of the network.

Global dynamics of the network
To investigate the global dynamics of the proposed network, electrical coupling weights will be used as bifurcation control parameters. Based on those electrical coupling weights, some two-parameter Lyapunov exponent graphs are provided in Fig. 2. These diagrams have been computed when increasing (left side) and decreasing (right side) both parameters. This computational technique represents a nice way to cover a wide range of the parameters and highlight hysteretic dynamics. From Fig. 2(a) and (b), the twoparameter diagrams on the left-and right-hand sides are almost the same, which excludes the possibility of hysteretic dynamics. In contrast, from Fig. 2(c), a large region of discrepancy is observed between both diagrams. That discrepancy is justified by the hysteretic nature of the network for the set of the electrical coupling weights used for the investigations. As it can be seen on those two-parameter Lyapunov exponents, the dynamics of the proposed network is able to vary from periodic firing activities k max ¼ 0 ð Þto chaotic firing activities k max [ 0 ð Þ based on the sign of the maximum Lyapunov exponent. Remark that for each firing activity observed in those diagrams, energy, also called Hamiltonian, is consumed. So the Hamilton energy associated with each electrical activity on the left-hand side of Fig. 2 is provided in Fig. 3. That energy is obtained by varying two electrical synaptic weights simultaneously and recoding the equivalent energy function as given in Eq. (8). From those diagrams in Fig. 3, it is clear the energy of that proposed small network of three neurons strongly depends on the values of the various electrical synaptic weights that enable it to couple neurons.
To see which type of firing activity can occur during the transition between electrical activities, the bifurcation diagrams of Fig. 4 and its corresponding graph of the maximum Lyapunov exponent have been computed when varying m 32 . As it can be seen in those diagrams, the dynamics of the small network of three neurons varies from periodic to periodic firing activities through several windows of chaotic firing activities. Recall that the bifurcation diagram is computed by varying a control parameter (electrical coupling weight) and recording the maxima of the membrane potential of any neuron at each iteration of the control parameter. It can also be obtained by superimposing two diagrams when increasing or decreasing the control parameter. For this computational method, the final states of the network for an iteration are used as initial conditions for the next iteration. For some discrete values of the electrical coupling weight m 32 , phase portraits are displayed in Fig. 5, showing the phenomenon of the period-doubling bifurcation exhibited by the considered network.

Multistable behavior of the small network
Previously, it has been obtained from Fig. 2(c) that the network of three neurons considered in this work was able to display hysteretic dynamics. That phenomenon is characterized by the coexistence of several firing activities for the same set of electrical coupling weights, but using different initials for the membrane potentials as well as retrieval variables. This phenomenon of multistability involving the coexistence of bifurcations is analyzed in this work using the bifurcation diagrams of Fig. 6. From these diagrams, three sets of data, namely blue, red, and black, associated with the coexistence of the firing activities are found. The details of the choice of the initial conditions and sweeping direction used to compute these diagrams are given in Table 1.
For a discrete value of the electrical coupling weight m 32 ¼ 0:994, three types of coexisting firing activities are reported in Fig. 7 using the state variables of all the neurons in our proposed network. Among the three coexisting firing activities, two are periodic, while one is chaotic. The bursting nature of that chaotic coexisting activity is supported using the time series of Fig. 8. From those time series, it is obvious that the membrane potential of the first and third neurons evolves as the fast variables, whereas the membrane potential of the second neuron evolves as the slow variable. Remember that bursting neuron activity occurs when the activity of neurons swaps between a resting regime and redundant spiking [14]. Originally, that phenomenon was exploited to characterize the electrical activity of neurons. When the model is divided into two subsystems with two different time scales, these burstings in neurons or coupled neurons may occur.
The fast subsystem contains fast variables, and the slow subsystem contains slow variables. Then, the model can display two time scales. Initial conditions are very important for the analysis of neuronal systems. In nature, the initial conditions of those neuronal systems, including each category of neurons or coupled neurons, can be influenced by the effects of light, temperature or fields that can be either electrical or magnetic. So, when varying the initial conditions of the small network of three neurons considered in this work, the basins of attraction in Fig. 9 are obtained. These basins of attraction are separated by colors that correspond to the set of the initial conditions that enable us to obtain each coexisting firing activity in Fig. 7. From those attraction basins, it is obvious that one of the chaotic firing activities (bursting oscillations) is less than the other one of the two periodic firing activities.

Circuit realization of the proposed small network
In this section, the objective is to confirm the hidden behavior of the small network of three neurons proposed 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. 10. The power supply of the operational amplifiers is AE15V symmetric. Remember that the original goal of neuromorphic engineering was to create analog integrated circuits based on the natural neuronal network (NNN) organizational principle to mimic some of the brain's functions. [1]. Using the well-known Kirchhoff's electrical circuit law, the circuit equations of the proposed network are given as: In the set of the nonlinear differential equations obtained in Eq. (13), X i and Y i represent the state evolution of the membrane potential and the recovery variable, respectively. Considering the comparison between Eqs. (13) and (1) yields: When performing PSPICE simulations, the phenomenon of the period-doubling bifurcation route to the chaos was also obtained for some discrete values of the control parameter R m 32 . For example, period-1 limit cycle was obtained for R m 32 ¼ 10kX, and period-2 limit cycle was obtained for R m 32 ¼ 10:52kX. The period-4 limit cycle was obtained for R m 32 ¼ 10:83kX and finally a chaotic attractor for R m 32 ¼ 11:52kX as its shown in Fig. 11. The coexisting behaviors found in the previous section have been also reported. Figure 12 shows the phase portraits of the three coexisting stable states exhibited by each neuron of the proposed small network using three different initial conditions. Finally, the chaotic bursting behavior found among the coexisting stable states has been further supported using the time series of Fig. 13.
There, it is obvious that state variables V X1 ð Þ and V X3 ð Þ stand for the fast subsystem with state variable V X2 ð Þdenoting the slow subsystem. Good accordance  Þplan for discrete value of m32: a m 32 ¼ 1 is a limit cycle of period-1, b m 32 ¼ 0:95 is a limit cycle of period-2, c m 32 0:923 is a limit cycle of period-4 and d m 32 Þ0:868 is chaotic attractor