Open AccessArticle

Nonstandard Nearly Exact Analysis of the FitzHugh–Nagumo Model

Shahid

^1,2

Mujahid Abbas

and

Eddy Kwessi

^3,*

Department of Mathematics, University of Karachi, Karachi 75270, Pakistan

Abdus Salam School of Mathematical Sciences, Government College University, Lahore 54000, Pakistan

Department of Mathematics, TrinityUniversity, 1 Trinity Place, San Antonio, TX 78212, USA

Author to whom correspondence should be addressed.

Symmetry 2024, 16(5), 585; https://doi.org/10.3390/sym16050585

Submission received: 12 February 2024 / Revised: 21 March 2024 / Accepted: 15 April 2024 / Published: 9 May 2024

(This article belongs to the Special Issue Nonlinear Symmetric Systems and Chaotic Systems in Engineering)

Download

Browse Figures

Versions Notes

Abstract

The FitzHugh–Nagumo model has been used empirically to model certain types of neuronal activities. It is also a non-linear dynamical system applicable to chemical kinetics, population dynamics, epidemiology and pattern formation. In the literature, many approaches have been proposed to study its dynamics. In this paper, initially, we have employed cutting-edge tools from discrete dynamics for discretization and fixed points. It has been proven that an exact discrete scheme exists for this paradigm. This project also considers the phase space and integral surfaces of these evolutionary equations. In addition, it carries out a thorough symmetry analysis of this reaction diffusion system to find equivalent systems. Moreover, steady-state solutions are obtained using ansatzes for traveling wave solutions. The existence of infinite traveling wave solutions has also been proven. Yet again, this investigation establishes the potential of symmetry methods to unravel non-linearity. Finally, singular perturbation theory has been employed to obtain analytical approximations and to study stability in different parameter regimes.

Keywords:

fix points; symmetry analysis; traveling wave solutions; reaction diffusion system; discrete dynamics; steady states

1. Introduction

Embedded within the cells of nervous system exists a composition of a special classification termed neurons. These types of neurons are a highly significant aspect acting as the main worker units in the nervous system. Neurons are designed in a very unique way in the communication system of gland cells and muscles, thus enabling them, importantly, to act as carriers of information between them. It is very natural to think that any model for neurons would be complicated. The FitzHugh–Nagumo model arose as a special case of the famed Hodgkin–Huxley (Nobel Prize 1963) model. The pioneering work by Nagumo [1] opened several avenues for research in neuroscience. As a result, there exist many different variants of the FitzHugh–Nagumo model. It can be used to model how the voltage of neurons changes as a function of time and space:

\{\begin{matrix} \frac{\partial u}{\partial t} & = D \frac{\partial^{2} u}{\partial x^{2}} - v + g (u) \\ \frac{\partial v}{\partial t} & = ε [- β v + c + u] \end{matrix},

(1)

where

g (u) = u - \frac{u^{3}}{3}

and

D, ε, β

and c are real parameters. The holy grail is the integral surfaces

u (x, t)

and

v (x, t)

. These quantities are known as the voltage potential (fast variable) and recovery variable (slow variable). This is a well-established neurodynamic model. A cursory Google search produces more than 200,000 results for FitzHugh–Nagumo. Despite its importance, the variant under consideration has not been solved exactly. Although, if we change

g (u)

slightly, introduce a diffusion term with a Laplacian in the second equation, or impose some extra conditions on the model, then it becomes very easy to work out numerous exact solutions. In this article, we focus on (1). In the natural and social sciences, the atomistic view of systems is deemed inadequate. This research paper integrates comprehensively discrete and continuous, quantitative and qualitative, automated and manual, numerical and exact, and geometrical and analytical methodologies towards this study to give vitality and relevance to the FitzHugh–Nagumo paradigm. The research study’s aim is underscored in recent mathematical discoveries in 2021, which were applied and implemented towards a better understanding of the FitzHugh–Nagumo model. In a nutshell, this is an incisive analytical examination of the FitzHugh–Nagumo model using contemporary mathematics to revitalize its significance in the 21st century.

Solutions of ordinary differential equations are integral curves, but at times, it is enough to obtain a solution of their corresponding difference equations. Similarly, although the solution of a pde could be a surface, at times, it is enough to obtain a lattice that approximates an integral surface. This lattice is the solution of a difference scheme corresponding to a given pde (Figure 1).

In the context of discrete dynamics, Mickens [2] revolutionized scientific computing by introducing nonstandard finite difference schemes (NSFDs). Nonstandard finite difference schemes have led to another structural-preserving modification proposed by Eddy [3] in their theoretical framework termed the nearly exact discretization scheme (NEDS). Using this scheme, which preserves the dynamical properties of the continuous model, Eddy [4] studied the FitzHugh–Nagumo model. This current paper is an extension of [4], as it considers the diffusion term missing in that paper. The collective resources available on the FitzHugh–Nagumo model without diffusion significantly outnumber those on the FitzHugh–Nagumo model with diffusion. Most non-linear dynamical systems are demanding. Unlike linear systems, in the absence of generalized methods, each non-linear system demands special attention. In discretization, utmost care is also needed. The scheme must respect the bifurcation, stability and equilibria of the continuous model. This problem is more striking for the FitzHugh–Nagumo model with diffusion. With an additional independent spatial variable, fixed points show a completely different behavior. In the beginning, when steady states (10) were computed for discretized models, they showed dependence on the spatial variable. This was later confirmed by an exact analysis ((24), (26) and (27)) of the continuous model. These fixed points on a continuum make it difficult to perform stability analyses. With the inclusion of the diffusion term in the model, another phenomenon occurs, as now, the parametric space is bigger, so the singular points of (1) behave differently; that is to say, we now have different bifurcations. Luckily, exact isolated singular points have also been computed. So, at least in the immediate vicinity of these isolated equilibrium points, the mathematics of neurons can be accurately observed. Never before have these fixed points been computed. Therefore, bifurcation and Turing patterns can now be studied exactly. With four natural parameters in the model, chaos can now be monitored exactly with these knobs adjustments.

This paper is organized as follows: Section 2 explores the latest development of tools in the field of discrete dynamical systems. In particular, it addresses the theory of nonstandard finite difference schemes. An intuitive argument for the existence of exact schemes has been presented. For a huge number of coupled systems of first-order equations, it has been shown that there will be an exact scheme. After that, we recall key ideas from the theory of second-order linear difference equations. By using a nearly exact discretization scheme and the results from the theory of difference equations, discrete steady states of the discretized FitzHugh–Nagumo model with diffusion have been determined. Furthermore, we visualize these fixed points in space. Lastly, we identify the limitation of perturbation method to carry out a stability analysis. Section 3 makes the transition from discrete to continuous neural dynamics. It also examines the phase space and integral surfaces of this well-established model from neurodynamics. There is limited research available on the similarity analysis of the FitzHugh–Nagumo model with diffusion, although some articles focus on the conditional symmetry [5] of the fast equation. The slow variable is always neglected. Section 3 also carries out a rigorous symmetry analysis of this reaction diffusion system. Up-to-date tools in similarity analyses have been used to attack the problem. Manipulating different ansatzes, traveling wave solutions are sought. Seeing a limited success with built-in utilities in Maple, Matlab and Mathematica, new algorithms for tanh, exp and Ricatti methods for traveling wave solutions have been generated. Incorporation of the latest modifications proposed by scholars working with these ansatzes has been attempted. Section 4 employs singular perturbation theory to obtain analytical approximations. It is extremely difficult to perform an asymptotic analysis of a system of partial differential equations. An attempt is made to combine lie symmetry analysis with asymptotics. In this last section, it has been demonstrated that infinite traveling wave solutions and an exact discrete scheme exist for (1). A new analytical geometric perspective of traveling waves is given. Even with the spatial variable present in the FitzHugh–Nagumo model with diffusion, the

u v

phase plane has been projected out to study the dynamical system.

2. Discrete Dynamics

The first step in this analysis of (1) involved discretizing the system. Since the available research includes some work on discretization, such as in [6], an innovative technique has been attempted. Pott [7] identifies the subtleties of the numerator and the denominator in the standard definition of a derivative and generalizes (2) to find that: (i) the rate at which the numerator and the denominator terms reach the base point need not be uniform and (ii) any increment or decrement therein need not be strictly linear initially.

\frac{d u}{d t} = {lim}_{h \to 0} \frac{u [t + ϕ_{1} (h)] - u (t)}{ϕ_{2} (h)}

(2)

Mickens [2] capitalizes on Potts’ generalization to propose a radical change in the numerical analysis of differential equations and presents an argument to prove the existence of an exact discretization scheme. The analysis in the current work demonstrates that there will always be an exact discretization scheme to find solution curves. For a wide number of coupled systems of first-order ordinary differential equations, the existence and uniqueness of solution curves have been well established. In principle, integral curves form a congruence, and they constitute orbits. The degree of freedom on these solution curves is one. Therefore, considering this as a time or a step size parameter, it is demonstrable that first-order ordinary differential equation will always have an exact scheme. However, while Mickens [2] assumed that this scheme is unique and the denominator function is always nontrivial, there may be examples of a non-unique exact scheme and a trivial denominator as well for some problems.

Sketching a simple geometric proof for the existence of an exact scheme is possible. The first derivative signifies the slope of the tangent line to the solution curves, which are smooth enough. The slope of this tangent line can be computed by the quotient of the rise over run. However, as these curves are one-dimensional manifolds, it is possible to calculate the exact slope by controlling the denominator only (3).

\frac{\partial u (x_{i}, t_{n})}{\partial t} \approx \frac{U_{i}^{(n + 1)} - U_{i}^{(n)}}{ϕ_{t}}

(3)

This computation can be carried out continuously along the entire solution curve. However, the difficulty is in finding the denominator function

ϕ_{t}

that makes this approximation an equality. The solution to this exact scheme, which is a string of discrete numbers, lies precisely on the solution curves of the differential equation. Further, this exact scheme works independently of the step size.

2.1. Review of Second-Order Linear Difference Equations

First, consider a second-order non-homogeneous linear difference equation:

X_{i + 1} + α_{i} X_{i} + β_{i} X_{i - 1} = γ_{i - 1}, for i = 1, 2, \dots .

(4)

where

α = \{α_{i}\}, β = \{β_{i}\}

, and

γ = \{γ_{i}\}

are sequences of real numbers with

γ \neq 0

. Also, consider the corresponding homogeneous linear difference equation:

X_{i + 1} + α_{i} X_{i} + β_{i} X_{i - 1} = 0, for i = 1, 2, \dots .

(5)

The following definition is an analogue of the Wronskian for differential equations.

Definition 1.

Let

x = \{x_{i}, i \in N\}

and

y = \{y_{i}, i \in N\}

be two solutions of the homogeneous linear difference Equation (5). The Casoratian of x and y, denoted by

C (x, y) = {C (x, y; i}, i \in N}

, is given by

C (x, y,; i) = |\begin{matrix} x_{i - 1} & y_{i - 1} \\ x_{i} & y_{i} \end{matrix}| = x_{i - 1} y_{i} - y_{i - 1} x_{i} for i = 1, 2, \dots .

Now, recall Heymann’s theorem on the relationship between the Casoratian and the linear independence of solutions to the homogeneous difference equation.

Theorem 1

(Heymann’s Theorem [8]). Let

x = \{x_{i}, i \in N\}

and

y = \{y_{i}, i \in N\}

be two solutions of the homogeneous linear difference Equation (5).

1.: The Casoratian $\{C_{i}, i \in N\}$ satisfies the linear homogeneous difference equation

$C_{i} - γ_{i - 1} C_{i - 1} = 0, for i = 1, 2, \dots .$
2.: x and y are linearly independent if and only if $C (x, y, i) \neq 0$ for all $i = 1, 2, \dots$ .

The last important theorem [8] concerns particular solutions of the non-homogeneous linear difference Equation (4).

Theorem 2.

Consider the second-order non-homogeneous linear difference Equation (4) and its homogeneous counterpart (5) with two linearly independent solutions

x = \{x_{i}, i \in N\}

and

y = \{y_{i}, i \in N\}

. Then, a particular solution of (4) is given as

x_{i, p} = \sum_{k = 1}^{i - 1} γ_{k - 1} \frac{|\begin{matrix} x_{k} & y_{k} \\ x_{i - 1} & y_{i - 1} \end{matrix}|}{|\begin{matrix} x_{k} & y_{k} \\ x_{k + 1} & y_{k + 1} \end{matrix}|}, i = 1, 2, \dots .

(6)

2.2. Discretization of the FitzHugh–Nagumo Model with Diffusion Using the Nearly Exact Scheme

Let

u = u (x, t)

and

v = v (x, t)

be functions of two variables. Consider (1). Using the nearly exact discretization scheme method [4], let

ϕ (x)

and

ψ (t)

be two functions of x and t, respectively. Let

U_{i}^{(n)} = U (x_{i}, t_{n})

, where

x_{i} = i ϕ (x)

and

t_{n} = n ψ (t)

. Let f be a three-times differentiable function. Using Taylor’s polynomial,

\begin{matrix} f (x + h) = f (x) + h f^{'} (x) + \frac{h^{2}}{2} f^{''} (x) + \frac{h^{3}}{6} f^{'''} (ε_{1}) \\ for x < ε_{1} < x + h . \end{matrix}

(7)

Likewise,

\begin{matrix} f (x - h) = f (x) - h f^{'} (x) + \frac{h^{2}}{2} f^{''} (x) - \frac{h^{3}}{6} f^{'''} (ε_{2}) \\ for x - h < ε_{2} < x . \end{matrix}

(8)

Therefore, adding the two equalities above together,

\begin{matrix} f (x + h) - 2 f (x) + f (x - h) = h^{2} f^{''} (x) - \frac{h^{3}}{6} [f^{'''} (ε_{1}) + f^{'''} (ε_{2})] . \end{matrix}

(9)

Thus, using Equation (7) with

h = t_{n + 1} - t_{n} = ϕ (t)

and

f (t) = u (x, t)

\frac{\partial u (x_{i}, t_{n})}{\partial t} \approx \frac{u (x_{i}, t_{n + 1}) - u (x_{i}, t_{n})}{t_{n + 1} - t_{n}} = \frac{U_{i}^{(n + 1)} - U_{i}^{(n)}}{ϕ (t)} .

Likewise, using (9) with

h = ϕ (x)

and

f (x) = u (x, t)

\begin{matrix} \frac{\partial^{2} u (x_{i}, t_{n})}{\partial x^{2}} \approx & \frac{u (x_{i + 1}, t_{n}) - 2 u (x_{i}, t_{n}) + u (x_{i - 1}, t_{n})}{ϕ {(x)}^{2}} \\ = \frac{U_{i + 1}^{(n)} - 2 U_{i}^{(n)} + U_{i - 1}^{(n)}}{ϕ {(x)}^{2}} . \end{matrix}

Given two functions

ϕ_{1} (t)

and

ϕ_{2} (t)

of t, the FitzHugh–Nagumo model with diffusion in Equation (1) becomes

\{\begin{matrix} \frac{U_{i}^{(n + 1)} - U_{i}^{(n)}}{ϕ_{1} (t)} = D \frac{U_{i + 1}^{(n)} - 2 U_{i}^{(n)} + U_{i - 1}^{(n)}}{ϕ {(x)}^{2}} - V_{i}^{(n)} + g (U_{i}^{n}), \\ \frac{V_{i}^{(n + 1)} - V_{i}^{(n)}}{ϕ_{2} (t)} = ε [- β V_{i}^{(n)} + c + U_{i}^{(n)}] . \end{matrix}

This system can be written in a simplified form as

\{\begin{matrix} \frac{U_{i}^{(n + 1)} - U_{i}^{(n)}}{ϕ_{1} (t)} = α W_{i - 1}^{(n)} - V_{i}^{(n)} + g (U_{i}^{(n)}), \\ \frac{V_{i}^{(n + 1)} - V_{i}^{(n)}}{ϕ_{2} (t)} = ε [- β V_{i}^{(n)} + c + U_{i}^{(n)}] . \end{matrix}

where

α = \frac{D}{ϕ {(x)}^{2}}

and

W_{i - 1}^{(n)} = U_{i + 1}^{(n)} - 2 U_{i}^{(n)} + U_{i - 1}^{(n)}

It follows the discrete FitzHugh–Nagumo model with diffusion:

\{\begin{matrix} U_{i}^{(n + 1)} = U_{i}^{(n)} + ϕ_{1} (t) (α W_{i - 1}^{(n)} - V_{i}^{(n)} + g (U_{i}^{(n)})), \\ V_{i}^{(n + 1)} = V_{i}^{(n)} + ϕ_{2} (t) (ε [- β V_{i}^{(n)} + c + U_{i}^{(n)}]) . \end{matrix}

It follows that the fixed (in time n) solutions of this FitzHugh–Nagumo model with diffusion are given by

U_{i} : = U_{i}^{(n + 1)} = U_{i}^{(n)}

and

V_{i} : = V_{i}^{(n + 1)} = V_{i}^{(n)},

for

i = 0, 1 \dots

. This implies

\begin{matrix} V_{i} & = & \frac{1}{β} (c + U_{i}), \\ W_{i - 1} & = & \frac{1}{α} [V_{i} - g (U_{i})] x = \frac{1}{α} [\frac{1}{β} (c + U_{i}) - g (U_{i})] . \end{matrix}

Our next result concerns the fixed solutions of the FitzHugh–Nagumo model with diffusion.

Theorem 3.

Consider the discrete FitzHugh–Nagumo model with diffusion. Suppose

(U_{i}, V_{i})

is given for

i = 0, 1

. Let

(U_{i}, V_{i})

be its fixed solution for

i = 2, 3, \dots

. Then, we have the recursive relations

\begin{matrix} U_{i} & = & U_{0} (i - 1) + i U_{1} + \sum_{k = 1}^{i} \frac{i - k - 1}{α} [\frac{1}{β} (c - U_{k}) - g (U_{k})], \\ V_{i} & = & \frac{1}{β} (c - U_{i}) . i = 2, 3, \dots \end{matrix}

(10)

Proof.

Here, we consider the second-order difference equation in i

U_{i + 1} - 2 U_{i} + U_{i - 1} = W_{i - 1} .

(11)

First, let

U_{i, h}

be a solution of the homogeneous equation

U_{i + 1} - 2 U_{i} + U_{i - 1} = 0 .

(12)

Let

U_{i, p}

be a particular solution of Equation (11) and

U_{i, g}

be its general solution. As a function of i, first solve the homogeneous equation. Since the coefficients are constants, the solution would be of the form

U_{i} = {(q)}^{i}

for a constant q. We obtain the characteristic polynomial

q^{2} - 2 q + 1 = 0 ⟹ q = 1 .

So, the general solutions of the homogeneous system would be of the form

U_{i, h} = A + i B

, where A and B are two constants for

i = 0

U_{0} = A

and for

i = 1

U_{1} = A + B

. Therefore, the solution of the inverse value problem (IVP) is then given as

U_{i, h} = U_{0} (1 - i) + i U_{1}

. The general solution of (11) is then given as

U_{i, g} = U_{0} (1 - i) + i U_{1} + U_{i, p} .

A particular solution

U_{i, p} = U_{i, p} (W_{i})

to Equation (11) remains to be found.

Given two sequences

P = \{P_{i}, i \in N\}

and

Q = \{Q_{i}, i \in N\}

of linearly independent solutions to the homogeneous Equation (12), a particular solution is given as

U_{i, p} = \sum_{k = 1}^{i} W_{k - 1} \frac{|\begin{matrix} P_{k} & Q_{k} \\ P_{i - 1} & Q_{i - 1} \end{matrix}|}{|\begin{matrix} P_{k} & Q_{k} \\ P_{k + 1} & Q_{k + 1} \end{matrix}|}, i = 1, 2, . . .

To choose sequences of linearly independent solutions of Equation (11), we select

P_{i} = A

and

Q_{i} = i B

, where A and B are non-zero constants. By Heymann’s theorem, it suffices to show that their Casoratian is not equal to zero. Indeed,

C (P, Q, i) = |\begin{matrix} P_{i} & Q_{i} \\ P_{i + 1} & Q_{i + 1} \end{matrix}| = P_{i} Q_{i + 1} - Q_{i} P_{i + 1} = (i + 1) A B - i A B = A B \neq 0 .

With these choices of P and Q,

|\begin{matrix} P_{k} & Q_{k} \\ P_{k + 1} & Q_{k + 1} \end{matrix}| = |\begin{matrix} A & k B \\ A & (k + 1) B \end{matrix}| = A B .

Likewise,

|\begin{matrix} P_{k} & Q_{k} \\ P_{i - 1} & Q_{i - 1} \end{matrix}| = |\begin{matrix} A & k B \\ A & (i - 1) B \end{matrix}| = A B (i - k - 1) .

And given that

W_{i - 1} = \frac{1}{α} [\frac{1}{β} (c + U_{i}) - g (U_{i})]

from above, we have

\begin{matrix} U_{i, g} & = & U_{0} (1 - i) + i U_{1} + \sum_{k = 1}^{i} W_{k - 1} (i - k - 1) \\ = & U_{0} (1 - i) + i U_{1} + \sum_{k = 1}^{i} \frac{i - k - 1}{α} [\frac{1}{β} (c + U_{k}) - g (U_{k})] . \end{matrix}

□

2.2.1. Geometry of Steady States

When an integral surface

u (x, t)

(Figure 2) in three-space is hit by the plane

x = i

, the trace is curve

u (i, t)

. Now, although

u = t

is a plane in three-space but when restricted to the plane

x = i

, its trace is a line. Now, the fixed points are precisely the intersection of this line

u = t

and the curve

u (i, t)

It is desirable for the users of neural networks to see the explicit relation between the fast and the slow variables. In the next subsection, we identify the limitation of the geometric method to help them view a phase plane for voltage potential and recovery variable.

2.2.2. Limitation of Geometry

We observe that the discrete FitzHugh–Nagumo model with diffusion

\{\begin{matrix} U_{i}^{(n + 1)} & = U_{i}^{(n)} + ϕ_{1} (t) (α W_{i - 1}^{(n)} - V_{i}^{(n)} + g (U_{i}^{(n)})) \\ V_{i}^{(n + 1)} & = V_{i}^{(n)} + ϕ_{2} (t) (ε [- β V_{i}^{(n)} + c + U_{i}^{(n)}]) \end{matrix},

is a non-linear dynamical system of the form

S_{n + 1} = F (S_{n})

, where for a fixed

i \in N

, we have

S_{n} = (U_{i}^{n}, V_{i}^{n})

and

F (S_{n}) = (f_{1} (S_{n}), f_{2} (S_{n}))

with

\begin{matrix} \{\begin{matrix} f_{1} (S_{n}) & = f_{1} (U_{i}^{n}, V_{i}^{n}) = U_{i}^{(n)} + ϕ_{1} (t) (α W_{i - 1}^{(n)} - V_{i}^{(n)} + g (U_{i}^{(n)})) \\ f_{2} (S_{n}) & = f_{2} (U_{i}^{n}, V_{i}^{n}) = V_{i}^{(n)} + ϕ_{2} (t) (ε [- β V_{i}^{(n)} + c + U_{i}^{(n)}]) \end{matrix} . \end{matrix}

Its Jacobian matrix is therefore given as

J F (u_{i}, v_{i}) = (\begin{matrix} \frac{\partial f_{1}}{\partial u_{i}} & \frac{\partial f_{1}}{\partial v_{i}} \\ \frac{\partial f_{2}}{\partial u_{i}} & \frac{\partial f_{2}}{\partial v_{i}} \end{matrix}),

or after calculation

J F (u_{i}, v_{i}) = (\begin{matrix} 1 + ϕ_{1} (t) (α \frac{\partial W_{i - 1}}{\partial u_{i}} + g^{'} (u)) & - ϕ_{1} (t) \\ ε ϕ_{2} (t) & 1 - ε β ϕ_{2} (t) \end{matrix}) .

(13)

The transformation

S_{n + 1} = F (S_{n})

occurs in the plane

x = i

. F is non-linear. The tangent space under consideration to the manifold around

(u_{i}^{n}, v_{i}^{n})

is the

x = i

plane. The tangent space to the manifold around

(u_{i}^{n + 1}, v_{i}^{n + 1})

is the

x = i

plane again. When we are trying to linearize the FitzHugh–Nagumo model with diffusion, the differential of F, that is, the Jacobian matrix (13), has

u_{i + 1}^{n}

and

u_{i - 1}^{n}

terms. So, we do not remain in the

x = i

plane. Therefore, we need to devise some other method. Recall that for the FitzHugh–Nagumo model [3], it is possible to perform a qualitative analysis, and thus the visualization of the

u v

phase plane is possible. But now, with the diffusion term, it appears that this is not possible.

The following examples of continuous dynamical systems will illustrate this problem. For any generic system of autonomous ODEs,

\{\begin{matrix} \frac{d u}{d t} = f (u, v), \\ \frac{d v}{d t} = g (u, v) . \end{matrix}

with the solution

\begin{matrix} u = \tilde{f} (t), \\ v = \tilde{g} (t) . \end{matrix}

If we eliminate t from this system of solutions, we obtain the phase curves

H (u, v) = 0

. Hence, for a system of coupled ODEs, it is possible to observe the dynamical system in the

u v

phase plane with t covertly as the running parameter for trajectories.

Example 1.

Consider the following system of ODEs:

\{\begin{matrix} \frac{d u}{d t} = v, \\ \frac{d v}{d t} = u . \end{matrix}

This system has the phase space diagram given in Figure 3 below.

We can visualize in Figure 3 the phase curves

u^{2} - v^{2} = k

in the phase plane.

Now, such a treatment is not always possible for a system of PDEs as there are two independent variables. But we can still see the

u v

phase planes for the fixed spatial plane

x = i

For a system of PDEs:

\{\begin{matrix} \frac{\partial u}{\partial t} = f (u, v), \\ \frac{\partial v}{\partial t} = g (u, v) . \end{matrix}

with the solution

u = \tilde{f} (x, t),

v = \tilde{g} (x, t) .

We cannot visualize a dynamical system, but if we fix the plane

x = k

, then we get

u = F (t),

v = G (t) .

Now, by eliminating t, we can see

H (u, v) = 0

in the

u v

plane.

Example 2.

Consider the following system of PDEs:

\{\begin{matrix} \frac{\partial u}{\partial t} = v \\ \frac{\partial v}{\partial t} = u \end{matrix},

for which the general solution could be given as

u (x, t) = F (x) \cdot e^{t} - G (x) \cdot e^{- t}, v (x, t) = F (x) \cdot e^{t} + G (x) \cdot e^{- t}

. We can visualize the

u v

plane if we fix the plane

x = i

for the phase curves

u^{2} - v^{2} = f (i) .

3. Continuous Dynamics

A reason for switching from a discrete FitzHugh–Nagumo model with diffusion to a continuous one in this research is that steady states have a completely different behavior in the case of the FitzHugh–Nagumo model with diffusion when compared to the FitzHugh–Nagumo model without diffusion. While the Hartman–Grobman theorem [8] was considered for the qualitative analysis of this non-linear dynamical system, its applicability is limited to isolated singular points. On the other hand, the current work has established that steady states are not all isolated for (1). Secondly, a celebrated paper on discrete dynamical systems by R. May [9] has indicated that even a simple difference equation can have highly complicated dynamics. Therefore, adding another foreign parameter, the step size, into the system would add to the existing complexity of reaction diffusion systems, which are already known for their chaotic behavior. In addition, real number system was constructed to measure physical quantities on a continuous scale. Numerical analysis never actualizes the potential of real numbers.Above all, if the laws of the universe are written in terms of differential equations, the solutions must firstly be sought in the space of continuous functions.

3.1. Steady States for the Continuous FitzHugh–Nagumo Model with Diffusion

Next, we ask ourselves: what is the geometry of the continuous model

\{\begin{matrix} \frac{\partial u}{\partial t} = D \frac{\partial^{2} u}{\partial x^{2}} - v + g (u) \\ \frac{\partial v}{\partial t} = ε [- β v + c + u] \end{matrix} .

Integral surfaces u and v are expected to be smooth manifolds because they are solutions to the above system of differential equations. There could be an infinite parametric representation for these integral surfaces. The most natural one is obtained by assuming the independent variables

x, t

as parameters. On the constant x-curves, the ordinary time derivatives give us slope of the lineal elements. So, this is the most natural way to view the slope field in the three-space for the solution surfaces u and v. Similarly, in the

v x t

space, we are looking for a surface v so that the lineal elements of this slope field at each point become tangent to the surface v. For (1), u nullclines are

\frac{\partial^{2} u}{\partial x^{2}} - v + g (u) = 0,

(14)

and v nullclines are

ε [- β v + c + u] = 0 .

(15)

Now, by solving (14) and (15) simultaneously, steady states (16) for (1) are obtained.

\begin{matrix} u = RootOf (- \sqrt{6} (\int^{Z} \frac{β D}{\sqrt{β D (f^{4} β + 6 F (t) β D - 6 f^{2} β + 6 f^{2} + 12 f c)}} d f) + x + G (t)), \\ v = \frac{c + RootOf (- \sqrt{6} (\int^{Z} \frac{β D}{\sqrt{β D (f^{4} β + 6 F (t) β D - 6 f^{2} β + 6 f^{2} + 12 f c)}} d f) + x + G (t))}{β} \end{matrix}

(16)

After substituting v from (15) into (14), we get

D \frac{\partial^{2} u}{\partial x^{2}} - \frac{c + u}{β} + g (u) = 0 .

(17)

This is an ordinary differential equation, so the solutions are integral curves. In the

u x t

space, these integral curves are cylindrical surfaces. Similarly, putting solutions of (17) into (15) will result in surfaces in the

v x t

space. The intersection of these surfaces is a curve. The upshot of all these arguments is that (1) also has a continuum of steady states, unlike the FitzHugh–Nagumo model without diffusion. Also note that the steady states of (1) depend on x, as already proven in Section 3. Fixed points for the space-independent FitzHugh–Nagumo model were found in [3]. Some fixed points for (1) are presented in the next section. A good exercise would be to compare discrete steady states (10) with (26).

3.2. Symmetry Analysis

A qualitative analysis of differential equations become mostly meaningless if the exact solution to the system is known. Lie symmetry analysis is one of the most powerful tools available to solve differential equations exactly and systematically, especially for non-linear differential equations, where no general methods are available.

(1) has two obvious Lie symmetries due to the spatio-temporal translation invariance. These two symmetries do not have any differential invariant, so they are not enough to run quadrature. G. Bluman [10] extended the idea of Lie symmetries by introducing the idea of non-classical symmetries.

3.2.1. Non-Classical Symmetry Analysis

Most of the time, it is difficult to tackle a system of differential equations. The two equations in Equation (1) were merged together and as a result, a third-order differential equation was obtained:

D u_{t x x} - u_{t t} + ε β D u_{x x} - ε β u_{t} - ε u + u_{t} (- u^{2} + 1) + ε β (u - \frac{u^{3}}{3}) - ε c = 0 .

(18)

When manual calculations of conditional symmetries [5] were carried out on (18), the only resulting symmetries were the spatio-temporal translations. Usually, in the similarity analysis of reaction diffusion systems, one encounters two second-order partial differential equations. The second equation of (1) is a first-order equation. Considering the differential consequences, the following equivalent system (19) of four equations was obtained. Again, manual calculations of Q-conditional symmetries [11] generated the same two symmetries.

\begin{matrix} u_{t} & = u_{x x} - v + u - u^{3}, \\ v_{t} & = - v + u, \\ u_{x} & = v_{t x} + v_{x}, \\ u_{t} & = v_{t t} + v_{t} . \end{matrix}

(19)

Cherniha and Vasyl [12] have extended this idea of non-classical symmetry. Their method of Q-conditional symmetries also produced only these two symmetries. They have further extended the idea to a so-called no-go area [13]. The result was the same. Recently, Bluman [14] has designed a method based on differential invariants to solve differential equations with known symmetries. As there are no conservation laws for (1), the discovery of this algorithm is significant. Using the translation in space symmetry, (1) is transformed to a locally related system (20). Technically, this locally related system [15] is just an inverted system, so the difficulty of solving the new system is the same as that for (1). Interchanging u and x in (1), we obtain

\begin{matrix} x_{t} & = \frac{x_{u u} - (u - \frac{u^{3}}{3} - v) x_{u}^{3}}{x_{u}^{2}}, \\ \frac{v_{t} x_{u} - x_{t} v_{u}}{2 x_{u}} & = ε (- β v + c + u) . \end{matrix}

(20)

This system (20) possesses an obvious non-useful solution

x = F (t); v = G (t)

. Using the differential invariant of space symmetry, it is possible to find a non-locally related system (21)

\begin{matrix} a_{u} - b_{t} & = 0, \\ b_{u} - \frac{2 a b^{2}}{D} + \frac{4 (v - u + \frac{u^{3}}{3}) b^{3}}{D} & = 0, \\ \frac{v_{t} b - a v_{u}}{2 b} + β ε v - ε u - ε c & = 0 . \end{matrix}

(21)

to (1). Now, there is a systematic relation between the symmetries of this non-locally related system and (1). Hence, there is a procedural relationship between the solutions of this non-locally related system and (1). At this stage, the solutions of (1) appear within reach, as even higher symmetries of one system can have a local correspondent in the other system, but this non-locally related system showed no new symmetries. Although there is a trivial solution

a = F (t), b = 0, v = G (t)

to this non-locally related system, this trivial solution corresponds to the steady states of (1). Steady states were already found using several other techniques. In fact, just two straight quadratures resulted in a singular solution. Symmetries were also exploited to obtain steady states. Moreover, using u and v nullclines, other steady states were determined. The integrability condition of this non-locally related system can further be used to reduce this non-locally related system to a potential system (22). This potential system (22) that is non-locally related to (1) is not easy to solve either.

\begin{matrix} b_{t} & = {(\frac{D b_{u} - 4 (- v + u - \frac{u^{3}}{3}) b^{3}}{2 b^{2}})}_{u}, \\ v_{t} & = ε (- β v + c + u) . \end{matrix}

(22)

Working with the determining equations, it was found that (1) can be solved exactly if a slight change in the structure of the cubic non-linearity is made in the fast variable in the first voltage equation. However, the previous analysis was conducted using a cubic voltage potential in the kinetics of (1); it was decided to look for some other tools rather than altering this cubic non-linearity. There are several modifications of the FitzHugh–Nagumo model with diffusion in the literature. Some of them also consider a diffusion term for the recovery variable. For such FitzHugh–Nagumo models, exact solutions are rather easy to find using symmetries [11]. (1), on the other hand, can be solved exactly if either (20), (21) or (22) can be solved. In fact, a general solution for (1) can be computed if any of these related systems can be solved using the transformations developed in [14,15]. Solving (1) via (18), (19) or (20) is rather difficult, as they are just equivalent to (1).

3.2.2. Similarity Solution

Point symmetry can be used to reduce the number of independent variables in a partial differential equation. The constant c in (1) is mostly neglected [16,17]. This constant c does not affect the symmetry. By ignoring this constant, the following FitzHugh–Nagumo Equation (23) is obtained.

\{\begin{matrix} \frac{\partial u}{\partial t} & = D \frac{\partial^{2} u}{\partial x^{2}} - v + u - \frac{u^{3}}{3} \\ \frac{\partial v}{\partial t} & = ε [- β v + u] \end{matrix},

(23)

Using similarity solutions, a verified group invariant steady state solution (24) was obtained for (23).

\begin{matrix} u (x, t) = c_{2} \sqrt{6} \sqrt{\frac{β - 1}{β c_{2}^{2} + 5 β - 6}} \\ JacobiSN (\frac{(6 c_{1} D β + \sqrt{6} \sqrt{D β (- 6 + 5 β)} x) \sqrt{6} \sqrt{\frac{β - 1}{β c_{2}^{2} + 5 β - 6}}}{6 D β}, \frac{c_{2} \sqrt{5 β^{2} - 6 β}}{- 6 + 5 β}) . \end{matrix}

(24)

3.3. Looking for Traveling Wave Solutions

Spatio-temporal symmetry confirms and Matlab plots [17] indicate the existence of traveling wave solutions to (1). Let us embark on a new journey to search for traveling wave solutions. Although built-in commands in Maple and Mathematica can solve the Nagumo Equation [5] for traveling wave solutions, they cannot evaluate (1). Our handwritten Maple codes for tanh, exp and the Ricatti method were not able to solve (1) completely either. This may be because the second equation of (1) does not include a spatial derivative term. These ansatzes seek soliton-type solutions which propagate in space whilst maintaining their shape over time. This is because the effects of PDE dispersion and dissipative terms tend to cancel out. In the tanh method, a higher exponent of the linear term did not balance out the non-linear term. Wazwaz [18] has suggested a substitution that can convert the rational obtained to an integer as a result of balancing exponents. This substitution does not work for (1). The exponents in the exp method did not balance out either and Maple’s TWSolutions was not able to return even a single fixed point, while our handwritten codes produced steady states. Using these tailored codes, steady-state solutions for (1) were obtained. It is possible that some steady states obtained are the same solutions in disguise.

3.3.1. Exp Method for Traveling Wave Analysis

Although the built-in functions in CAS have an algorithm for the exp method, they do not generate any traveling wave solutions for (1). The exponential method was tried manually, with Mn:=2; Md:=0 set as the order of approximation. Letting

ξ = k (x - b t)

in (18) resulted in the following ode.

\begin{matrix} ode = D (\frac{d^{3}}{d ξ^{3}} U (ξ)) b k^{3} - (\frac{d^{2}}{d ξ^{2}} U (ξ)) b^{2} k^{2} + ϵ β D (\frac{d^{2}}{d ξ^{2}} U (ξ)) k^{2} + ϵ β (\frac{d}{d ξ} U (ξ)) b k \\ - ϵ U (ξ) - (\frac{d}{d ξ} U (ξ)) b k (- U {(ξ)}^{2} + 1) + ϵ β (U (ξ) - \frac{1}{3} U {(ξ)}^{3}) - ϵ c = 0 \end{matrix}

\begin{matrix} F 1 = D (\frac{d^{3}}{d Y^{3}} F (Y)) b k^{3} - (\frac{d^{2}}{d Y^{2}} F (Y)) b^{2} k^{2} + ϵ β D (\frac{d^{2}}{d Y^{2}} F (Y)) k^{2} + ϵ β (\frac{d}{d Y} F (Y)) b k \\ - ϵ F (Y) - (\frac{d}{d Y} F (Y)) b k (- F {(Y)}^{2} + 1) + ϵ β (F (Y) - \frac{1}{3} F {(Y)}^{3}) - ϵ c = 0 \end{matrix}

F (Y) = a_{- 2} e^{- 2 Y} + a_{- 1} e^{- Y} + a_{0} + a_{1} e^{Y} + a_{2} e^{2 Y}

Unknowns = \{b, k, a_{- 2}, a_{- 1}, a_{0}, a_{1}, a_{2}\}

Three solution sets were found for the above unknowns. Two solutions generate complex steady states. Technically, algebraic methods should have produced a steady state, but it is this last solution set that generates a novel singular solution (25).

\begin{matrix} {b = b, k = k, a_{- 2} = 0, a_{- 1} = 0, a_{0} = \frac{1}{2} \frac{{((- 12 c + 4 \sqrt{\frac{9 c^{2} β - 4 β^{3} + 12 β^{2} - 12 β + 4}{β}}) β^{2})}^{1 / 3}}{β} \\ + \frac{2 (β - 1)}{{((- 12 c + 4 \sqrt{\frac{9 c^{2} β - 4 β^{3} + 12 β^{2} - 12 β + 4}{β}}) β^{2})}^{1 / 3}}, a_{1} = 0, a_{2} = 0} \end{matrix}

\{\begin{matrix} u (t, x) & = \frac{{((- 12 c + 4 \sqrt{- \frac{4 β^{3} - 9 c^{2} β - 12 β^{2} + 12 β - 4}{β}}) β^{2})}^{\frac{1}{3}}}{2 β} \\ + \frac{2 (β - 1)}{{((- 12 c + 4 \sqrt{- \frac{4 β^{3} - 9 c^{2} β - 12 β^{2} + 12 β - 4}{β}}) β^{2})}^{\frac{1}{3}}}, \\ v (t, x) & = - \frac{1}{β} [- c - \frac{{((- 12 c + 4 \sqrt{- \frac{4 β^{3} - 9 c^{2} β - 12 β^{2} + 12 β - 4}{β}}) β^{2})}^{\frac{1}{3}}}{2 β}] \\ + \frac{1}{β} [\frac{2 (β - 1)}{{((- 12 c + 4 \sqrt{- \frac{4 β^{3} - 9 c^{2} β - 12 β^{2} + 12 β - 4}{β}}) β^{2})}^{\frac{1}{3}}}] . \end{matrix}

(25)

Stability analysis, chaos, bifurcation, and Turing pattern can now be studied exactly. Once the voltage potential and blocking mechanism obtain this value, the neuron is not going to fire, and there will be no pulse, no wave and no exchange of ions across channels. Geometrically, these are two horizontal planes in the

u x t

and

v x t

space.

3.3.2. Ricatti Method

The Riccatti method of traveling wave solutions for (23), with M:=2 as the order of approximation,

F (Y) = a_{2} Y^{2} + a_{1} Y + a_{0}

and

ξ = k (x - c t)

, when applied to

\begin{matrix} D (\frac{\partial^{3}}{\partial t \partial x^{2}} u) - \frac{\partial^{2}}{\partial t^{2}} u + ε β D (\frac{\partial^{2}}{\partial x^{2}} u) - ε β (\frac{\partial}{\partial t} u) - \\ ε u + (\frac{\partial}{\partial t} u) (- u^{2} + 1) + ε β (u - \frac{1}{3} u^{3}) = 0 . \end{matrix}

produced the following ode.

\begin{matrix} o d e = 2 ϵ β D k^{2} a_{2} {(\frac{d}{d ξ} U (ξ))}^{2} + ϵ β c k a_{1} (\frac{d}{d ξ} U (ξ)) + 5 c k a_{2}^{2} U {(ξ)}^{4} (\frac{d}{d ξ} U (ξ)) a_{1} + \\ 4 c k a_{2}^{2} U {(ξ)}^{3} (\frac{d}{d ξ} U (ξ)) a_{0} + 4 c k a_{2} U {(ξ)}^{3} (\frac{d}{d ξ} U (ξ)) a_{1}^{2} + 2 c k a_{2} U (ξ) (\frac{d}{d ξ} U (ξ)) a_{0}^{2} + \\ 2 c k a_{1}^{2} (\frac{d}{d ξ} U (ξ)) U (ξ) a_{0} - 2 ϵ β a_{2} U {(ξ)}^{3} a_{1} a_{0} - 6 D c k^{3} a_{2} (\frac{d}{d ξ} U (ξ)) (\frac{d^{2}}{d ξ^{2}} U (ξ)) \\ - 2 D c k^{3} a_{2} U (ξ) (\frac{d^{3}}{d ξ^{3}} U (ξ)) + ϵ β D k^{2} a_{1} (\frac{d^{2}}{d ξ^{2}} U (ξ)) - ϵ a_{0} - ϵ a_{2} U {(ξ)}^{2} + \\ ϵ β a_{0} - \frac{ϵ β a_{0}^{3}}{3} - ϵ a_{1} U (ξ) - D c k^{3} a_{1} (\frac{d^{3}}{d ξ^{3}} U (ξ)) - 2 c^{2} k^{2} a_{2} U (ξ) (\frac{d^{2}}{d ξ^{2}} U (ξ)) + \\ 2 c k a_{2}^{3} U {(ξ)}^{5} (\frac{d}{d ξ} U (ξ)) - 2 c k a_{2} U (ξ) (\frac{d}{d ξ} U (ξ)) + c k a_{1}^{3} (\frac{d}{d ξ} U (ξ)) U {(ξ)}^{2} + \\ c k a_{1} (\frac{d}{d ξ} U (ξ)) a_{0}^{2} - ϵ β a_{2}^{2} U {(ξ)}^{5} a_{1} - ϵ β a_{2}^{2} U {(ξ)}^{4} a_{0} - ϵ β a_{2} U {(ξ)}^{4} a_{1}^{2} - \\ ϵ β a_{2} U {(ξ)}^{2} a_{0}^{2} - ϵ β a_{1}^{2} U {(ξ)}^{2} a_{0} - ϵ β a_{1} U (ξ) a_{0}^{2} - c k a_{1} (\frac{d}{d ξ} U (ξ)) + \\ ϵ β a_{2} U {(ξ)}^{2} + ϵ β a_{1} U (ξ) - \frac{ϵ β a_{2}^{3} U {(ξ)}^{6}}{3} - \frac{ϵ β a_{1}^{3} U {(ξ)}^{3}}{3} - c^{2} k^{2} a_{1} (\frac{d^{2}}{d ξ^{2}} U (ξ)) - \\ 2 c^{2} k^{2} a_{2} {(\frac{d}{d ξ} U (ξ))}^{2} + 2 ϵ β c k a_{2} U (ξ) (\frac{d}{d ξ} U (ξ)) + 6 c k a_{2} U {(ξ)}^{2} (\frac{d}{d ξ} U (ξ)) a_{1}_a_{0} + \\ 2 ϵ β D k^{2} a_{2} U (ξ) (\frac{d^{2}}{d ξ^{2}} U (ξ)) = 0 \end{matrix}

U n k n o w n s = {c, k, a_{0}, a_{1}, a_{2}}

For this set of unknowns, twenty-four real steady states were computed. Two of them (26) and (27) are as follows. The steady-state curve (26) in Figure 4 for

ε : = 0.1; β = 0.3; x_{0} = 0; D = 0.7

in terms of the spatial variable is actually an integral surface. This time-independent curve can be viewed as a cylindrical surface in

u x t

space, see Figure 5.

u = - \frac{tanh (\frac{\sqrt{2} \sqrt{\frac{β - 1}{D β}} (x + x_{0})}{2}) \sqrt{3} \sqrt{β - 1}}{\sqrt{β}} .

(26)

u = \frac{tanh (\frac{\sqrt{2} \sqrt{\frac{β - 1}{D β}} (x + x_{0})}{2}) \sqrt{3} \sqrt{β - 1}}{\sqrt{β}} .

(27)

4. Eclectic Dynamics

No theorem can prove the existence or uniqueness of solutions of (1). Hence, it cannot be said that an exact discrete scheme exists for (1). In this section, it has been proven that exact discrete schemes actually exist for (1). Furthermore, the existence of infinite traveling wave solutions for (1) is proven. The problem of viewing the

u v

phase plane for the FitzHugh–Nagumo model with diffusion will also be solved. Since an eclectic approach has been taken to tackle these problems, a new nonstandard term, eclectic dynamics, is coined. Tools from different branches of dynamical systems are used together.

Although the solution surfaces for (1) are in the

u x t

and

v x t

space, the indication that traveling wave solutions exist, which in any case has been proven here, made us transform the coordinate system. Let us put space-time translation symmetry to good use. If a given three-space is cut with a non-horizontal plane, with a direction cosine which is yet to be determined, then the problem can be reduced to ODEs. Now, (1) has been reduced to a system of one second-order and one first-order coupled ordinary differential equations. One algebraic and one differential substitution will change our new system to a third-order ordinary differential equation. When this ODE is transformed to system of first-order ODEs, the system is autonomous (28). Alternatively, let us give a new perspective on traveling wave solutions. Let us assume that the image of the whole line

ξ = x - m t

is a constant u. Moreover, for every possible translation of this line, i.e., for every

ξ

intercept, u assumes constant values continuously. Now, the first question is ’does there exist any slope m to satisfy this condition?’ The answer is yes. Steady-state solutions have already been computed. Hence, we are justified in assuming the existence of such a line in the domain plane. Since (1) does not carry a space or time variable in the reaction part, substitution of this traveling wave variable in (1) transformed it in the system of ODEs (28).

\begin{matrix} \{\begin{matrix} \frac{d}{d ξ} u (ξ) = w (ξ) \\ \frac{d}{d ξ} w (ξ) = m \cdot D \cdot w (ξ) - v (ξ) + u (ξ) - \frac{u {(ξ)}^{3}}{3} \\ \frac{d}{d ξ} v (ξ) = \frac{ϵ}{m} (- β \cdot v (ξ) + u (ξ) + c) \end{matrix} . \end{matrix}

(28)

Notice that w in (28) is just an auxiliary variable. So, by using a command, like scene in Maple, one can isolate the

u v

phase plane from the

u v w

phase space. Ultimately, visualization (Figure 6) of the

u v

phase plane become possible of course, not with time as the running parameter, but

ξ

. Discrete analysis of the FitzHugh–Nagumo model with diffusion did not allow us to view the

u v

phase plane directly.

Theorem 4.

There are infinite traveling wave solutions for (1).

Proof.

Equation (28), which is equivalent to (1), involves only polynomials. So, the existence and uniqueness theorem [19] proves (for

m \neq 0

) the existence of a unique solution for any given initial value. For

m = 0

, steady states have already been computed. □

The right-hand side of the evolution Equation (28) only has polynomial non-linearity. Another important thing to notice is that (28) possesses unique solutions for any non-zero m. So, we can claim that there are infinite traveling wave solutions for (1). The cherry on the top is that the Picard iteration can now be run to obtain solutions to (1) with a priori error analysis. The third-order ode obtained using the traveling wave coordinates and substitution is too sturdy but nevertheless is autonomous. Using

ξ

translation symmetry, (1) has now been reduced to a second-order ode (29). This second-order ordinary differential equation is not solvable exactly. It admits no Lie symmetry. For

c_{1} \neq 0

and

W^{'} = \frac{d}{d T} W

, we have

\begin{matrix} \begin{matrix} 3 D W W^{″} c_{1} c_{2}^{2} - 9 D W^{2} c_{1} c_{2}^{2} + 3 W^{2} (D c_{2}^{2} β ε - c_{1}^{2}) W^{'} \\ + W^{4} (((T^{3} - 3 T) β + 3 T + 3 c) ε W + 3 c_{1} (T^{2} + β ε - 1)) = 0 \end{matrix} . \end{matrix}

(29)

Letting

c_{2} = 0

in the above equation reduces (1) to an unsolvable first-order ordinary differential equation of Abel type. This completes our journey for an analytical solution of (1). Although (1) has now been reduced merely to a first-order ordinary differential equation, the lack of existing analytical methods to solve this Abel-type equation in terms of elementary or special functions reveals the problem in solving (1) exactly.

In the absence of an exact solution for last sixty-two years, the dynamics of (1) have been observed numerically (Section 2). If (28) can be solved exactly, then traveling wave solutions of (1) can be found. This evolutionary system is autonomous and only has

ξ

translation symmetry. Of course, our steady-state solution (25) is also a verified solution to (28). So, this is a tremendous breakthrough. Stability analyses can now be performed without the problem we faced in Section 2. As a future project, the nearly exact discretization scheme could be tested on this system of three first-order coupled autonomous ordinary differential equations (28). Actually, now we are going to prove that an exact discrete scheme exists for (1). The existence and uniqueness of integral curves for (28) prove that these solution curves are orbits of a one-parameter group. Hence, a step-size denominator function can always be found which will never let the output of an iteration leave the solution trajectory no matter what step size is taken. Hence, an exact discrete scheme exists for the FitzHugh–Nagumo model with diffusion.

Asymptotic Analysis

Considering the limited success with exact methods, the presence of an almost negligible parameter

ε

in (1), and a sharp spike in the solutions, singular perturbation theory appears as the most promising approach. Asymptotic analysis of (1) can be carried out in several ways. The simplest technique is to feed the second-order ode (29) into Mathematica^®. Recently, the AsymptoticDsolveValue command has been introduced by Wolfram. While this feature is only applicable to second-order odes, the current work has already reduced (1) to a second-order ode (29). In addition, recent additions to Maple^® also include a similar command, ByPerturbation, for second-order odes. Neither AsymptoticDsolveValue nor ByPerturbation can evaluate (29). Singular perturbation theory can be directly applied to (18). As part of the current work, Maple-supported manual calculations were performed, providing an analytic approximation (30). This heuristic result can be made increasingly more precise by assuming smaller epsilons and possibly higher-order approximations. However, the plotting (see Figure 7 below) and the presence of the tanh function in the closed-form solution are already significant achievements.

\begin{matrix} u (x, t) = \{\begin{matrix} - \frac{tanh (c_{3} t - \frac{\sqrt{3} \sqrt{3 D - \sqrt{- 12 D^{2} c_{3}^{2} + 9 D^{2}}} x}{6 D} + c_{1}) \sqrt{2} \sqrt{D (3 D - \sqrt{- 12 D^{2} c_{3}^{2} + 9 D^{2}})}}{2 D} \\ + \frac{\sqrt{2} \sqrt{D (3 D - \sqrt{- 12 D^{2} c_{3}^{2} + 9 D^{2}})} (\frac{3 D - \sqrt{- 12 D^{2} c_{3}^{2} + 9 D^{2}}}{6 D} - 1)}{2 D c_{3}} + ε \frac{(t + x) (x - t)}{20736 D^{6} c_{3}^{3}} \\ - (\frac{{(x - t)}^{2}}{41472 D^{6} c_{3}^{3}} - \frac{{(t + x)}^{2}}{41472 D^{6} c_{3}^{3}}) ε \end{matrix} \end{matrix}

(30)

5. Discussion

A critique of the nonstandard finite difference scheme was presented. Thence, the FitzHugh–Nagumo model with diffusion was discretized using the nearly exact discretization scheme. Subsequently, fixed points were found. These fixed points were then visualized in space. The underlying problems in carrying out stability analysis were identified. After that, using tools from symmetry analysis, locally related non-locally related and potential systems corresponding to the FitzHugh–Nagumo model with diffusion were discovered. These equivalent systems are alternative routes to solving the FitzHugh–Nagumo model with diffusion analytically. A proven group-invariant steady-state solution was presented. Using symmetry analysis, the FitzHugh–Nagumo model with diffusion was reduced to a second-order and first-order ode. The underlying difficulty in solving the FitzHugh–Nagumo model analytically was identified. Employing different ansatzes, verified steady states were computed. With these results in hand, it is very easy to study the bifurcation, chaos, patterns and stability of the FitzHugh–Nagumo model with diffusion. Then, the existence of infinite traveling wave solutions was proven. It has been proven that there exists an exact discrete scheme corresponding to the FitzHugh–Nagumo model with diffusion. Later, using Maple coding, the relation between the voltage potential u and recovery variable v was determined. In the end, an analytic approximation was presented for the FitzHugh–Nagumo model with diffusion. As a remark, we mention that all the exact solutions in this article were verified using Maple’s pdetest.

Author Contributions

Conceptualization, S. and E.K.; methodology, S. and E.K.; software, S.; validation, S. and E.K.; formal analysis, S. and E.K.; investigation, S. and E.K.; resources, E.K.; writing—original draft preparation, S.; writing—review and editing, E.K.; visualization, S. and E.K.; supervision, E.K. and M.A.; project administration, E.K. and M.A.; funding acquisition, E.K. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by The American Mathematical Society and Simmons Foundation, grant number AMS-SIMMONS-PUI-23028GRT.

Data Availability Statement

The following resources are available on request from the first author: Matlab code for the FitzHugh–Nagumo model with/without diffusion; manual calculations of the related system using symmetry analysis; Maple code for exp and Ricatti methods of traveling wave solutions; the remaining nineteen steady-state solutions of the FitzHugh- Nagumo model; the Maple interface to observe the explicit relation between the voltage potential and recovery variable for the FitzHugh–Nagumo model with diffusion; and Maple-based calculations for the perturbed solution of the FitzHugh–Nagumo model with diffusion.

Acknowledgments

The authors express their gratitude to G. Desouza, who weaves a tapestry of unity among us, through a thread of inspiration. We are also grateful to the reviewers for their suggestions and corrections.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

References

Nagumo, J.; Arimoto, S.; Yoshizawa, S. An active pulse transmission line simulating nerve axon. Proc. IRE 1962, 50, 2061–2070. [Google Scholar] [CrossRef]
Mickens, R.E. Nonstandard Finite Difference Models of Differential Equations; World Scientific: Singapore, 1994. [Google Scholar]
Kwessi, E.; Elaydi, S.; Dennis, B.; Livadiotis, G. Nearly exact discretization of single species population models. Nat. Resour. Model. 2018, 31, e12167. [Google Scholar] [CrossRef]
Kwessi, E.; Edwards, L.J. A Nearly Exact Discretization Scheme for the FitzHugh–Nagumo model. Differ. Equ. Dyn. Syst. 2021, 32, 253–275. [Google Scholar] [CrossRef]
Shih, M.; Momoniat, E.; Mahomed, F.M. Approximate conditional symmetries and approximate solutions of the perturbed FitzHugh-Nagumo equation. J. Math. Phys. 2005, 46, 023503. [Google Scholar] [CrossRef]
Chapwanya, M.; Jejeniwa, O.A.; Appadu, A.R.; Lubuma, J.S. An explicit nonstandard finite difference scheme for the FitzHugh-Nagumo equations. Int. J. Comput. Math. 2018, 96, 1993–2009. [Google Scholar] [CrossRef]
Potts, R.B. Differential and Difference Equations. Am. Math. Mon. 1982, 89, 402–407. [Google Scholar] [CrossRef]
Kulenovic, M.R.S.; Merino, O. Discrete Dynamical Systems and Difference Equations with Mathematica, 1st ed.; Chapman and Hall/CRC: London, UK, 2002. [Google Scholar]
May, R. Simple Mathematical models With Very Complicated Dynamics. Nature 1976, 261, 459–467. [Google Scholar] [CrossRef] [PubMed]
Bluman, G.W. Construction of Solutions to Partial Differential Equations by the Use of Transformation Groups. Ph.D. Thesis, Caltech, Pasadena, CA, USA, 1967. [Google Scholar]
Cherniha, R.; Davydovych, V. Nonlinear Reaction-Diffusion Systems: Conditional Symmetry, Exact Solutions and Their Applications in Biology; Springer International Publishing: Berlin/Heidelberg, Germany, 2017. [Google Scholar]
Cherniha, R.; Davydovych, V. Conditional symmetries and exact solutions of nonlinear reaction-diffusion systems with non-constant diffusivities. Commun. Nonlinear Sci. Numer. Simul. 2012, 17, 3177–3188. [Google Scholar] [CrossRef]
Cherniha, R.; Davydovych, V. A hunter-gatherer–farmer population model: Lie symmetries, exact solutions and their interpretation. Eur. J. Appl. Math. 2018, 30, 338–357. [Google Scholar] [CrossRef]
Bluman, G.W.; de la Rosa, R.; Santos Bruzón, M.; Luz Gandarias, M. Differential invariant method for seeking nonlocally related systems and nonlocal symmetries. I: General theory and examples. Proc. R. Soc. 2021, 477, 20200909. [Google Scholar] [CrossRef]
Bluman, G.W.; Yang, Z. A symmetry-based method for constructing nonlocally related partial differential equation systems. J. Math. Phys. 2013, 54, 093504. [Google Scholar] [CrossRef]
Battista, N.A.; Miller, L.A. Bifurcations in valveless pumping techniques from a coupled fluid–structure-electrophysiology model in heart development. arXiv 2017, arXiv:1709.05264. [Google Scholar] [CrossRef]
Schiesser, W.E.; Griffiths, G. Traveling Wave Analysis of Partial Differential Equations: Numerical and Analytical Methods with Matlab and Maple; Elsevier Science: Amsterdam, The Netherlands, 2010. [Google Scholar]
Wazwaz, A.-M. The tanh method for traveling wave solutions of nonlinear equations. Appl. Math. Comput. 2004, 154, 713–723. [Google Scholar] [CrossRef]
Strogatz, S.H. Nonlinear Dynamics and Chaos: With Application to Physics, Biology, Chemistry and Engineering; Perseus Books Publishing: Cambridge, UK, 2024. [Google Scholar]

Figure 1. Solution surface for some partial differential equation and a solution set of the corresponding difference equation.

Figure 2. For the integral surface, two fixed points in the

x = i

plane are shown. The planes

u = t

and

x = i

meet on the line shown. The trace of

u (x, t)

in the

x = i

plane is shown in the dotted curve.

Figure 2. For the integral surface, two fixed points in the

x = i

plane are shown. The planes

u = t

and

x = i

meet on the line shown. The trace of

u (x, t)

in the

x = i

plane is shown in the dotted curve.

Figure 3. Evolution of trajectories in the phase plane. Only upper branches of these rectangular hyperbolas are shown here, with t as a covert variable.

Figure 4. A steady-state solution of the FitzHugh-Nagumo model with only the spatial profile shown.

Figure 5. Cylindrical surface view of the steady-state solution of the FitzHugh-Nagumo model.

Figure 6.

u v

phase curve from Equation (28); [u(0) = 0, w(0) = 2, v(0) = 0]; c = −0.001293;

ε

= 0.008; m = 10;

β

= 2.54; D = 0.03.

Figure 6.

u v

phase curve from Equation (28); [u(0) = 0, w(0) = 2, v(0) = 0]; c = −0.001293;

ε

= 0.008; m = 10;

β

= 2.54; D = 0.03.

Figure 7. Analytic approximation Equation (30) with

D = 0.1; c_{3} = 0.03; c_{1} = 0; ε = 0.001

Figure 7. Analytic approximation Equation (30) with

D = 0.1; c_{3} = 0.03; c_{1} = 0; ε = 0.001

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 MDPI and/or the editor(s). MDPI 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.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Shahid; Abbas, M.; Kwessi, E. Nonstandard Nearly Exact Analysis of the FitzHugh–Nagumo Model. Symmetry 2024, 16, 585. https://doi.org/10.3390/sym16050585

AMA Style

Shahid, Abbas M, Kwessi E. Nonstandard Nearly Exact Analysis of the FitzHugh–Nagumo Model. Symmetry. 2024; 16(5):585. https://doi.org/10.3390/sym16050585

Chicago/Turabian Style

Shahid, Mujahid Abbas, and Eddy Kwessi. 2024. "Nonstandard Nearly Exact Analysis of the FitzHugh–Nagumo Model" Symmetry 16, no. 5: 585. https://doi.org/10.3390/sym16050585

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Nonstandard Nearly Exact Analysis of the FitzHugh–Nagumo Model

Abstract

1. Introduction

2. Discrete Dynamics

2.1. Review of Second-Order Linear Difference Equations

2.2. Discretization of the FitzHugh–Nagumo Model with Diffusion Using the Nearly Exact Scheme

2.2.1. Geometry of Steady States

2.2.2. Limitation of Geometry

3. Continuous Dynamics

3.1. Steady States for the Continuous FitzHugh–Nagumo Model with Diffusion

3.2. Symmetry Analysis

3.2.1. Non-Classical Symmetry Analysis

3.2.2. Similarity Solution

3.3. Looking for Traveling Wave Solutions

3.3.1. Exp Method for Traveling Wave Analysis

3.3.2. Ricatti Method

4. Eclectic Dynamics

Asymptotic Analysis

5. Discussion

Author Contributions

Funding

Data Availability Statement

Acknowledgments

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI