Dynamics of the FitzHugh-Nagumo equation with numerical methods

Research Article
Open access

Dynamics of the FitzHugh-Nagumo equation with numerical methods

Erjia Fu 1*
  • 1 Wuhan University    
  • *corresponding author 2020300002101@whu.edu.cn
Published on 17 November 2023 | https://doi.org/10.54254/2753-8818/10/20230339
TNS Vol.10
ISSN (Print): 2753-8818
ISSN (Online): 2753-8826
ISBN (Print): 978-1-83558-131-5
ISBN (Online): 978-1-83558-132-2

Abstract

In this paper, our research focuses on investigating the behavior of the FitzHugh-Nagumo (FHN) equation dynamics, a mathematical model that describes the spiking behavior of neurons. We begin by presenting the mathematical formulation of the FHN equation as a system of two coupled ordinary differential equations. We then apply methods from the theory of dynamical systems to analyze the behavior of the system, including stability analysis and phase-plane analysis. We also discuss the bifurcations that occur in the system as its parameters are varied. In addition, we present numerical methods for solving the FHN equation, including explicit and implicit methods. We compare the accuracy and efficiency of these methods and discuss their suitability for different types of problems.Our results show that the FHN equation exhibits a rich range of dynamical behaviors, including periodic and chaotic solutions, bistability, and hysteresis.Additionally, FHN equation can be used to simulate the behavior of networks of excitable cells, such as neural circuits or cardiac tissue,and FHN models can also be used to study the effects of drugs or other interventions on the behavior of excitable cells. The numerical methods we present provide a powerful tool for studying the FHN equation and its variants, allowing us to explore the parameter space and predict the behavior of the system under different conditions.

Keywords:

FitzHugh-Nagumo equation, dynamical systems, numerical methods,bifurcation analysis

Fu,E. (2023). Dynamics of the FitzHugh-Nagumo equation with numerical methods. Theoretical and Natural Science,10,179-185.
Export citation

1. Introduction

The mathematical model known as the FitzHugh-Nagumo (FHN) equation is utilized to describe the excitability of neurons, which was developed by Richard FitzHugh as a simplified version of Hodgkin and Huxley's model of action potential generation in the giant axon of a squid [1]. Subsequently, Nagumo and his colleagues designed and analyzed an equivalent electric circuit for this model [2].

The study of the firing of neurons has been greatly advanced by the work of Hodgkin and Huxley, who developed a four-dimensional system of differential equations that describe the electrochemical transmission of neuronal signals along the cell membrane, for which they received the Nobel Prize (Hodgkin and Huxley, 1952). This system is similar to those found in electrical circuits, with the neuron consisting of a cell body that receives electrical stimuli, which are then conducted along the axon to other neurons via synapses. However, the current is not made up of electrons but rather ions, primarily sodium and potassium.

The Hodgkin-Huxley system is difficult to work with due to its highly nonlinear nature. R. FitzHugh and Nagumo J [1-2] showed that a breakthrough was achieved by Fitzhugh and Nagumo, who produced a simpler model that captured the essential behavior of nerve impulses, including the phenomenon of excitability. The FitzHugh-Nagumo model consists of a set of two coupled nonlinear differential equations [3]:

v' = v - \( \frac{{v^{3}}}{3} \) - w + I

(1)

w' = \( \frac{ε(v + a - bw)}{τ} \)

(2)

Here, the equation involves several variables, including v that indicates the neuron's membrane potential, w that tracks the refractory period through a recovery variable, I that represents the input current, and ε, τ, and a that are constants. The FitzHugh-Nagumo model has been used extensively to study a variety of phenomena, including pattern formation, synchronization, and bifurcation in neuronal networks. M. Jinyan ZhangAnd[4] introduces some skills to illustrate some background information about Hopf bifurcation as follows:

Consider the Hopf bifurcation:

\( \begin{cases} \begin{array}{c} \dot{x}=y+λ(x-x³+I) \\ \dot{y}=-x+λ(a-by) \end{array} \end{cases} \)

(3)

Denote function:

\( Φ(A)=\int _{0}^{2π}x(x-\frac{{x^{2}}}{3}+I)+y(a-by) \)

(4)

Let \( x=Acos{t} \) , \( y=Asin{t} \) ,

\( Φ(A)=\int _{0}^{2π}Acos{t}(Acos{t}-A\frac{{A^{3}}{cos^{3}}{t}}{3}+I)+ Asin{t}(a-bAsin{t})dt \)

\( =\int _{0}^{2π}[{A^{2}}{cos^{2}}{t}-\frac{{A^{4}}{cos^{4}}{T}}{3}+IAcos{t}+Asin{t}a-b{A^{2}}{sin^{2}}{t}]dt \)

\( \int _{0}^{2π}[{A^{2}}(\frac{cos{2}t+1}{2})-{A^{4}}(\frac{π}{8}+\frac{1}{12}sin{2}t+\frac{1}{96}sin{4}t)-b{A^{2}}(\frac{1-cos{2}t}{2})]dt \)

= \( {A^{2}}π-\frac{{A^{4}}}{4}π-b{A^{2}}π \)

Let \( Φ(A)=0 \) , A=0 or \( A=±2\sqrt[]{1-b} \)

\( {Φ^{ \prime }}(A)=2Aπ-{A^{3}}π-2bAπ \)

\( {Φ^{ \prime }}(-2\sqrt[]{-b})=-4\sqrt[]{1-b}π-8{(1-b)^{\frac{3}{2}}}+4b\sqrt[]{1-b}π \lt 0 \)

So, we can derive from the above conclusion that the Hopf

bifurcation \( \begin{cases} \begin{array}{c} \dot{x}=y+λ(x-\frac{{x^{3}}}{3}+I) \\ \dot{y}=-x+λ(a-by) \end{array} \end{cases} \) have limit cycle near \( \begin{cases} \begin{array}{c} x=2\sqrt[]{1-b}cos{t} \\ y=2\sqrt[]{1-b}sin{t} \end{array} \end{cases} \)

when \( λ \gt 0 \) it is stable, when \( λ \lt 0 \) it is unstable.

Our investigation focuses on examining the behavior of the FHN equation dynamics using the framework of ordinary differential dynamical systems. We begin by presenting the mathematical formulation of the FHN equation and its basic properties. We then apply methods from the theory of dynamical systems to analyze the behavior of the system, including stability analysis, phase-plane analysis, and bifurcation analysis.Meanwhile,we present some numerical methods for solving the FHN equation, including explicit and implicit time-stepping methods and shooting methods for boundary value problems. We compare the accuracy and efficiency of these methods and discuss their suitability for different types of problems.Numerical methods for solving the FitzHugh-Nagumo equation include explicit and implicit methods, as well as numerical schemes based on finite difference, finite element, and spectral methods. These methods can be used to obtain accurate solutions to the equation and to study the behavior of its solutions under different initial conditions and parameter values.

2. Background

2.1. Introduction to some basic numerical methods and its appliance in FitzHugh-Nagumo equation

The explicit Euler method is a first-order method for numerically approximating the solution of an ordinary differential equation (ODE) at discrete time steps. The first-order ordinary differential equation can be solved using the explicit Euler method, which has a general form:

\( {y_{n+1}}={y_{n}}+hf({t_{n}},{y_{n}}) \)

(5)

where \( {y_{n}} \) is the numerical approximation of the solution at time \( {t_{n}} \) , h is the time step, and \( f({t_{n}},{y_{n}}) \) is the right-hand side of the ODE evaluated at time \( {t_{n}} \) and the numerical approximation \( {y_{n}} \) .

Using the explicit Euler method to the FitzHugh-Nagumo Equation, we can update the numerical approximations \( {v_{n}} \) and \( {w_{n}} \) at each time step as follows:

\( {v_{n+1}}={v_{n}}+h({v_{n}}-(v_{n}^{3})/3-{w_{n}}+I) \)

(6)

\( {w_{n+1}}={w_{n}}+h(({v_{n}}+a-b{w_{n}})/tau) \)

(7)

where h is the time step.

It's worth noting that explicit time-stepping methods like the explicit Euler method may not be stable for certain values of the time step h and system parameters. In such cases, implicit time-stepping methods may be more suitable.

The implicit Euler method is a first-order method for numerically approximating the solution of an ODE at discrete time steps. The general form of the implicit Euler method for a first-order ODE is:

\( {y_{n+1}}={y_{n}}+hf({t_{n+1}},{y_{n+1}}) \)

(8)

Here, \( {y_{n}} \) represents the numerical estimation of the solution at time \( {t_{n}} \) , h denotes the time interval, and \( f({t_{n+1}},{y_{n+1}}) \) corresponds to the right-hand side of the ordinary differential equation computed at time \( {t_{n+1}} \) and the numerical approximation \( {y_{n+1}} \) .

Using the implicit Euler method to the FitzHugh-Nagumo Equation,we can update the numerical approximations v_(n+1)and w_(n+1) at each time step as follows:

\( {v_{n+1}}={v_{n}}+h({v_{n+1}}-(v_{n+1}^{3})/3-{w_{n+1}}+I) \)

(9)

\( {w_{n+1}}={w_{n}}+h(({v_{n+1}}+a-b{w_{n+1}})/tau) \)

(10)

where h is the time step.

Note that the implicit Euler method requires the numerical approximation at the next time step, \( {y_{n+1}} \) , to be solved using a nonlinear equation that involves \( {y_{n+1}} \) on both sides of the equation. This nonlinear equation can be solved iteratively using methods such as Newton's method or the fixed-point iteration method.

Implicit time-stepping methods like the implicit Euler method are generally more stable for stiff systems like the FitzHugh-Nagumo Equation, but they may require more computation time per time step.

The shooting method is a numerical technique for solving boundary value problems (BVPs) that involves iteratively guessing the boundary conditions and adjusting them until the solution satisfies the differential equation.

The shooting method proceeds as follows:

Choose an initial guess for the value of one of the boundary conditions, say V ( \( {t_{0}} \) )= \( {v_{0}}. \)

Integrate the FitzHugh-Nagumo equations numerically from t0 to \( {t_{f}} \) using a suitable ODE solver, such as the fourth-order Runge-Kutta method.

Check whether the other boundary condition is satisfied at \( {t_{f}} \) . If it is, we have found a solution to the BVP. If not, adjust the initial guess for v0 and repeat steps 2 and 3 until the desired accuracy is achieved.

The shooting method can be further refined by using a root-finding algorithm, such as Newton's method, to improve the accuracy of the solution.

Note that the choice of initial guess for \( {v_{0}} \) can have a significant impact on the convergence of the shooting method. It is often helpful to use a combination of analytical and numerical techniques to determine a suitable initial guess.

2.2. Some advanced and modern methods to solve FitzHugh-Nagumo equation

In addition to the basic numerical methods mentioned in the previous section, there are several advanced and modern methods for solving the FitzHugh-Nagumo equation.

Spectral methods: Spectral methods are a class of numerical methods that use orthogonal functions, such as Fourier or Chebyshev polynomials, to approximate the solution of a differential equation. Spectral methods have high accuracy and can be used to solve nonlinear problems, making them well-suited for the FitzHugh-Nagumo equation.

Finite element methods: Finite element methods (FEM) are a popular class of numerical methods for solving partial differential equations (PDEs). FEM involves dividing the domain of the PDE into smaller subdomains and approximating the solution using a set of basic functions defined over each subdomain. FEM can be used to solve the FitzHugh-Nagumo equation with complex geometries or boundary conditions.

Neural network methods: Neural network methods, such as the artificial neural network (ANN) and deep learning models, have recently been applied to solve differential equations. These methods use the flexibility and parallelism of neural networks to approximate the solution of a differential equation. Neural network methods have shown promising results in solving the FitzHugh-Nagumo equation and other complex systems.

Reduced-order modeling: Reduced-order modeling (ROM) is a technique used to approximate the solution of a high-dimensional system using a low-dimensional representation. ROM can be used to solve the FitzHugh-Nagumo equation by approximating the solution using a reduced set of basic functions. ROM has the advantage of reducing computational cost while maintaining accuracy.

Hybrid methods: Hybrid methods combine different numerical methods to solve a differential equation. For example, a hybrid method may use the spectral method to approximate the solution in the interior of the domain and use the shooting method to satisfy the boundary conditions. Hybrid methods can be used to solve the FitzHugh-Nagumo equation with high accuracy and efficiency.

3. Recent researches

Hoff A, Celestino A and Medeiros E S [5-7] present a numerical study of the dynamical behavior of two coupled FitzHugh-Nagumo (FHN) oscillators and explores the system's response to different parameter values, including both single and double coupling. The FHN oscillator model mathematically describes neural impulse conduction within neurons, encompassing the dynamics of cell membrane potential, the excitable membrane's resting bias, and the threshold for action potential.

The authors use numerical methods to construct Lyapunov and isochronous charts, as well as bifurcation diagrams for the two types of coupling, to calculate the Lyapunov exponents and the number of local maxima of variables within one cycle. The results show the existence of Arnold tongues in the network of two single-coupled oscillators, which self-organize sequentially in a branch of the Stern-Brocot tree. The bifurcation diagrams demonstrate the relationship between these Arnold tongues and other periodic structures in the Lyapunov chart. Additionally, the system exhibits multistability, which is demonstrated on the plane of attraction basins.

The paper also presents the network's dynamic characteristics as a function of the coupling strength and other system parameters. Through numerical solution using a fixed time step and a 500×500 value grid via the Runge-Kutta method, the authors find that the system displays numerous steady states and exhibits Arnold tongues in a network of two unidirectional coupled oscillators. The dynamic characteristics of the network are described by analyzing the changes in coupling strength and other system parameters.

Overall, this study contributes to the understanding of the dynamics of coupled FHN oscillators and the effects of different parameters on their behavior. The findings have implications for the study of neural systems and could inform the design of new models for simulating neural activity. Future research could expand on this work by investigating the behavior of larger networks of FHN oscillators and exploring the effects of different coupling topologies.

Liu F and Yu Q [8-9] show the use of a fractional FitzHugh-Nagumo monodomain model, which employs fractional derivatives to capture spatial heterogeneities and connectivities in the extracellular domain. This approach leads to improved simulation accuracy of cardiac tissue heterogeneity. Describing the membrane potential function and its coupling with the system, the model is composed of a set of fractional Riesz space nonlinear reaction-diffusion models and a discrete differential equation, and it extends the standard monodomain model to describe the propagation of potential in heterogeneous cardiac tissue. Investigating the impact of coupling on the system's behavior, the authors also present a numerical technique for solving the fractional coupled monodomain model. The model is represented as a two-dimensional fractional Riesz space nonlinear reaction-diffusion model subject to initial and boundary conditions, and decoupled numerical techniques are employed for its solution.To solve the model's nonlinear reaction-diffusion equations, the authors use an operator splitting method that is applicable to locally Lipschitz continuous functions. The numerical simulation of the equations is carried out in a finite domain using the shifted Grünwald-Letnikov scheme for discretization. Subsequently, the authors prove the stability and convergence of an implicit numerical method for a continuous problem, demonstrating two theorems. The first theorem proves the stability of the implicit numerical method, which is unconditionally stable and limits the error of each iteration to a constant C, independent of the step sizes hx, hy, and τ. The proof demonstrates that the error of each iteration is limited by a factor of C times the error of the previous iteration. The second theorem proves the convergence of the implicit numerical method, with the error between the solution of the implicit numerical method and the true solution of the continuous problem limited by a constant C multiplied by the step sizes hx, hy, and τ, assuming a smooth solution to the continuous problem. The proof uses the concept of error and the fact that the error of each iteration is limited by a constant C. The proof also uses a technique called reverse error analysis to show that the error between the solution of the implicit numerical method and the true solution of the continuous problem is limited by a constant C multiplied by the step sizes hx, hy, and τ. Subsequently, the authors present numerical results to verify the correctness of the theoretical analysis. Finally, the authors discuss the rationale for action potential duration (APD) as an important biomarker of arrhythmia risk in cardiac electrophysiology.

Bonaventura L [10-12] demonstrate the importance of synchronization in neuronal networks and the challenges faced in numerical simulations. The dynamics of neuronal networks are complex, exhibiting nonlinearity and multiple time scales, and the synchronicity properties vary with the size and structure of the network. Efficient simulation of neuronal networks necessitates specialized techniques that can handle stiff problems that arise from the slow-fast nature of the system's dynamics. The article proposes an economic version of the standard implicit ODE solver specifically designed for efficient and accurate simulation of neuronal networks, simplifying the algebraic system to be solved during each time interval. The developed method is tested on networks based on three different single-cell models and on sparse, medium, and dense networks to demonstrate its efficiency. The article concludes by discussing future research directions in the field.The author describes three different single-neuron models used to construct networks. The first model is the classic FitzHugh-Nagumo model which incorporates two variables that respectively represent the membrane potential and slow recovery variable. The second model is an extension of the FHN model, which includes an additional variable that represents the neuron ICC. This model also involves the incorporation of sigmoid activation function parameters for the ICC variable. The third model is the Hodgkin-Huxley (HH) model consisting of three variables representing the membrane potential, the ionic transport across the membrane through ion channels, and an external current parameter. The equations for each model are provided, and the meanings and ranges of their parameters are explained.Then, author describes the construction of two different neuronal network models: the FitzHugh-Nagumo neuron network and the intracellular calcium concentration neuron network. The FitzHugh-Nagumo neuron network is a model representing a group of electrically coupled neurons connected by symmetric connections. The model has two variables, x and y, defined by functions f and g. The equation for x describes the membrane potential of each neuron, while the equation for y describes the recovery variable of the neuron. The model equation can be represented as a vector and is given by the matrix A multiplied by the state vector, where A is a matrix depending on the coupling strength between neurons. The intracellular calcium concentration neuron network is a model representing a group of electrically coupled neurons connected by symmetric connections. The model has three variables, x, y, and z, representing the membrane potential, the recovery variable, and the intracellular calcium ion concentration of each neuron, respectively. The equations for x and y are similar to those of the FitzHugh-Nagumo model, while the equation for z describes the dynamics of the calcium ion concentration within the neuron. The model equation can also be represented as a vector and depends on the coupling strength between neurons.Meanwhile,the author describes a method for solving nonlinear systems of equations that arise in numerical differential equation simulations. Specifically, it describes a method that uses implicit discretization to solve these equations, which can be more accurate than explicit discretization but also more computationally expensive. The method involves using the implicit Euler method as a starting point and then extending it to more accurate multistep or multistage implicit methods. The crucial computational step involves reformulating the solution of the nonlinear system into the solution of a linear system, where the size of the linear system is equivalent to the size of a vector unknown, rather than its multiple. The method is demonstrated through two examples, one involving the FitzHugh-Nagumo neuron network and the other involving the intracellular calcium concentration neuron network. In both cases, the method involves using the Newton method to solve the nonlinear system, which requires solving a specific form of linear system. As a result of the unique structure of the linear system, it can be solved effectively by direct methods, which leads to the computational efficiency of the method.To conclusion, the author introduces a high-order implicit method called ESDIRK that can solve nonlinear systems of equations in numerical differential equation simulations. It can also handle stiff problems, making it suitable for simulating neuronal networks.

4. Challenge

One area for future research is the application of the FitzHugh-Nagumo equation to more complex systems. For example, while the equation is often used to model single neurons, it can also be applied to networks of neurons. Future research could focus on studying the behavior of these networks and the emergence of complex dynamics such as synchronization, which is important in understanding brain function.Moving forward, there is potential for the proposed approach to be applied to a broader range of large-scale neural networks with biological significance. Additionally, further advancements in computational efficiency could be achieved by exploring self-adjusting multirate extensions of these numerical methods.

Another area for future research is the development of new numerical methods for solving the FitzHugh-Nagumo equation. While there are many methods available, each has its own strengths and weaknesses, and there is still a need for more efficient and accurate methods. In particular, methods that can handle stochasticity and uncertainty in the system could be valuable.

In addition, future research could also explore the relationship between the FitzHugh-Nagumo equation and other models in mathematical neuroscience. For example, the FitzHugh-Nagumo equation is related to the Hodgkin-Huxley model, which describes the dynamics of ion channels in neurons. Comparing and contrasting these models could provide insights into the different levels of abstraction used in modeling neuronal behavior.

And it worth to mention that in studying the FitzHugh-Nagumo equation include the need to incorporate realistic biological parameters into the model, as well as the challenge of validating model predictions experimentally. In addition, the equation is often simplified from more complex biological systems, which can limit its applicability in certain contexts.

Overall, the FitzHugh-Nagumo equation continues to be an important model in mathematical neuroscience, with many opportunities for future research and challenges to overcome. Advances in numerical methods and the application of the equation to more complex systems could lead to a better understanding of brain function and new insights into the treatment of neurological disorders.

5. Conclusion

In conclusion, the FitzHugh-Nagumo equation is an important model in mathematical neuroscience that is widely used for modeling single neurons and neural networks. The equation's strengths lie in its simplicity and solvability, but it also faces challenges and limitations, such as the realism of model parameters and experimental validation of model predictions. Future research directions include exploring the relationship between the FitzHugh-Nagumo equation and other mathematical neuroscience models, developing new numerical methods to improve the accuracy and efficiency of the model, and applying the model to more complex systems. With deeper research on the FitzHugh-Nagumo equation, we can better understand the functioning of the nervous system and provide new ideas and methods for the treatment of neurological disorders.


References

[1]. R. FitzHugh. Impulses and physiological states in theoretical models of nerve membrane. Biophysical Journal, 1(6):445–466, 1961.

[2]. Nagumo J, Arimoto S, Yoshizawa S. An active pulse transmission line simulating nerve axon[J]. Proceedings of the IRE, 1962, 50(10): 2061-2070.

[3]. M. Morris W. Hirsch, Stephen Smale, Robert L. Devaney. Differential Equations, Dynamical Systems, and an Introduction to Chaos, [M]. Oxford, Elsevier, 2013

[4]. M. Jinyan Zhang, Beiye Feng, The Stability Theory of Ordinary Differential Equations and Bifurcation Problems, [M]. Beijing, Beijing University Press, 2008.

[5]. Hoff A, dos Santos J V, Manchein C, et al. Numerical bifurcation analysis of two coupled FitzHugh-Nagumo oscillators[J]. The European Physical Journal B, 2014, 87: 1-9.

[6]. Celestino A, Manchein C, Albuquerque H A, et al. Stable structures in parameter space and optimal ratchet transport[J]. Communications in Nonlinear Science and Numerical Simulation, 2014, 19(1): 139-149.

[7]. Medeiros E S, Medrano-T R O, Caldas I L, et al. Torsion-adding and asymptotic winding number for periodic window sequences[J]. Physics Letters A, 2013, 377(8): 628-631.

[8]. Liu F, Turner I, Anh V, et al. A numerical method for the fractional Fitzhugh–Nagumo monodomain model[J]. Anziam Journal, 2012, 54: C608-C629.

[9]. Yu Q, Liu F, Turner I, et al. Stability and convergence of an implicit numerical method for the space and time fractional Bloch–Torrey equation[J]. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2013, 371(1990): 20120150

[10]. Bonaventura L, Fernández-García S, Gómez-Mármol M. Efficient implicit solvers for models of neuronal networks[J]. arXiv preprint arXiv:2210.01697, 2022.

[11]. Bandera A, Fernández-García S, Gómez-Mármol M, et al. A multiple timescale network model of intracellular calcium concentrations in coupled neurons: Insights from ROM simulations[J]. Mathematical Modelling of Natural Phenomena, 2022, 17: 11.

[12]. Bonaventura L, Mármol M G. The TR-BDF2 method for second order problems in structural mechanics[J]. Computers & Mathematics with Applications, 2021, 92: 13-26.


Cite this article

Fu,E. (2023). Dynamics of the FitzHugh-Nagumo equation with numerical methods. Theoretical and Natural Science,10,179-185.

Data availability

The datasets used and/or analyzed during the current study will be available from the authors upon reasonable request.

Disclaimer/Publisher's Note

The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of EWA Publishing and/or the editor(s). EWA Publishing and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

About volume

Volume title: Proceedings of the 2023 International Conference on Mathematical Physics and Computational Simulation

ISBN:978-1-83558-131-5(Print) / 978-1-83558-132-2(Online)
Editor:Roman Bauer
Conference website: https://www.confmpcs.org/
Conference date: 12 August 2023
Series: Theoretical and Natural Science
Volume number: Vol.10
ISSN:2753-8818(Print) / 2753-8826(Online)

© 2024 by the author(s). Licensee EWA Publishing, Oxford, UK. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license. Authors who publish this series agree to the following terms:
1. Authors retain copyright and grant the series right of first publication with the work simultaneously licensed under a Creative Commons Attribution License that allows others to share the work with an acknowledgment of the work's authorship and initial publication in this series.
2. Authors are able to enter into separate, additional contractual arrangements for the non-exclusive distribution of the series's published version of the work (e.g., post it to an institutional repository or publish it in a book), with an acknowledgment of its initial publication in this series.
3. Authors are permitted and encouraged to post their work online (e.g., in institutional repositories or on their website) prior to and during the submission process, as it can lead to productive exchanges, as well as earlier and greater citation of published work (See Open access policy for details).

References

[1]. R. FitzHugh. Impulses and physiological states in theoretical models of nerve membrane. Biophysical Journal, 1(6):445–466, 1961.

[2]. Nagumo J, Arimoto S, Yoshizawa S. An active pulse transmission line simulating nerve axon[J]. Proceedings of the IRE, 1962, 50(10): 2061-2070.

[3]. M. Morris W. Hirsch, Stephen Smale, Robert L. Devaney. Differential Equations, Dynamical Systems, and an Introduction to Chaos, [M]. Oxford, Elsevier, 2013

[4]. M. Jinyan Zhang, Beiye Feng, The Stability Theory of Ordinary Differential Equations and Bifurcation Problems, [M]. Beijing, Beijing University Press, 2008.

[5]. Hoff A, dos Santos J V, Manchein C, et al. Numerical bifurcation analysis of two coupled FitzHugh-Nagumo oscillators[J]. The European Physical Journal B, 2014, 87: 1-9.

[6]. Celestino A, Manchein C, Albuquerque H A, et al. Stable structures in parameter space and optimal ratchet transport[J]. Communications in Nonlinear Science and Numerical Simulation, 2014, 19(1): 139-149.

[7]. Medeiros E S, Medrano-T R O, Caldas I L, et al. Torsion-adding and asymptotic winding number for periodic window sequences[J]. Physics Letters A, 2013, 377(8): 628-631.

[8]. Liu F, Turner I, Anh V, et al. A numerical method for the fractional Fitzhugh–Nagumo monodomain model[J]. Anziam Journal, 2012, 54: C608-C629.

[9]. Yu Q, Liu F, Turner I, et al. Stability and convergence of an implicit numerical method for the space and time fractional Bloch–Torrey equation[J]. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2013, 371(1990): 20120150

[10]. Bonaventura L, Fernández-García S, Gómez-Mármol M. Efficient implicit solvers for models of neuronal networks[J]. arXiv preprint arXiv:2210.01697, 2022.

[11]. Bandera A, Fernández-García S, Gómez-Mármol M, et al. A multiple timescale network model of intracellular calcium concentrations in coupled neurons: Insights from ROM simulations[J]. Mathematical Modelling of Natural Phenomena, 2022, 17: 11.

[12]. Bonaventura L, Mármol M G. The TR-BDF2 method for second order problems in structural mechanics[J]. Computers & Mathematics with Applications, 2021, 92: 13-26.