Dynamics of an array mechanical arms coupled each to a Fitzhugh-Nagumo neuron

In this work, an array of electromechanical systems driven by an electrical line of Fitzhugh-Nagumo neuron is analyzed. It is shown that a single electromechanical system can display diﬀerent dynamical behaviors such as single and multiple pulse generation, transient chaos, permanent chaos, and antimonotonicity according to the system parameters. In the case of an array of the electromechanical system constituted of a series of coupled discrete Fitzhugh-Nagumo neuron, the numerical simulation shows that as the action potential ﬂows in the discrete array, each electromechanical system executes a pulse-like motion coming at each resting state as the electrical signal passes the node. The electromechanical system analyzed can be seen as a model for multi-periodic actuation processes or a leg model in a millipede system. Furthermore, this line can also carry an envelope of action potential and can be useful for various kinds of information processing systems.


Introduction
The study of wave propagation in systems of excitable biological cells is an important aspect of today's cardio-physiology and neurophysiology. One typical example is propagation failure, leading to failure of these systems; this can be fatal in the case of the cardiac action potential. Biological neurons intervene as an electrical input-output membrane voltage models to produce the spikes due to their membranes which have voltage dependent ionic channels [1]. It is shown that biological neurons can function either excitatory or inhibitory depending on the excitability threshold of the ion channels [2,3]. Inspired by the biological process, scientists have made efforts in neuroscience to understand through the numerical simulation the electrical behavior of a real biological neuron. Among these models used, Hodgkin-Huxley [4,5], Hindmarsh Rose [6,7], Morris Lecar [8,9], and Fitzhugh Nagumo models [10,11], to name a few, exhibit an interesting dynamics investigation neurons. Being most used because of its simplicity and small equivalent circuit [13], Fitzhugh Nagumo model is a relaxation oscillator [14]. It is also well known as a generalization of the Van der pol equations that are used to model the behavior of excitable systems [15]. The excitability represents the fact that from his resting state if a stimulus (small perturbation) is applied on the excitable cell for a brief time interval, the cell membrane voltage reaches a threshold value and generates an action potential before returning to rest [16,17,18]. Generally, in their dynamical behavior, Fitzhugh Nagumo displays a solution that converges to a fixed point or a limit cycle [19,20]. The latter representing the phase space is described in the framework neuron as the periodic signal of the membrane potential as it evolves in time.
The investigation of the behavior of the limit cycle has become an ongoing research topic [21,22,23,24]. For the actuation tasks, we will analyze a system where the Fitzhugh Nagumo neuron is used to control a mechanical arm. This is similar to the coupling between a biological neuron and a muscle. In fact, when the brain decides to put in motion a part of the body and gives the command to the motor neurons to execute this movement, it is the muscles at the end of the chain of command that ultimately contract to move the body part concerned. The axons of these motor neurons, emerging from the spinal cord, form a nerve that extends to the muscles to transmit this command. Where the tip of each axon comes into proximity with a muscle fibre, it forms a synapse with that fibre.
This device can be used to explain the mechanical movement of the arm or the leg for example. Similarly, the legs of millipedes or centipedes can be viewed as an array of Fitzhugh Nagumo neurons coupled each to an electromechanical system. It can also be used for the legs of artificial millipedes where they move their legs in a wave like motion from the front to the back. For this reason, the excitable Fitzhugh Nagumo neuron has been widely studied as an electrical nonlinear transmission line for the propagation wave phenomena [25,26,27,28] in which the nerve impulse during its propagation evolves from cell to cell by keeping its shape. It is mentioned that the electrical Fitzhugh Nagumo line can be viewed as a discrete medium [29,30,31] or as a continuous medium [32,33,34]. In our work, we will consider the discrete medium.
Inspired by the topics of ref [35,36], the main objective of this work is to analyze the dynamical behavior of an array of electromechanical systems powered by an array of Fitzhugh Nagumo neurons. The outline of the paper is organized as follows. Section 2 is devoted to the system constituted of one single Fitzhugh Nagumo neuron coupled magnetically to a mechanical arm. The corresponding equations and their dynamical behavior are presented. Section 3 is concerned with the array of electromechanical arms powered by an array of discrete excitable Fitzhugh Nagumo neurons. We analyzed the effect of the coupling strength on the mechanical arms displacement. The paper ends with the conclusions presented in section 4.
Dynamics of an array mechanical arms coupled each to a Fitzhugh-Nagumo neuron 3 2 Electromechanical arm powered by a Fitzhugh-Nagumo neuron

Description and mathematical model
We consider an equivalent circuit of Fitzhugh-Nagumo neuron coupled magnetically to a mechanical arm, as shown in Figure 1. In the equivalent circuit of Figure 1a), Is represents the external electrical stimulation current source, C is the membrane capacitance, G is the membrane conductance, L and r represent respectively the inductance and resistance of the membrane. In this model, the conductance G is the only nonlinear element, and its voltage-current characteristics is given in equation (1): where i G and v are the current through and voltage across the conductance respectively. µ 1 and µ 2 are respectively the threshold voltage and the diffusion potential of the neuron. Finally, α is a fitting parameter and is a function of the potentials of different ions present in the neuron. The Fitzhugh-Nagumo model of nonlinear conductance shown in equation (1) can be simulated by a different nonlinear electric circuit, using a tunnel diode or a nonlinear resistor with a smooth cubic v − i characteristic.
The block MS represents our mechanical arm and its internal structure is shown in Figure 1b). The coil is positioned in the air gap of the magnet and a beam is rigidly attached to the coil. In addition, a spring is added to avoid the movement of the mobile beam away from the balanced position established when the system was assembled.
The interaction of the current through the windings and the magnetic field produces mechanical vibrations of the mobile beam.
Two different strategies are responsible for electrical communication between neurons. One is the consequence of low resistance intercellular pathways, called "gap junctions". The second occurs in the absence of cell-to-cell contacts and is a consequence of the extracellular electrical fields generated by the electrical activity of neurons. In the same manner, the capacitance Cm is used to realize the coupling between the neuron and the electromechanical system. The application of the Kirchhoff laws to the circuit shown in Figure 1a) leads to the following differential equations: L di dt Lm and rm are respectively the inductance and resistance of the windings, v represents the membrane potential of the cell, namely the voltage across the membrane capacitance C, i is the current through the inductance L, and represents biologically the recovery variable related to the inactivation of the sodium channels. i m is the current flowing through the winding, u is the voltage across the capacitance C m and represents the electrical synaptic potential or coupling potential, and t is the time. The equation of motion of the mobile beam of mass m is given by: where β 0 is the damping coefficient, and K is the spring constant. is the total length of the conductor used in the winding. The term −B dx dt represents the induced voltage in the winding, while the term B i represents the Laplace force acting on the conducting wire in the magnetic field. Introducing the news variables and the following dimensionless parameters: , and q = 1 we obtain that the dimensionless equations governing the dynamics of the whole system are thus: The action of the electrical subsystem on the electromechanical part is visualized through the parameter σ while the effect of the electromechanical system in the electrical block is measured through the product dγ. If σ dγ, the electrical subsystem is viewed as a voltage source by the mechanical part. In contrast, if σ dγ, the power provided by the electrical part will be less than the power required by the mechanical part. The moderated values of σ and dγ are then required for the subsystems to influence each other.
Dynamics of an array mechanical arms coupled each to a Fitzhugh-Nagumo neuron 5

Stability analysis
Assume E (v 0 , w 0 , z 0 , u 0 , y 0 , x 0 ) is the equilibrium point of the system. From the system of equations (8), we have w 0 = v0 r , z 0 = 0, u 0 = v 0 , y 0 = 0, and x 0 = 0 while v 0 is the solution of the following third-order algebraic equation: To derive solutions of equation (9), let us consider the following parameters δ and ν and determinant ∆: (10) As odd degree polynomial, equation (9) can have one or three real solutions according to the values of the above parameters, and the sign of ∆.
If ∆ < 0, equation (9) has three real solutions, namely: The values of the parameters used in this work are the following: Using the parameters given above, the evolution of the equilibrium potential v 0 is plotted as a function of the stimulation current Is in Figure 2, and for two different values of the parameter α.
For the graph of Figure 2, we have considered = 1 m. The curve with a dashed line and the curve with a full line are obtained respectively for α = 1 and α = 2.
As predicted analytically, the system can present one equilibrium potential or three equilibrium potentials depending on the values of the parameters. As noticed from the graph, for α = 2, the case of three equilibrium points arrives here for values of the stimulation currents between 17.01 mA and 60.815 mA.
The stability of the system is determined by the eigenvalues of the Jacobian matrix at equilibrium points E given by : Its eigenvalues are solutions of the characteristic polynomial in λ of the linearization at E, namely: Here the coefficients a n , with 0 ≤ n ≤ 5 are given respectively as: We have used the Routh-Hurwitz criterion [37] to analyze the stability of the equilibrium points. Using the parameters given above, the stability boundary of the system is plotted in Figure 3a) and Figure 3b) for α = 1 and α = 2 respectively.
In Figure 3a), we have two colored regions, while in Figure 3b) we have three colored regions. If a coupled of values Is and is chosen in the black region, the system is stable around the single equilibrium point. However, if a coupled of values I s and is chosen in the white region, the system is unstable. In the blue region of Figure 3b), the system has three equilibrium points, and one of these is stable.

Limit cycle prediction
Although precise knowledge of the waveform of a limit cycle is usually not mandatory, knowledge of the existence of a limit cycle, and its approximate amplitude and Dynamics of an array mechanical arms coupled each to a Fitzhugh-Nagumo neuron In the Black area, the system is stable around the single equilibrium point while in the white area, the system is unstable. In the blue area, the system has one stable equilibrium point amount three. a) α = 1 and b) α = 2.
frequency, is a prerequisite to good system design. The limit cycle phenomenon deserves special attention since it is apt to occur in any nonlinear physical system. A limit cycle is desirable here since it provides the real situation where all the variables are not constant. We will apply the linear theory to the quasi-linearized system, and points of neutral stability are sought. Any undamped oscillations so arrived at are interpreted as limit cycles in the original nonlinear system.
If we consider the solution of equation (14) to be of the form of λ = jω, where j 2 = −1 and ω is the natural radian frequency of the system, then the following radian frequencies and the condition for self-starting are obtained: and ω 2 = a 3 + a 2 3 − 4a 1 a 5 2a 5 .
(16) During our investigations, we have found that for large values of the parameter , the system oscillates with the radian frequency ω 1 while it oscillates with the radian frequency ω 2 for small values of . With the stimulation current Is varies from 0 mA up to 200 mA, we have found the corresponding values of that satisfy the set of equations (16) and our result is presented in the graph of Figure 4a). The corresponding value of the frequency obtained from equation (15) is plotted in Figure 4b).
The curves of Figure 4 are obtained for α = 1 and the values of < 10 m have been considered. We can notice from the graph that, as the stimulation current increases, the value of that can satisfy the oscillation condition given by equation (16)

Oscillation boundaries
Using the Routh-Hurwitz coefficients, we have determined the oscillation boundaries of the system and our results are presented in Figure 3a). If the values of I s and are chosen in the black areas of Figure 3a), the system will converge towards the equilibrium point. Otherwise, the system will fall in the oscillation. For example, with = 1 m we find that the system will oscillate for 33.7 mA≤ Is ≤ 154.0 mA. To verify this result, we have consider the following three values of I s : I s = 30 mA, I s = 90 mA and I s = 160 mA. The corresponding time series of the membrane potential and the mechanical displacement are plotted respectively in Figure 5 and Figure 6. As shown on the graphs of Figure 5 and Figure 6, when the stimulation current is less than 33.7 mA or greater than 154.0 mA, the system converges toward the single stable equilibrium point. While if the stimulation current is between 33.7 mA and 154.0 mA, the system provides oscillating signals. The frequency of the signal plotted Dynamics of an array mechanical arms coupled each to a Fitzhugh-Nagumo neuron

Effect of some parameters
It is interesting at this level to see how some parameters affect the dynamics of the system, specially the amplitude and frequency of the mechanical displacement. This element can be important to fix the parameters of the system. We will first analyze the effect of the length of the conductor used in the winding. Therefore, we plot the amplitudes of the membrane potential and the mechanical displacement as a function of and this for two different values of the spring constant as presented in Figure 7. For both Figures 7a) and 7b), the full lines and the dashed lines are plotted respectively for K = 10 N/m and K = 20 N/m. As revealed by the graphs, when the system oscillates, the amplitude of the membrane potential is closed to unity while the amplitude of the mechanical displacement is a decreasing function of . Otherwise, the system converges to the equilibrium point. The mechanical displacement reaches its maximum value when 2 m. We can also notice from the graphs that the equilibrium points are not a function of the spring constant K and the length of the conductor used in the winding, and this has been observed during our analytical investigations presented previously. Finally, the graphs show that the amplitude of the mechanical displacement and the oscillation domain are decreasing functions of the spring constant. Furthermore, we analyze the dynamical behavior of the system using the stimulation current Is as the control parameter. The amplitude responses of electrical and mechanical subsystems are respectively shown in Figure 8a) and Figure 8b) for two values of . For both figures, the curves with the full line are plotted for = 1 m, while the curves with dashed lines are plotted for = 2 m. We can notice from the graphs that when the system falls in the oscillation, the amplitude of the membrane potential and the amplitude of the mechanical displacement are not affected by the input current I s . The figures also show that the amplitudes remain constant (no oscillation) as the stimulation current I s increases until a critical value from where both amplitudes increase abruptly.

Chaotic behavior
The characteristic dynamical behaviors are finally investigated by varying the type of the input current source. We have done many simulations, but a part of transient chaos [38,39], permanent chaotic behavior has not been found in the autonomous system.
Hence, the constant current source is replaced here by a sinusoidal current source of amplitude Im and frequency f defined as where I s represents the DC component and ω is the normalized radian frequency of the external source. We will study the dynamics of the system keeping constant the following parameters: In the bifurcation diagrams of Figures 9a), 9b), 9c) and 9d), chaotic states are observed. The system follows a period-doubling route to chaos and also undergoes a reverse period-doubling sequence. These forward and reverse period doubling sequences, as a parameter of the system increases in a monotone way, are called antimonotonicity. While in the bifurcation diagrams of Figures 9e) and 9f), only periodic states are observed. For I s = −12 mA, the system undergoes the sequence: For Is = −15 mA, the system undergoes the sequence: These bifurcation diagrams show the period bubble and the primer bubble respectively.
The study of the system (8) with the external current source given by equation (17) reveals the existence of chaotic behavior, following the period-doubling route to chaos. Similar results can be obtained using a current source that provides a square signal with DC component I s and frequency f . For illustration, phase portraits of the system in the chaotic state are shown in Figures 10a) and 10b). 3 Array of electromechanical arms powered by an array of Fitzhugh Nagumo neurons

Description of the model and equations
As mentioned in the introduction, the other goal of this work is to analyze the behavior of a system consisting of an electromechanical system fixed at each node of a discrete line of coupled Fitzhugh-Nagumo type oscillator. The coupling between two neighborhood elements is made by a resistance R c as illustrated in Figure 11. Dynamics of an array mechanical arms coupled each to a Fitzhugh-Nagumo neuron 13 As mentioned previously, the block MS represents the electromechanical system. Rc is an intercellular resistance that represents the coupling resistance between the cells. n = 1, 2, · · · , N , where N = 6000 is the total number of cells. Let v n and u n be respectively the voltages drop across the capacitances C and Cm in the n th cell. Still in the n th cell, i n and i mn are respectively the currents through the inductances L and L m and x n is the mechanical displacement of the mobile beam. The application of the Kirchhoff's and Newton laws to the complete system leads to the following differential equations: We will proceed as we have done in the previous section. Then, considering the dimensionless variables introduced in the previous section: v n , w n , z n , u n and x n , the overall system is described by the following nonlinear ordinary differential equations: The parameters used in the system of equations (19) have been defined in the above section and we will consider = 4 m and L = 240 mH.
Let us first recall some interesting results when the mechanical part is not connected and the membrane inductance is neglected (L = 0 H) in the equivalent circuit. The coupling term in the first equation of system (19) can be approximated with partial derivatives with respect to distance, x , assuming that the spacing between two adjacent units is small. If we assume that the voltage v varies slowly from one unit section to the other, the discrete spatial coordinate n can be replaced by a continuous one x , the network is then described by the following diffusion equation: We start the analysis by looking for wave solutions of equation (20) in the form of This assumption of a traveling wave converts the partial differential equation (20) into the second order ordinary differential equation In the above ordinary differential equation, the traveling wave speed V 0 is an unknown parameter that must be obtained by the analysis. The corresponding traveling wave solution and the traveling wave speed are given by where the new introduced parameters A and B are worth Next, let us consider the influence of the recovery term in system (19), where the inductance L and the mechanical part are included. There are no analytical expressions so far, and the results will be achieved in the following sections only by means of computer simulations. Nevertheless, the expressions obtained in equations (22) and (23) can provide some ideas. Looking at an understanding of the underlying physical processes and possible technical applications, we additionally study the influence of circuit parameters on the signal waveform.

Propagation of the signal and motion of the electromechanical arms
The set of discrete differential equations describes the propagation of electrical signal and the motion of N mechanical arms. To determine the propagation conditions, numerical solution is obtained using a fourth order Runge-Kutta algorithm with a time step ∆t = 10 −3 . The initial conditions are chosen as vn (0) = 0, wn (0) = 0, x n (0) = 0,ẋ n (0) = 0, z n (0) = 0 and u n (0) = 0 where n = 1, 2, · · · , N . The boundary conditions are considered for the first and last nodes as The excitation is performed with a constant current source applied just to the first cell, that is I sn = 0 mA for n = 1.
We will first fix the stimulation current I s1 = 10 mA and analyze the effect of the coupling resistance Rc. After multiple simulations, we have found that the action potential is created in the first cell and propagates through the other cells for 18.7 Ω Rc 123.7 Ω. Figure 12 displays the signals propagation in the discrete line and the behavior of the mobile beams when the action potential propagates in the line for Figures 12a) and 12b) show respectively the propagation of the action potentials v n and the displacement x n of the mobile beams as function of time. As shown in these Figures, the action potential shape is preserved during the propagation and each electromechanical subsystem exhibits a pulse-like behavior moving from its resting state, then increases to a maximal value and then decreases to the rest state. This observation is interesting since it indicates that the mobile beam executes an actuation work and then returns to its resting state. This profile is similar to the use of legs of artificial millipedes.
During our numerical investigations, we found that the traveling wave speed decreases as the coupling resistance increases. This fact is qualitatively in good agreement with the analytical expression of V 0 given in equation (22). For verification, the time evolution of the system at different cells are presented in Figure 13 for R c = 80 Ω.
Dynamics of an array mechanical arms coupled each to a Fitzhugh-Nagumo neuron The graphs in Figures 12 and 13 are plotted for the same values of the parameters except that the first are plotted for R c = 20 Ω and the second ones for R c = 80 Ω. We can notice that the traveling wave speed has decreased when the coupling resistance is increased from 20 Ω to 80 Ω. A part of that traveling wave speed, the curves of Figures 12 and 13 have the same behavior.

Effect of the stimulation current
We continue our investigations by analyzing the behavior of the discrete line when the stimulation current is increased to I s1 = 25 mA. We have then notice that, for certain values of the coupling resistance, an envelope of action potential propagates in the line. Our results plotted for Rc = 100 Ω are shown in Figure 14.  Figure 14a) and Figure 14b) show respectively the behavior of the nonlinear electrical line and the behavior of some electromechanical systems. A packet of three action potentials is provided by the first cell and propagates through the line. Each mechanical arm exhibits a packet of three pulse-like behavior before returning to rest.

Space-time evolution of mechanical arms
Finally, we present respectively in figure 15a) and 15b) the spatiotemporal evolution of the mechanical arm and the spatiotemporal variation between displacements of the arm of nearest-neighbor defined as shif t(t)=x i (t) − x i+1 (t).

Conclusion
We have firstly analyzed the dynamics of one single Fitzhugh-Nagumo neuron coupled to a mechanical arm. In this case, we have presented the oscillation boundaries where we have shown that for certain values of stimulation current, the system can converge towards the equilibrium point and for another values, the system can fall in oscillation or display a limit cycle behavior. The oscillation frequency and the oscillation conditions were determined. Forward period-doubling bifurcation sequences followed by reverse period-doubling sequences, as a parameter is varied in a monotone way, antimonotonicity is observed in the system. Secondly, the investigations of an array of electromechanical arms powered by an array of discrete excitable Fitzhugh-Nagumo have been done. When the first cell is excited, the nerve impulse propagates in the electrical line, and each mechanical arm displays a pulse-like behavior. It was found that the train of three pulses can also propagate through the line according to the coupling strength.