Dynamical Computations of the FitzHugh- Nagumo Equation

Okey Oseloka Onyejekwe; Okey Oseloka Onyejekwe

ISSN: 2641-3027

Archive of Biomedical Science and Engineering

Research Article Open Access Peer-Reviewed

Dynamical Computations of the FitzHugh- Nagumo Equation

Okey Oseloka Onyejekwe*

Author and article information

The Robnello Unit for Continuum Mechanics Applications and Nonlinear Dynamics, Umuagu, Oshimili South, Asaba, Delta State, Nigeria

*Corresponding author: Okey Oseloka Onyejekwe, The Robnello Unit for Continuum Mechanics Applications and Nonlinear Dynamics, Umuagu, Oshimili South, Asaba, Delta State, Nigeria, E-mail: okuzaks@yahoo.com

DOI: 10.17352/abse.000026

Received: 16 August, 2021 | Accepted: 31 August, 2021 | Published: 01 September, 2021

Keywords: Hodgin-huxley model; FitzHugh-nagumo equation; Nonlinear differential equations; Dynamical system; Bifurcation; Nonlinear stability analysis; Qualitative features; Excitation media

Cite this as

Onyejekwe OO (2021) Dynamical Computations of the FitzHugh- Nagumo Equation. Arch Biomed Sci Eng 7(1): 027-036. DOI: 10.17352/abse.000026

Abstract

The Hodgin-Huxley model is one of the most widely studied biological systems of nonlinear differential equations that is applied to explore nerve cells activities via electrical communications. In this paper we consider some numerical aspects of a simplified version of this model known as the FitzHugh-Nagumo (FHN) equation. Dynamical experiments conducted herein not only confirm those obtained from earlier studies but also facilitate a better understanding of the qualitative features of the FHN model especially those that initiate the behavior a threshold-triggered excitation media. To this end, methods of dynamical system analysis such as bifurcation and linear stability analysis are deployed to investigate the general qualitative features of an inhibitor-activator system which characterizes the FHN system of equations.

Main article text

Introduction

FitzHugh and Nagumo [1,2] were motivated by the Hodgkin and Huxley [3] theoretical description of a cell reaction based on equivalent circuit diagrams to construct a much simpler two variable model. Their two-coupled nonlinear ordinary differential equations replaced the four more realistic but cumbersome Hodgkin and Huxley (HH) type; and displayed almost the same dynamical properties and stability characteristics. On the whole, the FitzHugh-Nagumo (FHN) model falls under the broad category of reaction-diffusion; inhibitor-activator systems where the inhibitor variable is slow and the activator variable is fast. Though it was originally used to describe the action potential activity of the neurons; the coupling of excitation and diffusion make it amenable to handle the propagation of waves in excitable media that can arise in spiral wave dynamics [4]. This system of equations has therefore been found to provide a corner stone for other models with similar physical characteristics. By extension, its morphological features can also be extended to power generating systems involving water supply and irrigation or birth-death reaction models involving the spread of disease and micro-organisms, including robotic prostheses employing the transmission of signals between the central nervous systems and smart robots and more complex applications found in bio-engineering and neuronal physics [5].

The FHN model exhibits oscillations that are characteristic of strongly dissipative systems with time scale separation as is demonstrated by high amplitude response to periodic inputs in a selected frequency range. To some varying degrees, examples of these can be found in the RLC circuit, the Vander Pol equations [6] and the mass-spring-damper systems. This overriding characteristic, makes it a key ingredient and a starting point of understanding more complex phenomena as in the case of biological cells which maintain a difference of electric potential across their membrane. Many other models were developed after the HH model. These include those of Morris-Lecar (ML) and FHN [7] to mention just a few. Each of these vary in utility, efficiency and biophysical relevance as each was tailored to supplement the original HH model and achieve some specific objectives.

The inherent interdisciplinary nature of the FHN model serves as a key motivation for using qualitative techniques to understand many key ingredients of the underlying mathematical formulation. There has been a key interest in addressing questions of stability and long term response of systems of equations to external inputs and stimuli from the context of equilibrium dynamics. This has led to improvements from mere geometric characterizations to the design of more reliable and robust systems. In this work, we go a little bit further by re-establishing some key aspects of these attempts with the view of arriving at a better understanding of the nontrivial behaviors involving neuronal excitation and response given their applications in various fields.

van der Pol [8], observed that excitable membranes display similarities with self-sustained oscillations. The FHN model is an extension of this overriding concept. The injection of a small current through a cell membrane makes the cell to passively respond with a small voltage deviation. This continues until the attainment of a certain threshold after which the response becomes more prominent. This is called a spike or an action potential and is found mostly in neurons and muscle cells. The display of oscillations and nonlinear effects in living organisms and chaotic motions that may follow thereafter provide a fertile ground for the study of complex dynamical systems that sometimes go beyond the capacity of quantitative analysis alone. The involvement of nonlinear dynamical systems theory in this field, has been found rewarding not only from the viewpoint of providing explanations of this complex and fascinating behavior but also in providing a better understanding of similar activities away from equilibrium. This aspect can be very advantageous in predicting overall outcomes.

The phase space not only analyzes the space determined by all system variables in the space where all possible states of the system exists, but also proffers a rich qualitative information into the behavior of a system under different scenarios. This technique is defined by the so called space volume which for a dynamical system is [9].

\dot{V} = \int_{} \nabla . F (v, w) d w d v = \int (\frac{\partial F_{v}}{\partial v} + \frac{\partial F_{w}}{\partial w}) d w d v = \int (f^{'} (v) - ε) d w d v (1)

where v and w represent cell membrane voltage and recovery respectively.

The above equation yields a geometric representation of the trajectories and vector fields of a dynamical system in the phase plane. In addition an intuitive insight of the behavior of the system is obtained without actually solving the problem.

FitzHugh-nagumo model

Richard FitzHugh carried out an in depth study of the dynamics of the HH model. One of his key observations was that the model followed a fast or slow kinetics. This motivated the conversion of the original four- dimensional HH model to its two-dimensional analog. He initially obtained the following two-variable model

$\begin{array}{l} C \frac{d V}{d t} = I_{e x t} - {\bar{g}}_{N a} m_{\infty}^{3} (0.8 - n) (V - V_{N a}) - {\bar{g}}_{K} n^{4} (V - V_{K}) - g_{L} (V - V_{L}) (2 a) \\ \frac{d n}{d t} = \frac{m_{\infty} (V) - m}{τ_{n}} (2 b) \end{array}$

where V is the potential, t is the time $g_{N A}^{}, g_{k}$ , are the sodium and potassium conductances respectively $I_{e x t}, I_{L}$ , the exit and leakage currents, $g_{L}$ the combined conductance for leakage current, n is channel gating variable, $τ_{n}$ is the voltage sensitive time constant, $m_{\infty} (V)$ is the voltage-sensitive steady state activation function, m is the probability of an activation gate being open. $V_{N a} V_{K}$ are potentials for sodium and potassium ions. $V_{L}$ is ohmic leak potential. For more information concerning the origin of this equation, the reader can refer to Izikevich [10].

He further simplified these two nonlinear coupled and continuous two dimensional ODEs to a general autonomous dynamical system of ODEs that are put in a general form: $\dot{v} = v (t) (α - v (t)) (v (t) - 1) - w (t) + I (2 c)$ $\dot{w} = ε (β v (t) - γ w (t)) (2 d)$

where $ε > 0$ sets the time- scale for the w recovery dynamics, $α$ is an inhibitor threshold-like parameter, the quantity I corresponds to the current applied to the membrane and will represent the bifurcation parameter, the variable $v (t)$ is the voltage-like variable with a cubic nonlinearity that allows for physical regeneration and self-excitation, $w (t)$ is the recovery variable permitting a slower negative feedback and γ is the shunting variable.

Solution of equations 2a and 2b shows a slow collection and fast release of voltage, including nullclines of v and w that have cubic and linear shapes respectively and for certain parameter values, equations (2a) and (2c) become:
$\dot{v} = v - (1 / 3) v^{3} - w + I (3 a)$ $\dot{w} = ϕ (v - 0.8 w + b) (3 b)$

The parameters b and Ø are dimensionless and positive. The amplitude of Ø is the inverse of a time constant and measures how fast the voltage-like and recovery variables $v (t)$ and $w (t)$ change with respect to each other.

The FHN model represents the generalized version of the van der Pol relaxation oscillator. As simple as it looks, it satisfactorily approximates the HH model and retains some of its utilitarian properties like excitability, stable equilibrium and existence can qualitatively reproduce the major dynamical features of realistic neural models. Being dimensionless it may not be particularly relevant from a physiological consideration, yet it offers immense scope for dynamical analysis and paves the way for future study of various forms of excitable systems and biophysical models.

The self- activating and inhibition properties of equations (3a) and (3b) bring to mind the interaction of foxes and rabbits in the Lokta-Volterra equations except that the system is given an inhomogeneous growth by the presence of an applied current which initiates a slow collection and fast release of voltage. Neurons cannot collect voltage and fire immediately, they need time to rest. Moreover, if the stimulus current I does not exceed a certain threshold, then the voltage v does not possess the capacity to excite the system and it remains dormant. On the other hand, when the threshold is exceeded, the dynamics of the system suddenly changes and the neurons start firing with v displaying peak and trough profiles repeatedly. If the stimulus current I continues to increase, another threshold is exceeded and w comes to another steady state that shuts down the activity of v

Dynamical analysis of the FHN model

There is a need to define the general topology that comes with the FHN equations on the v-w plane by focusing on the geometrical depiction of specific orbits, the limit cycles and activities in the vicinity of equilibrium points and in particular.

Valuable information concerning a single point in space as well as details of certain aspects of topological properties can be arrived at by examining qualitatively the phase space spanned by w and v in equations (2a) and (2b). Linearization of the nonlinear system of equations near their fixed or equilibrium points is key to determining the stability property as well as classifying the space paths with the help of the eigenvalues of the Jacobian matrix of the system of equations. In addition, according to the Hartman-Grobman theorem [10], a nonlinear system is geometrically similar (topologically conjugate) to a linear system in a vicinity that is sufficiently close to its equilibrium point.

At steady state, equations (2a) and (2b) can be rewritten to read:

$\dot{v} = f (v, w) = - v (v - a) (v - 1) - w = v^{3} - v^{2} (1 - a) - v a - w = 0 (4 a)$ $\dot{w} = g (v, w) = ε (v - γ w) = b v - c w = 0 (4 b)$

Where $b = ε v, c = ε γ$ respectively. For $v > γ w a n d w > - v (v - a) (v - 1), \dot{w} > 0 a n d \dot{v} < 0$ , and there is a fixed point at the origin. In order to determine the nature of this fixed point, we have to start the linearization procedure by defining the Jacobian at the point, that is $((w, v) \approx (0, 0))$ : ${(\begin{matrix} \frac{\partial f}{\partial v} & \frac{\partial f}{\partial w} \\ \frac{\partial g}{\partial v} & \frac{\partial g}{\partial w} \end{matrix}) |}_{v = 0; w = 0} = (\begin{matrix} - 3 v_{0}^{2} - 2 v_{0} (1 - a) - a & - 1 \\ ε & - ε γ \end{matrix}) = (\begin{matrix} - a & - 1 \\ ε & - ε γ \end{matrix}) = (\begin{matrix} - a & - 1 \\ b & - c \end{matrix}) (5)$

The eigenvalues, in addition to the trace and determinant of equation (5) provide the tools for implementing a comprehensive qualitative analysis of the FHN model provided we have a non-zero determinant and a negative trace. Details of the steps involved can be found in several texts [11-14] . $Δ = \det J = ε (1 + a γ), = (a c + b), τ = t r a c e J = - (a + ε γ) = - (a + b)$ The eigenvalues are represented by: $λ_{1, 2} = \frac{τ \pm \sqrt{(τ^{2} - 4 Δ)}}{2} (6)$

Simulations and Discussions

In order to identify at which levels the FHN model produces excitations as well as make accurate statements concerning the stability of the equilibrium points, we vary the problem parameters and apply different magnitudes of the excitation current as inputs into the FHN model for equations (4a) and (4b). In this work, we make attempts to emphasize the physics of our observations concerning the self-activation and inhibition properties of this model. The time lag that occurs before one of the components of the FHN model reacts may be attributed to the fact that infinitesimal perturbations does not often respond linearly to a system. Though the system may eventually become chaotic as one of the problem parameters is increased, and perturbations grow accordingly yet there could be some points along the trajectory where the effects are negative. To see these effects clearly we have opted to plot the spectral abscissa against the bifurcation parameter for different values of problem parameters as can be seen from figures 3,4,9a and 10. This is accompanied by explanations.

Firstly we adopt the following parameter values: a=0.139, $ε = .008, γ = 2.54, I = 0$ . In Figure (1a) the action potential is indicated by the time evolution of the membrane voltage. Figure 1b compares the temporal evolution of both the membrane potential and the corresponding recovery response for a time interval. A significant variation in the profiles can be observed for $t \approx 140$ after which they both converge asymptotically to zero. Numerical calculations identify an equilibrium point at (0,0) and complex conjugate roots of $(- 0.07966 \pm 0.06659 i)$ which indicates a stable focus at the origin. With respect to equations (4a) and (4b) ,when $v > γ w, \dot{w} > 0,$ and for $w > - v (v - a) (v - 1), \dot{v} < 0$ ; there is clearly a fixed point at the origin as confirmed by Figures (1c) and (1d). Numerical experiments also indicate that most of the interesting dynamical activities happen within a stimulus current range of $0.035010 \leq I \geq 0.035143$ . For example Figure (2a) shows a computer generated v-w time evolution profiles surrounding the rest point, as well as an unstable limit cycle (Figure 2b) and nullclines (Figure 2c) for I=.0351320 with calculated equilibrium point (.0781,.0307). The solution is unstable because the real part of the computed eigenvalue is positive $(.0001 \pm 0871 i)$ . On the other hand a slight change in the value of the bifurcation parameter gives a stable focus as an example for I=.035010, the equilibrium points are (.0778,.0306) accompanied by a complex conjugate eigenvalues of $(- 0.0002 \pm .0871 i)$ . A slight increase of stimulus to I=.0351434 has a computed fixed point of (.0781, 0.0308) and complex conjugate eigenvalues $(0.0002 \pm .0871 i)$ which is indicative of an unstable focus.

To a large extent, the generation of cell activity in a neuron does not depend totally on the stimulus strength alone instead a strong response can be motivated as long as the stimulus exceeds a certain threshold level.

We hasten to comment that sudden dynamical changes of stabilities in nonlinear systems is indicative of Hopf Bifurcation. In between these two cases, there must exist a situation where the real component of the eigenvalues are either negligible or non-existent and the rotational component (the imaginary part) plays a major part in the overall dynamics by yielding an insight into the frequency of the linearized oscillation. Figures 3,4 show the time evolution of eigenvalues and an insight into the bifurcation of the model. The figure is far from uniform for a reason. While the straight or the linear regimes reflect areas where there are single dominant eigenvalues, the nonlinear dome-shaped regions are indicative of regions where the dominant eigenvalues are characterized by complex conjugate pairs.

We further explore the dynamics of the FHN model by turning our attention to equation (3) for the case where all the parameters are much smaller than one; $φ = 0.08 a n d b = 0.7$ . The Jacobian matrix is given as : ${(\begin{matrix} \frac{\partial f}{\partial v} & \frac{\partial f}{\partial w} \\ \frac{\partial g}{\partial v} & \frac{\partial g}{\partial w} \end{matrix}) |}_{V = V_{*}; W = W_{*}} = (\begin{matrix} 1 - V_{*}^{2} & - 1 \\ 0.08 & - .064 \end{matrix}) (7)$

where $V^{*}, W^{*}$ are values of the dependent variables at equilibrium points.

To illustrate this further, we start with I=0. Figure 5a shows the time evolution of the membrane potential and the recovery. It can be seen that the external stimulus factor for the membrane current is not large enough to initiate interaction between the recovery current and the voltage. Figure 5b displays the stable point location on the nullclines. The point of intersection of the nullclines or the equilibrium point is (-1.1994,-0.6243). The eigenvalues are a pair of complex conjugates :(-0.2473+0.2083i, -0.2473-0.2083i) and represents a spiral sink. The magnitude of the imaginary part of these eigenvalues, 0.2083i yields the frequency of the linearized oscillation and corresponds to a period of $\frac{2 π}{0.2083} = 30.16411573$ . Figure 5c depicts a phase portrait of the model for this specific case. The curve in the phase plane where $\dot{v} = 0$ is clearly visible. Away from this curve, $\dot{v}$ is at least an order of magnitude greater than $\dot{w}$ ; in addition, the values of the parameter chosen for this model are an order of magnitude smaller than one. Hence the flow moves rapidly towards the $\dot{v} = 0$ curve. Since there is only one point in the nullclines where the $\dot{v} = 0$ can be crossed by the $\dot{w} = 0$ line, the flow is forced to move slowly in this region and stay very close to the nullcline ( see the cusp very close to the equilibrium point). Most of the exciting topological features happen in this region as $t \to \infty$

Figure 6a shows the time evolution of membrane current and response for a stimulus of I = 0.32 . The equilibrium (v,w) is (-0.9769, -0.3461) and the corresponding complex eigenvalues are (-0.0052+0.2782i, -0.0052-.2782i,) . This again corresponds to a spiral sink. Figure 6b shows the limit cycle for this set of parameters. A more comprehensive picture is displayed in Figure 6c with the limit cycle imposed on the v-w nullclines. However some interesting features come into play when the stimulus current is slightly changed to I = 0.34 (Figure 6d). There is a slight change in equilibrium; (v, w) = (-0.9601,-0.3251) and the complex conjugate eigenvalues come out as (0.0111+0.2748i, 0.0111-0.2748i). The change in the signs of the real components of the complex conjugate pairs of eigenvalues indicate expanding and contracting spirals with respect to the equilibrium points. We note that for the previous case (I=0.32), the real part is negative hence it is a contracting spiral ;the phase paths approach the stable point and the system is stable. For this later case, the real components of the complex conjugate eigenvalues are positive and indicate instability.

The stimulus current that corresponds to the point of transition i.e. where the real component of the complex conjugate eigenvalues attain zero values and passes the imaginary axis in the complex plane was found to be I=0.3264; with a stable point of (-0.9716, -0.3395) and complex conjugate pair of (0000+0.2772i, 0000-0.2772i). The equilibrium point for this case is identified as a center because the spiral component clearly dominates the eigenvalues. As a fact, where one or more of the eigenvalues approach a vanishing real component either stability or instability is possible but non-conclusive. Further increase in the value of the stimulus current, at some stage, say at I = 1.45 results in the decoupling of the voltage potential and the recovery (Figure 6e) with stability coordinates located at (0.9993, 2.1166) . The eigenvalues (-0.0213+0.2807i, -0.0213-0.2807i) are complex conjugates, with a relatively strong rotation component. The corresponding limit cycle is shown in Figure 6f.

For the set of parameters chosen herein, we have seen the topological features displayed by the FHN model as well as the ubiquitous significance of the bifurcation parameter I. If it is below a certain threshold, then the membrane potential v does not come fast enough to excite the system. However when it exceeds this level, the neurons begin to fire; with v exhibiting peaks and troughs accompanied by a relatively slow acting recovery variable or current. As demonstrated, further increase in I propels the system through another threshold, which pushes the recovery current still further into another steady state with an eventual shutting down of the activation potential v. In order to have a clearer understanding of the physics of the dynamics, I is allowed to grow slowly in time by admitting small random noisy perturbations of its magnitude and the overall effect studied.

Figure 7a shows a plot of the recovery variable w versus the input current (bifurcation parameter). As the input current is increased, the system displays a transition from one qualitative type of dynamics to another. This transition involves movements away from an equilibrium point and away from the limit cycle. This is accompanied by an increase in the magnitude negative eigenvalue to a point where it becomes zero with a disappearance of equilibrium. The whole system moves from quiescence to periodic spiking. Or as often as is the case, two complex-conjugate eigenvalues with negative real parts approach the imaginary axis and become totally imaginary. These are hallmarks of Hopf-bifurcation. Right at the beginning, there is a transition from resting to periodic spiking (subcritical) and towards the end, the profile completely switches off. Thus we observe two major qualitative events which are emblematic of Hopf bifurcations. A 3-D representation of the essential features of the prevailing dynamics are depicted in Figure 7b. The two ends of this image reveal that there must be a pair of bifurcations for the specified values of the model parameters as equilibrium is de-stabilized.

Figure 7c shows the time evolution of voltage at fixed points versus the bifurcation parameter. This offers a clue concerning the range of equilibria the system undergoes for that period. Obviously for points before the Hopf bifurcation, the stable fixed point solution will have equal maximum and minimum. This is immediately altered as the values of the stimulus current that achieves bifurcation are attained and the gradient of the profile significantly changed.

Figure 8 shows the limit cycles of FHN model. It clearly indicates the values of the membrane voltage and recovery currents for neuron switch off and on. Further information concerning the dynamics of the Hopf bifurcation can be obtained by exploring information contained in the eigenvalues of the Jacobian matrix. To achieve this, we plot the maximum real parts of the eigenvalues of the Jacobian matrix the so called spectral abscissa versus the stimulus current [15]. The plot, as expected, is not smooth (Figure 9a). It breaks into four distinct regions namely regions where the dominant eigenvalues of the Jacobian matrix are a complex conjugate pair (Hopf bifurcation) as well as those where there is a single dominant eigenvalue. On the whole, they characterize the flow patterns initiated by the range of the stimulus current adopted for the model. And a transition from one qualitative state to another.

Figure 9b is an attempt to complement Figure 9a by giving more information concerning stability or the existence of a limit cycle attractor. Like before we employ various values of the bifurcation parameter I. But for each of the values of I, we neglect the transient period and plot the evolution of minimum and maximum membrane voltages for each value of I. It is observed that when I is small, the solutions converge to the stable equilibrium and both minimum and maximum voltages are equal to the resting voltage. However with an increase in I both values begin to diverge. This signals the existence of a limit cycle attractor whose amplitude increases as I does.

Figure 9c is a plot of the real and imaginary components of the eigenvalues of the FHN model. Basically it displays three regions of interest. The first region is marked by a resting state where the imaginary component is zero and there is little or no rotation. This stage continues until the real component attains a certain negative value and loses its stability effect. There is a dramatic change in profiles for the second region as the imaginary component diverges and takes on positive and negative values. The change in the magnitude of the real component is less pronounced and still remains negative as the system tries to retain its stability. As the real part attains a zero value, the transition from stable to unstable begins to happen and the real component becomes more positive (the second part of region 2 within the elliptical shape). The third region corresponds to a stage where the system becomes more unstable and less rotational. The overall observation is in consonance with the discussion carried out herein.

In Figure 10 The Lyapunov exponents defined by the real part of the Jacobian eigenvalues are plotted against the bifurcation parameter. The plot is divided into four regions which correspond to transitions between regimes of different dynamical activities. Regions marked by straight lines are those that where there is a single dominant real eigenvalue. Whereas the dome shaped region is where the dominant eigenvalues of the Jacobian matrix are a complex conjugate pair.

Conclusion

The FitzHugh-Nagumo (FHN) model is a direct demonstration of a membrane excitation and a simplified version of the Hodgin-Huxley (HH) model. It essentially simulates the generation of action potential on an excitable cell membrane arising from a stimulus current I. In the work presented herein we employ an FHN two first-order ODEs to generate some essential features of neuron messaging and reception as well as the propagation of impulses in multi-cellular organisms. Considerable efforts have been placed in emphasizing the qualitative aspects of the numerical results as a way of better understanding the underlying physics. Interestingly we took advantage of the two-ODE FHN model by using a computational-geometric technique known as the phase plane analysis as an exploratory tool. With this, we were able to compute the steady state values of the solutions (referred to as fixed points) before classifying the dynamics and stability of these solutions. A good many of the issues arising from the ensuing results may have been treated differently elsewhere; nevertheless they still constitute a key area of research in other areas of real-life applications involving the FHN-like models. As mentioned earlier, we have been able to gain a better insight into the underlying physics motivating the self activation-inhibitor response of the FHN model to various input parameters until a threshold value is attained. As can be seen from the graphs of figs. 9 and 10 the values of the numerical abscissa dip below zero for some of those parameter values. This signals a dissipation before a Hopf-Bifurcation is attained.

References

Copyright

© 2021 Onyejekwe OO. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.