Novel Hopf Bifurcation Exploration and Control Strategies in the Fractional-Order FitzHugh–Nagumo Neural Model Incorporating Delay

Abstract

In this article, we propose a new fractional-order delay-coupled FitzHugh–Nagumo neural model. Taking advantage of delay as a bifurcation parameter, we explore the stability and bifurcation of the formulated fractional-order delay-coupled FitzHugh–Nagumo neural model. A delay-independent stability and bifurcation conditions for the fractional-order delay-coupled FitzHugh–Nagumo neural model is acquired. By designing a proper $P{D}^p$ controller, we can efficaciously control the stability domain and the time of emergence of the bifurcation phenomenon of the considered fractional delay-coupled FitzHugh–Nagumo neural model. By exploiting a reasonable hybrid controller, we can successfully adjust the stability domain and the bifurcation onset time of the involved fractional delay-coupled FitzHugh–Nagumo neural model. This study shows that when the delay crosses a critical value, a Hopf bifurcation will arise. When we adjust the control parameter, we can find other critical values to enlarge or narrow the stability domain of the fractional-order delay-coupled FitzHugh–Nagumo neural model. In order to check the correctness of the acquired outcomes of this article, we present some simulation outcomes via Matlab 7.0 software. The obtained theoretical fruits in this article have momentous theoretical significance in running and constructing networks.

1. Introduction

The human brain consists of thousands of neurons that are connected with each other in a very complex network of neurons [1]. In order to grasp the inherent essence of neurons, many researchers have shown significant interest in establishing appropriate neural network models to uncover the laws of interaction among them. Generally speaking, the following canonical models are often used to describe the interaction among different neurons: the Hodgkin–Huxley model [2], the Morris–Lecar model [3], the Hindmarsh–Rose model [4], the FitzHugh–Nagumo model [5], and the integrate-and-fire model [6,7]. In particular, the FitzHugh–Nagumo model is an important model that displays strong oscillation and bifurcation under appropriate parameters. The FitzHugh–Nagumo model comes from the simplified version of the canonical Hodgkin–Huxley model [1] that depicts the dynamics of neurons and the dynamics of excitable systems in various areas, such as solid-state physics and chemical reaction kinetics [8]. The FitzHugh–Nagumo model has been greatly used in computational neuroscience and nonlinear dynamics [1]. Nowadays, many valuable works on the dynamics of various FitzHugh–Nagumo models have been carried out and many excellent fruits have been reported. For example, Demina and Kudryashov [9] dealt with the meromorphic solution of a FitzHugh–Nagumo system via Nevanlinna’s theory of differential systems. He et al. [10] explored the existence and stability of order-1 and order-2 periodic solutions to a FitzHugh–Nagumo neuron system, with state-dependent impulses via the fixed point of Poincaré mapping and geometric analysis skills. Zheng and Shen [11] carried out an analysis of Turing instability caused by the random network in the FitzHugh–Nagumo system. Gani and Ogawa [12] focused on the instability of periodic traveling wave solutions to a FitzHugh–Nagumo system in excitable media. For more related topics on this aspect, one can see [13,14,15,16,17].

In 2015, Anderson Hoff et al. [18] dealt with the numerical bifurcation issue of the following coupled FitzHugh–Nagumo neural model:

$$\left\{\begin{array}{c}{\displaystyle \frac{d{w}_1\left(t\right)}{dt}=a\left(w_2+w_1-\frac{1}{3}{w}_1^3\right)+b(w_1-w_3),}\hfill \\ {\displaystyle \frac{d{w}_2\left(t\right)}{dt}=-\frac{1}{a}(w_1-\alpha +\beta w_2),}\hfill \\ {\displaystyle \frac{d{w}_3\left(t\right)}{dt}=a\left(w_4+w_3-\frac{1}{3}{w}_3^3\right)+b(w_3-w_1),}\hfill \\ {\displaystyle \frac{d{w}_4\left(t\right)}{dt}=-\frac{1}{a}(w_3-\alpha +\beta w_4),}\hfill \end{array}\right.$$

(1)

where $w_1\left(t\right)\mathrm{and}{w}_3\left(t\right)$ denote the voltage across the cell membrane of a neuron at time t; $w_2\left(t\right)\mathrm{and}{w}_4\left(t\right)$ denote the recovery state of the resting membrane of a neuron at time t; $\alpha ,\beta ,a$ denote parameters; and b denotes the coupling strength between the network elements. For more details, one can see [18,19].

In many cases, time delay often exists in FitzHugh–Nagumo neural systems due to the postponement of signal propagations of different neurons. Based on this viewpoint, Jia et al. [20] proposed the following delay-coupled FitzHugh–Nagumo neural model:

$$\left\{\begin{array}{c}{\displaystyle \frac{d{w}_1\left(t\right)}{dt}=a\left(w_2+w_1-\frac{1}{3}{w}_1^3\right)+b(w_1-w_3(t-\theta )),}\hfill \\ {\displaystyle \frac{d{w}_2\left(t\right)}{dt}=-\frac{1}{a}(w_1-\alpha +\beta w_2),}\hfill \\ {\displaystyle \frac{d{w}_3\left(t\right)}{dt}=a\left(w_4+w_3-\frac{1}{3}{w}_3^3\right)+b(w_3-w_1(t-\theta )),}\hfill \\ {\displaystyle \frac{d{w}_4\left(t\right)}{dt}=-\frac{1}{a}(w_3-\alpha +\beta w_4),}\hfill \end{array}\right.$$

(2)

where $\theta$ stands for the delay. By discussing the characteristic equation of system (2), Jia et al. [20] set up the stability and bifurcation conditions for system (2). Moreover, due to the normal form theorem and center manifold theory, the Hopf bifurcation nature is analyzed.

Note that all the works explored above merely deal with the first-order FitzHugh–Nagumo neural models. At present, fractional calculus is a very valuable tool used in depicting the memory traits and hereditary peculiarities of different materials and evolution processes [21,22,23,24,25,26,27,28]. Fractional calculus has been used in many areas, such as mechanics, different kinds of physical waves, finance, neural networks, secure cryptography, auto-control, etc. [29,30]. Recently, abundant studies on fractional dynamical models have been published. For instance, Li and Yan [31] discussed the Hopf bifurcation control issue in a fractional delayed predator–prey system with disease and cannibalism. Kao and Li [32] explored the asymptotic multistability and local S-asymptotic $\omega$-periodicity in fractional neural network models with impulse. Luo et al. [33] handled the fixed-time control issue in a fractional chaotic model due to the backstepping approach. Jin et al. [34] set up delay-dependent and order-dependent criteria on stability and stabilization in fractional memristive neural network models with delay. For more details, one can see [35,36,37,38,39].

Relying on the exploration above and in order to depict the memory trait and hereditary peculiarities of the voltage across the cell membrane of a neuron and the recovery state of the resting membrane of a neuron, we modify model (2) as the fractional-order delay-coupled FitzHugh–Nagumo neural model, as follows:

$$\left\{\begin{array}{c}{\displaystyle \frac{d^q{w}_1\left(t\right)}{d{t}^q}=a\left(w_2+w_1-\frac{1}{3}{w}_1^3\right)+b(w_1-w_3(t-\theta )),}\hfill \\ {\displaystyle \frac{d^q{w}_2\left(t\right)}{d{t}^q}=-\frac{1}{a}(w_1-\alpha +\beta w_2),}\hfill \\ {\displaystyle \frac{d^q{w}_3\left(t\right)}{d{t}^q}=a\left(w_4+w_3-\frac{1}{3}{w}_3^3\right)+b(w_3-w_1(t-\theta )),}\hfill \\ {\displaystyle \frac{d^q{w}_4\left(t\right)}{d{t}^q}=-\frac{1}{a}(w_3-\alpha +\beta w_4),}\hfill \end{array}\right.$$

(3)

$0<q\le 1$ denotes a constant. All other parameters of model (3) have identical implications to those in model (1). Clearly, model (3) and model (1) own the identical equilibrium point.

Periodic oscillation is a vital dynamical peculiarity in delayed FitzHugh–Nagumo neural models. The delay-driven Hopf bifurcation is a special periodic oscillation. The delay-driven Hopf bifurcation plays a vital role in designing neural networks. Thus, it has received great attention from many scholars. However, many works on delay-driven Hopf bifurcation are only concerned with integer-order dynamical models. The study on the delay-driven Hopf bifurcation of fractional-order dynamical models is reversely rare. Recently, some articles have focused on the delay-driven Hopf bifurcation of fractional-order dynamical systems. For example, Amine et al. [40] analyzed the stability and Hopf bifurcation of a fractional tumor virotherapy system with delay. Xu et al. [41] revealed the impact of both delays on Hopf bifurcation in fractional four-dimensional neural networks. Huang et al. [42] considered the stability and bifurcation behaviors of fractional-order ring-structured neural networks. Alidousti [43] investigated the stability and delay-driven Hopf bifurcation in a fractional-order predator–prey scavenger system. For more details, one can see [44,45,46,47,48,49,50].

While some authors have investigated the delay-driven Hopf bifurcation in fractional-order dynamical models, numerous unresolved questions still remain to be addressed. For example, the challenge of designing a suitable controller to adjust the stability domain and the timing of bifurcation generation in fractional-order delayed FitzHugh–Nagumo neural models stands out. This motivates us to seek the delay-driven Hopf bifurcation and its control issue of fractional-order FitzHugh–Nagumo neural models.

In this work, we explore the following key aspects:

(i)

The stability and bifurcation phenomena of system (3).

(ii)

The $P{D}^p$ controller used to control the stability domain and the timing of the bifurcation generation of system (3).

(iii)

A hybrid controller to adjust the stability domain and the timing of the bifurcation generation of system (3).

The elementary framework in this article is as follows. Section 2 provides a basic theory about fractional dynamical equations. Section 3 explores the stability and bifurcation problem of model (3). Section 4 explores the Hopf bifurcation control issue of model (3) due to the $P{D}^p$ controller. Section 5 investigates the Hopf bifurcation control issue of model (3) due to the hybrid controller. Section 6 carries out software simulations to verify the correctness of the acquired key outcomes of this article. Section 7 concludes this article.

2. Preliminaries

In this segment, we list some necessary definitions and lemmas on fractional calculus, which will be applied in the following proof. Denote $R_{+}=\left\{y\right|y\ge 0,y\in R\}.$

Definition 1 

([51]). Define the Caputo fractional-order derivative as follows:

$${\mathcal{D}}^{q}l\left(\chi \right)=\frac{1}{\Gamma (k-q)}{\int }_{{\chi }_0}^{\chi }\frac{l^{\left(k\right)}\left(s\right)}{{(\chi -s)}^{q-k+1}}ds,$$

where $l\left(\chi \right)\in ([{\chi }_0,\infty ),R),\Gamma \left(s\right)={\int }_0^{\infty }{\chi }^{s-1}{e}^{-\chi }d\chi ,$ $\chi \ge {\chi }_0$ and $k\in Z^{+},k-1\le q<k.$

The Laplace transform of the Caputo fractional-order derivative is given by

$$\mathcal{L}\{{\mathcal{D}}^{q}l\left(t\right);s\}=s^q\mathcal{L}\left(s\right)-\sum _{j=0}^{m-1}{s}^{q-j-1}{l}^{\left(j\right)}\left(0\right),m-1\le q<m\in Z^{+},$$

where $\mathcal{L}\left(s\right)=\mathcal{L}\left\{l\right(t\left)\right\}.$ Especially, if $l^{\left(j\right)}\left(0\right)=0,j=1,2,\cdots ,m,$ then $\mathcal{L}\{{\mathcal{D}}^{q}l\left(t\right);s\}=s^q\mathcal{L}\left(s\right).$

Definition 2 

([52]). $(w_{1\ast },w_{2\ast },w_{3\ast },w_{4\ast })$ is said to be an equilibrium point of system (3) if

$$\left\{\begin{array}{c}{\displaystyle a\left(w_{2\ast }+w_{1\ast }-\frac{1}{3}{w}_{1\ast }^3\right)+b(w_{1\ast }-w_{3\ast })=0,}\hfill \\ {\displaystyle -\frac{1}{a}(w_{1\ast }-\alpha +\beta w_{2\ast })=0,}\hfill \\ {\displaystyle a\left(w_{4\ast }+w_{3\ast }-\frac{1}{3}{w}_{3\ast }^3\right)+b(w_{3\ast }-w_{1\ast })=0,}\hfill \\ {\displaystyle -\frac{1}{a}(w_{3\ast }-\alpha +\beta w_{4\ast })=0}\hfill \end{array}\right.$$

(4)

Lemma 1 

([53,54]). Consider the fractional-order dynamical system as follows:

$${\displaystyle \frac{d^{q}u\left(t\right)}{d{t}^q}=l(t,u\left(t\right)),u\left(0\right)=u_0,}$$

(5)

where $q\in (0,1]$ and $l(t,u\left(t\right)):R^{+}\times R^n\to R^n$. The equilibrium point of system (5) is locally asymptotically stable if every eigenvalue ω of $\frac{\partial l(t,u)}{\partial u}$ evaluated near the equilibrium point obeys $\left|arg\left(\omega \right)\right|>\frac{q\pi }{2}$.

Lemma 2 

([55]). Consider the following fractional-order dynamical system as follows:

$$\left\{\begin{array}{c}{\displaystyle \frac{d^{q_1}{\mathcal{V}}_1\left(t\right)}{d{t}^{q_1}}=l_{11}{\mathcal{V}}_1(t-{\theta }_{11})+l_{12}{\mathcal{T}}_2(t-{\theta }_{12})+\cdots +l_{1n}{\mathcal{V}}_n(t-{\theta }_{1n}),}\hfill \\ {\displaystyle \frac{d^{q_2}{\mathcal{V}}_2\left(t\right)}{d{t}^{q_2}}=l_{21}{\mathcal{V}}_1(t-{\theta }_{21})+l_{22}{\mathcal{T}}_2(t-{\theta }_{22})+\cdots +l_{2n}{\mathcal{V}}_n(t-{\theta }_{2n}),}\hfill \\ {\displaystyle ⋮}\hfill \\ {\displaystyle \frac{d^{q_n}{\mathcal{V}}_n\left(t\right)}{d{t}^{q_n}}=l_{n1}{\mathcal{V}}_1(t-{\theta }_{n1})+l_{n2}{\mathcal{T}}_2(t-{\theta }_{n2})+\cdots +l_{nn}{\mathcal{V}}_n(t-{\theta }_{nn}),}\hfill \end{array}\right.$$

(6)

where $0<q_i<1(i=1,2,\cdots ,n)$, the initial values ${\mathcal{V}}_k\left(t\right)={\omega }_k\left(t\right)\in C[-{max}_{k,l}{\theta }_{kl},0],$ $t\in [-{max}_{k,l}{\theta }_{kl},0],$ $k,l=1,2,\cdots ,n.$ Set

$$\Delta \left(\nu \right)=\left[\begin{array}{cccc}{\nu }^{q_1}-l_{11}{e}^{-\nu {\theta }_{11}}& -l_{12}{e}^{-\nu {\theta }_{12}}& \cdots & -l_{1n}{e}^{-\nu {\theta }_{1n}}\\ -l_{21}{e}^{-\nu {\theta }_{12}}& {\nu }^{q_2}-l_{22}{e}^{-\nu {\theta }_{22}}& \cdots & -l_{2n}{e}^{-\nu {\theta }_{2n}}\\ ⋮& ⋮& \ddots & ⋮\\ -l_{n1}{e}^{-\nu {\theta }_{n1}}& -l_{n2}{e}^{-\nu {\theta }_{n2}}& \cdots & {\zeta }^{q_n}-l_{nn}{e}^{-\nu {\theta }_{nn}}\end{array}\right],$$

(7)

then the null solution of Equation (7) is said to be Lyapunov-asymptotically stable if every root of $det(\Delta (\nu \left)\right)=0$ admits a negative real part.

3. Bifurcation Issue

In this segment, we explore the stability trait and Hopf bifurcation phenomenon of system (3). According to [20], we know that if

$$\left(A_1\right)\alpha >\sqrt{3},\beta \in (0,1)$$

holds, then system (3) admits a unique positive equilibrium point $W_0(w_{1\ast },w_{2\ast },w_{3\ast },w_{4\ast })$, where $w_{1\ast },w_{2\ast },$ $w_{3\ast },w_{4\ast }$ obey

$$\left\{\begin{array}{c}{\displaystyle a\left(w_{2\ast }+w_{1\ast }-\frac{1}{3}{w}_{1\ast }^3\right)+b(w_{1\ast }-w_{3\ast })=0,}\hfill \\ {\displaystyle -\frac{1}{a}(w_{1\ast }-\alpha +\beta w_{2\ast })=0,}\hfill \\ {\displaystyle a\left(w_{4\ast }+w_{3\ast }-\frac{1}{3}{w}_{3\ast }^3\right)+b(w_{3\ast }-w_{1\ast })=0,}\hfill \\ {\displaystyle -\frac{1}{a}(w_{3\ast }-\alpha +\beta w_{4\ast })=0.}\hfill \end{array}\right.$$

(8)

The linear system of Equation (8) near $W_0(w_{1\ast },w_{2\ast },w_{3\ast },w_{4\ast })$ can be expressed as

$$\left\{\begin{array}{c}{\displaystyle \frac{d^q{w}_1\left(t\right)}{d{t}^q}=(a+b-a{w}_{1\ast }^2)w_1\left(t\right)+a{w}_2\left(t\right)-b{w}_3(t-\theta ),}\hfill \\ {\displaystyle \frac{d^q{w}_2\left(t\right)}{d{t}^q}=-\frac{1}{a}{w}_1\left(t\right)-\frac{\beta }{a}{w}_2\left(t\right),}\hfill \\ {\displaystyle \frac{d^q{w}_3\left(t\right)}{d{t}^q}=(a+b-a{w}_{3\ast }^2)w_3\left(t\right)+a{w}_4\left(t\right)-b{w}_1(t-\theta ),}\hfill \\ {\displaystyle \frac{d^q{w}_4\left(t\right)}{d{t}^q}=-\frac{1}{a}{w}_3\left(t\right)-\frac{\beta }{a}{w}_4\left(t\right).}\hfill \end{array}\right.$$

(9)

The associated characteristic equation of (9) reads as follows:

$$det\left[\begin{array}{cccc}{s}^q-(a+b-a{w}_{1\ast }^2)& -a& b{e}^{-s\theta }& 0\\ \frac{1}{a}& s^q+\frac{\beta }{a}& 0& 0\\ b{e}^{-s\theta }& 0& s^q-(a+b-a{w}_{3\ast }^2)& -a\\ 0& 0& \frac{1}{a}& s^q+\frac{\beta }{a}\end{array}\right]=0.$$

(10)

Applying (10), one obtains

$$s^{4q}+a_1{s}^{3q}+a_2{s}^{2q}+a_3{s}^q+a_4+(a_5{s}^{2q}+a_6{s}^q+a_7)e^{-2s\theta }=0,$$

(11)

where

$$\left\{\begin{array}{c}{\displaystyle a_1=\frac{2\beta }{a}-2a-2b+a(w_{1\ast }^2+w_{3\ast }^2),}\hfill \\ {\displaystyle a_2=\frac{{\beta }^2}{a^2}-\frac{2\beta }{a}(a+b-a{w}_{1\ast }^2)-\frac{2\beta }{a}(a+b-a{w}_{3\ast }^2)}\hfill \\ {\displaystyle +(a+b-a{w}_{1\ast }^2)(a+b-a{w}_{3\ast }^2)+2,}\hfill \\ {\displaystyle a_3=\frac{2\beta }{a}-(a+b-a{w}_{1\ast }^2)-\frac{{\beta }^2}{a^2}a+b-a{w}_{1\ast }^2)}\hfill \\ {\displaystyle +\frac{\beta }{a}(a+b-a{w}_{1\ast }^2)(a+b-a{w}_{3\ast }^2)-\frac{{\beta }^2}{a^2}(a+b-a{w}_{3\ast }^2)}\hfill \\ {\displaystyle +\frac{\beta }{a}(a+b-a{w}_{1\ast }^2)(a+b-a{w}_{3\ast }^2)-(a+b-a{w}_{3\ast }^2),}\hfill \\ {\displaystyle a_4=1-\frac{\beta }{a}(a+b-a{w}_{1\ast }^2)+\frac{{\beta }^2}{a^2}(a+b-a{w}_{1\ast }^2)(a+b-a{w}_{3\ast }^2)}\hfill \\ {\displaystyle -\frac{\beta }{a}(a+b-a{w}_{3\ast }^2),}\hfill \\ {\displaystyle a_5=-b^2,}\hfill \\ {\displaystyle a_6=\frac{2\beta b^2}{a},}\hfill \\ {\displaystyle a_7=\frac{{\beta }^2}{a^2}.}\hfill \end{array}\right.$$

(12)

Assuming that $s=i\varphi =\varphi \left(\cos \frac{\pi }{2}+isin\frac{\pi }{2}\right)$ is a root of (11), we obtain

$$\begin{array}{ccc}& & {\varphi }^{4q}(\cos 2q\pi +i\sin 2q\pi )+a_1{\varphi }^{3q}\left(\cos \frac{3q\pi }{2}+i\sin \frac{3q\pi }{2}\right)\hfill \\ & & +a_2{\varphi }^{2q}(\cos q\pi +i\sin q\pi )+a_3{\varphi }^q\left(\cos \frac{q\pi }{2}+i\sin \frac{q\pi }{2}\right)\hfill \\ & & +a_4+\left[a_5{\varphi }^{2q}(\cos q\pi +i\sin q\pi )+a_6{\varphi }^q\left(\cos \frac{q\pi }{2}+i\sin \frac{q\pi }{2}\right)+a_7\right]\hfill \\ & & \times (\cos 2\varphi \theta -i\sin 2\varphi \theta )=0,\hfill \end{array}$$

(13)

which leads to

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{A}}_1\left(\varphi \right)\cos 2\varphi \theta +{\mathcal{A}}_2\left(\varphi \right)\cos 2\varphi \theta =-{\mathcal{A}}_3\left(\varphi \right),}\hfill \\ {\displaystyle {\mathcal{A}}_2\left(\varphi \right)\cos 2\varphi \theta -{\mathcal{A}}_1\left(\varphi \right)\sin 2\varphi \theta =-{\mathcal{A}}_4\left(\varphi \right),}\hfill \end{array}\right.$$

(14)

where

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{A}}_1\left(\varphi \right)=b_1{\varphi }^{2q}+b_2{\varphi }^q+b_3,}\hfill \\ {\displaystyle {\mathcal{A}}_1\left(\varphi \right)=b_4{\varphi }^{2q}+b_5{\varphi }^q,}\hfill \\ {\displaystyle {\mathcal{A}}_1\left(\varphi \right)=b_6{\varphi }^{4q}+b_7{\varphi }^{3q}+b_8{\varphi }^{2q}+b_9{\varphi }^q+b_{10},}\hfill \\ {\displaystyle {\mathcal{A}}_1\left(\varphi \right)=b_{11}{\varphi }^{4q}+b_{12}{\varphi }^{3q}+b_{13}{\varphi }^{2q}+b_{14}{\varphi }^q,}\hfill \end{array}\right.$$

(15)

where

$$\left\{\begin{array}{c}{\displaystyle b_1=a_5\cos q\pi ,}\hfill \\ {\displaystyle b_2=a_6\cos \frac{2q\pi }{2},}\hfill \\ {\displaystyle b_3=a_7,}\hfill \\ {\displaystyle b_4=a_5\sin q\pi ,}\hfill \\ {\displaystyle b_5=a_6\sin \frac{2q\pi }{2},}\hfill \\ {\displaystyle b_6=\cos 2q\pi ,}\hfill \\ {\displaystyle b_7=a_1\cos \frac{3q\pi }{2},}\hfill \end{array}\right.$$

(16)

and

$$\left\{\begin{array}{c}{\displaystyle b_8=a_2\cos q\pi ,}\hfill \\ {\displaystyle b_9=a_3\cos \frac{2q\pi }{2},}\hfill \\ {\displaystyle b_{10}=a_4,}\hfill \\ {\displaystyle b_{11}=\sin 2q\pi ,}\hfill \\ {\displaystyle b_{12}=a_1\sin \frac{3q\pi }{2},}\hfill \\ {\displaystyle b_{13}=a_2\sin q\pi ,}\hfill \\ {\displaystyle b_{14}=a_3\sin \frac{2q\pi }{2}.}\hfill \end{array}\right.$$

(17)

It follows from (14) that

$${\mathcal{A}}_1^2\left(\varphi \right)+{\mathcal{A}}_2^2\left(\varphi \right)={\mathcal{A}}_3^2\left(\varphi \right)+{\mathcal{A}}_4^2\left(\varphi \right),$$

(18)

which results in

$$c_1{\varphi }^{8q}+c_2{\varphi }^{7q}+c_3{\varphi }^{6q}+c_4{\varphi }^{5q}+c_5{\varphi }^{4q}+c_6{\varphi }^{3q}+c_7{\varphi }^{2q}+c_8{\varphi }^q+c_9=0,$$

(19)

where

$$\left\{\begin{array}{c}{\displaystyle c_1=b_6^2+b_{11}^2,}\hfill \\ {\displaystyle c_2=2(b_6{b}_7+b_{11}{b}_{12}),}\hfill \\ {\displaystyle c_3=b_7^2+b_{12}^2+2(b_6{b}_8+b_{11}{b}_{13}),}\hfill \\ {\displaystyle c_4=2(b_6{b}_9+b_7{b}_8+b_{11}{b}_{14}+b_{12}{b}_{13}),}\hfill \\ {\displaystyle c_5=b_8^2+b_{13}^2-b_1^2-b_4^2+2(b_6{b}_{10}+b_7{b}_9+b_{12}{b}_{13}),}\hfill \\ {\displaystyle c_6=2(b_7{b}_{10}+b_8{b}_9+b_{13}{b}_{14}-2b_1{b}_2-2b_4{b}_5),}\hfill \\ {\displaystyle c_7=b_9^2+b_{14}^2-b_5^2-b_2^2+2(b_8{b}_{10}-b_1{b}_3),}\hfill \\ {\displaystyle c_8=2(b_9{b}_{10}-b_2{b}_3),}\hfill \\ {\displaystyle c_9=b_{10}^2-b_3^2.}\hfill \end{array}\right.$$

(20)

Denote

$$Q\left(\varphi \right)=c_1{\varphi }^{8q}+c_2{\varphi }^{7q}+c_3{\varphi }^{6q}+c_4{\varphi }^{5q}+c_5{\varphi }^{4q}+c_6{\varphi }^{3q}+c_7{\varphi }^{2q}+c_8{\varphi }^q+c_9$$

(21)

and

$$S\left(\varphi \right)=c_1{\varphi }^8+c_2{\varphi }^7+c_3{\varphi }^6+c_4{\varphi }^5+c_5{\varphi }^4+c_6{\varphi }^3+c_7{\varphi }^2+c_8\varphi +c_9$$

(22)

Lemma 3. 

(i)

Suppose that $a_4+a_7\ne 0$ and $c_k>0(k=1,2,\cdots ,9)$, then Equation (11) does not admit any roots with a zero real part.

(ii)

Suppose that $c_9>0$ and ∃${\varsigma }_0>0$, obeying $S\left({\varsigma }_0\right)<0$, then Equation (11) admits at least two couples of purely imaginary roots.

Proof. 

(i)

In view of (21), we obtain the following:

$$\begin{array}{ccc}& & \frac{dQ\left(\varphi \right)}{d\varphi }=8q{c}_1{\varphi }^{8q-1}+7q{c}_2{\varphi }^{7q-1}+6q{c}_3{\varphi }^{6q-1}+5q{c}_4{\varphi }^{5q-1}\hfill \\ & & +4q{c}_5{\varphi }^{4q-1}+3q{c}_6{\varphi }^{3q-1}+2q{c}_7{\varphi }^{2q-1}+q{c}_8{\varphi }^{q-1}.\hfill \end{array}$$

(23)

Since $c_l>0(l=,2,\cdots ,8)$, one obtains $\frac{dQ\left(\varphi \right)}{d\varphi }>0$,∀ $\varphi >0$. Apart from $Q\left(0\right)=c_9>0$, one knows that Equation (19) admits no positive real root. In view of $a_4+a_7\ne 0$, one understands that $s=0$ is not the root of (11). This concludes the proof of (i).

(ii)

Clearly, $S\left(0\right)=c_9>0,S\left({\varsigma }_0\right)<0({\varsigma }_0>0)$ and ${lim}_{\varphi \to +\infty }\frac{S\left(\varphi \right)}{d\varphi }=+\infty$; then, there exist ${\varsigma }_1\in (0,{\varsigma }_0)$ and ${\varsigma }_2\in ({\varsigma }_0,+\infty )$ obeying $S\left({\varsigma }_1\right)=S\left({\varsigma }_2\right)=0,$ and Equation (19) admits at least both positive real roots. So, (11) admits at least two couples of purely imaginary roots. This ends the the proof of (ii).

□

Assume that Equation (19) admits eight positive real roots ${\varphi }_i(i=1,2,\cdots ,8)$. It follows from (14) that

$${\theta }_i^k=\frac{1}{2{\varphi }_i}\left[\arccos \left(-\frac{{\mathcal{A}}_1\left({\varphi }_i\right){\mathcal{A}}_3\left({\varphi }_i\right)+{\mathcal{A}}_2\left({\varphi }_i\right){\mathcal{A}}_4\left({\varphi }_i\right)}{{\mathcal{A}}_1^2\left({\varphi }_i\right)+{\mathcal{A}}_2^2\left({\varphi }_i\right)}\right)+2k\pi \right],$$

(24)

where $k=0,1,2,\cdots ,i=1,2,\cdots ,8.$ Let

$${\theta }_0=\underset{i=1,2,\cdots ,8}{min}\left\{{\theta }_i^0\right\},{\varphi }_0{=\varphi |}_{\theta ={\theta }_0}.$$

(25)

Next, the following assumption is needed:

($A_2)$ $U_{11}{U}_{21}+U_{12}{U}_{22}>0,$ where

$$\left\{\begin{array}{c}{\displaystyle U_{11}=4q{\varphi }_0^{4q-1}\cos \frac{(4q-1)\pi }{2}+3q{a}_1{\varphi }_0^{3q-1}\cos \frac{(3q-1)\pi }{2}}\hfill \\ {\displaystyle +2q{a}_2{\varphi }_0^{2q-1}\cos \frac{(2q-1)\pi }{2}+q{a}_3{\varphi }_0^{q-1}\cos \frac{(q-1)\pi }{2}}\hfill \\ {\displaystyle +\left(2q{a}_5{\theta }_0^{2q-1}\cos \frac{(2q-1)\pi }{2}+q{a}_6{\theta }_0^{q-1}\cos \frac{(q-1)\pi }{2}\right)\cos 2{\varphi }_0{\theta }_0}\hfill \\ {\displaystyle +\left(2q{a}_5{\theta }_0^{2q-1}\sin \frac{(2q-1)\pi }{2}+q{a}_6{\theta }_0^{q-1}\sin \frac{(q-1)\pi }{2}\right)\sin 2{\varphi }_0{\theta }_0,}\hfill \\ {\displaystyle U_{12}=4q{\varphi }_0^{4q-1}\sin \frac{(4q-1)\pi }{2}+3q{a}_1{\varphi }_0^{3q-1}\sin \frac{(3q-1)\pi }{2}}\hfill \\ {\displaystyle +2q{a}_2{\varphi }_0^{2q-1}\sin \frac{(2q-1)\pi }{2}+q{a}_3{\varphi }_0^{q-1}\sin \frac{(q-1)\pi }{2}}\hfill \\ {\displaystyle -\left(2q{a}_5{\varphi }_0^{2q-1}\cos \frac{(2q-1)\pi }{2}+q{a}_6{\varphi }_0^{q-1}\cos \frac{(q-1)\pi }{2}\right)\sin 2{\varphi }_0{\theta }_0}\hfill \\ {\displaystyle +\left(2q{a}_5{\varphi }_0^{2q-1}\sin \frac{(2q-1)\pi }{2}+q{a}_6{\varphi }_0^{q-1}\sin \frac{(q-1)\pi }{2}\right)\cos 2{\varphi }_0{\theta }_0,}\hfill \\ {\displaystyle U_{21}=2{\varphi }_0\left(a_5{\varphi }_0^{2q}\cos q\pi +a_6{\varphi }_0^q\cos \frac{q\pi }{2}+a_7\right)\sin 2{\varphi }_0{\theta }_0}\hfill \\ {\displaystyle +2{\varphi }_0\left(a_5{\varphi }_0^{2q}\sin q\pi +a_6{\varphi }_0^q\sin \frac{q\pi }{2}\right)\cos 2{\varphi }_0{\theta }_0,}\hfill \\ {\displaystyle U_{22}=-2{\varphi }_0\left(a_5{\varphi }_0^{2q}\cos q\pi +a_6{\varphi }_0^q\cos \frac{q\pi }{2}+a_7\right)\cos 2{\varphi }_0{\theta }_0}\hfill \\ {\displaystyle +2{\varphi }_0\left(a_5{\varphi }_0^{2q}\sin q\pi +a_6{\varphi }_0^q\sin \frac{q\pi }{2}\right)\sin 2{\varphi }_0{\theta }_0.}\hfill \end{array}\right.$$

(26)

Lemma 4. 

Let $s\left(\theta \right)={\zeta }_1\left(\theta \right)+i{\zeta }_2\left(\theta \right)$ be the root of Equation (11) near $\theta ={\theta }_0$, obeying ${\zeta }_1\left({\theta }_0\right)=0,{\zeta }_2\left({\theta }_0\right)={\varphi }_0$, then $Re\left[\frac{ds}{d\theta }\right]{|}_{\theta ={\theta }_0,\varphi ={\varphi }_0}>0.$

Proof. 

Due to (11), one obtains the following:

$$\begin{array}{ccc}& & \left(4q{s}^{4q-1}+3q{a}_1{s}^{3q-1}+2q{a}_2{s}^{2q-1}+q{a}_3{s}^{q-1}\right)\frac{ds}{d\theta }\hfill \\ & & +\left(2q{a}_5{s}^{2q-1}+q{a}_6{s}^{q-1}\right)\frac{ds}{d\theta }{e}^{-2s\theta }\hfill \\ & & -2e^{-2s\theta }\left(\frac{ds}{d\theta }\theta +s\right)(a_5{s}^{2q}+a_6{s}^q+a_7)=0.\hfill \end{array}$$

(27)

Using (27), we obtain the following:

$${\left[\frac{ds}{d\theta }\right]}^{-1}=\frac{U_1\left(s\right)}{U_2\left(s\right)}-\frac{\theta }{s},$$

(28)

where

$$\left\{\begin{array}{c}{\displaystyle U_1\left(s\right)=\left(4q{s}^{4q-1}+3q{a}_1{s}^{3q-1}+2q{a}_2{s}^{2q-1}+q{a}_3{s}^{q-1}\right)}\hfill \\ {\displaystyle +\left(2q{a}_5{s}^{2q-1}+q{a}_6{s}^{q-1}\right)e^{-2s\theta },}\hfill \\ {\displaystyle U_2\left(s\right)=2s{e}^{-2s\theta }\left(a_5{s}^{2q}+a_6{s}^q+a_7\right).}\hfill \end{array}\right.$$

(29)

Then,

$$\mathrm{Re}\left\{\frac{ds}{d\theta }\right\}{|}_{\theta ={\theta }_0,\varphi ={\varphi }_0}=\mathrm{Re}\left\{\frac{U_1\left(s\right)}{U_2\left(s\right)}\right\}{|}_{\theta ={\theta }_0,\varphi ={\varphi }_0}=\frac{U_{11}{U}_{21}+U_{12}{U}_{22}}{U_{21}^2+U_{22}^2}.$$

(30)

By $\left(A_2\right)$, we have the following:

$$\mathrm{Re}\left\{{\left[\frac{ds}{d\theta }\right]}^{-1}\right\}{|}_{\theta ={\theta }_0,\varphi ={\varphi }_0}>0.$$

(31)

□

Next, we provide the following assumption:

($A_3)$ The following inequalities hold:

$$\left\{\begin{array}{c}{\displaystyle {\Pi }_1=a_4>0,}\hfill \\ {\displaystyle {\Pi }_2=det\left[\begin{array}{cc}{a}_4& 1\\ a_3+a_6& a_2+a_5\end{array}\right]>0,}\hfill \\ {\displaystyle {\Pi }_3=det\left[\begin{array}{ccc}{a}_4& 1& 0\\ a_3+a_6& a_2+a_5& a_4\\ 0& a_4+a_7& a_3+a_6\end{array}\right]>0,}\hfill \\ {\displaystyle {\Pi }_4=(a_3+a_6){\Pi }_3>0.}\hfill \end{array}\right.$$

(32)

Lemma 5. 

If $\theta =0$ and $\left(A_3\right)$ is satisfied, then system (3) remains locally asymptotically stable.

Proof. 

When $\theta =0$, then (11) becomes the following:

$${\lambda }^4+a_1{\lambda }^3+(a_2+a_5){\lambda }^2+(a_3+a_6)\lambda +a_4+a_7=0.$$

(33)

Due to $\left(A_3\right)$, one knows that every root ${\lambda }_i$ of (33) obeys $|\arg \left({\lambda }_i\right)|>\frac{q\pi }{2}(i=1,2,3,4).$ Thus, Lemma 5 is correct. □

Related to the discussion above, the following outcome is lightly acquired as follows:

Theorem 1. 

If $\left(A_1\right)$–$\left(A_3\right)$ are fulfilled, then the positive equilibrium point $W_0(w_{1\ast },w_{2\ast },w_{3\ast },w_{4\ast })$ of system (3) remains locally asymptotically stable if θ falls into the range of $[0,{\theta }_0)$. Moreover, system (3) generates a Hopf bifurcation around the positive equilibrium point $W_0(w_{1\ast },w_{2\ast },w_{3\ast },w_{4\ast })$ when $\theta ={\theta }_0.$

4. Bifurcation Control via the ${\mathit{PD}}^{\mathit{p}}$ Controller

In this segment, we handle the Hopf bifurcation control aspect of system (3) via the $P{D}^p$ controller. In view of the idea in [56,57], we design the following $P{D}^p$ controllers:

$$\rho \left(t\right)={\rho }_p(w_2\left(t\right)-w_{2\ast })+{\rho }_d\frac{d^q(w_2\left(t\right)-w_{2\ast })}{d{t}^q},$$

(34)

and

$$\varrho \left(t\right)={\rho }_p(w_4\left(t\right)-w_{4\ast })+{\rho }_d\frac{d^q(w_4\left(t\right)-w_{4\ast })}{d{t}^q},$$

(35)

where ${\rho }_d\ne 1$ and ${\rho }_p$ stand for the derivative control parameter and the proportional control parameter, respectively. Adding the two controllers, $\rho \left(t\right)$ and $\varrho \left(t\right)$, to the second equation and the fourth equation of system (3), respectively, we have the following:

$$\left\{\begin{array}{c}{\displaystyle \frac{d^q{w}_1\left(t\right)}{d{t}^q}=a\left(w_2+w_1-\frac{1}{3}{w}_1^3\right)+b(w_1-w_3(t-\theta )),}\hfill \\ {\displaystyle \frac{d^q{w}_2\left(t\right)}{d{t}^q}=-\frac{1}{a}(w_1-\alpha +\beta w_2)+{\rho }_p(w_2\left(t\right)-w_{2\ast })+{\rho }_d\frac{d^q(w_2\left(t\right)-w_{2\ast })}{d{t}^q},}\hfill \\ {\displaystyle \frac{d^q{w}_3\left(t\right)}{d{t}^q}=a\left(w_4+w_3-\frac{1}{3}{w}_3^3\right)+b(w_3-w_1(t-\theta )),}\hfill \\ {\displaystyle \frac{d^q{w}_4\left(t\right)}{d{t}^q}=-\frac{1}{a}(w_3-\alpha +\beta w_4)+{\rho }_p(w_4\left(t\right)-w_{4\ast })+{\rho }_d\frac{d^q(w_4\left(t\right)-w_{4\ast })}{d{t}^q}.}\hfill \end{array}\right.$$

(36)

Clearly, system (36) has the same equilibrium point as that of system (3). The linear system of Equation (36) near $W_0(w_{1\ast },w_{2\ast },w_{3\ast },w_{4\ast })$ can be expressed as follows:

$$\left\{\begin{array}{c}{\displaystyle \frac{d^q{w}_1\left(t\right)}{d{t}^q}=(a+b-a{w}_{1\ast }^2)w_1\left(t\right)+a{w}_2\left(t\right)-b{w}_3(t-\theta ),}\hfill \\ {\displaystyle \frac{d^q{w}_2\left(t\right)}{d{t}^q}=-\frac{1}{a(1-{\rho }_d)}{w}_1\left(t\right)-\frac{1}{1-{\rho }_d}\left(\frac{\beta }{a}-{\rho }_p\right)w_2\left(t\right),}\hfill \\ {\displaystyle \frac{d^q{w}_3\left(t\right)}{d{t}^q}=(a+b-a{w}_{3\ast }^2)w_3\left(t\right)+a{w}_4\left(t\right)-b{w}_1(t-\theta ),}\hfill \\ {\displaystyle \frac{d^q{w}_4\left(t\right)}{d{t}^q}=-\frac{1}{a(1-{\rho }_d)}{w}_3\left(t\right)-\frac{1}{1-{\rho }_d}\left(\frac{\beta }{a}-{\rho }_p\right)w_4\left(t\right).}\hfill \end{array}\right.$$

(37)

The associated characteristic equation of (37) reads as follows:

$$det\left[\begin{array}{cccc}{s}^q-(a+b-a{w}_{1\ast }^2)& -a& b{e}^{-s\theta }& 0\\ \frac{1}{a(1-{\rho }_d)}& s^q+\frac{\frac{\beta }{a}-{\rho }_p}{1-{\rho }_d}& 0& 0\\ b{e}^{-s\theta }& 0& s^q-(a+b-a{w}_{3\ast }^2)& -a\\ 0& 0& \frac{1}{a(1-{\rho }_d)}& s^q+\frac{\frac{\beta }{a}-{\rho }_p}{1-{\rho }_d}\end{array}\right]=0.$$

(38)

Applying (38), one obtains the following:

$$s^{4q}+d_1{s}^{3q}+d_2{s}^{2q}+d_3{s}^q+d_4+(d_5{s}^{2q}+d_6{s}^q+d_7)e^{-2s\theta }=0,$$

(39)

where

$$\left\{\begin{array}{c}{\displaystyle d_1=\frac{2}{a(1-{\rho }_d)}-2a-2b+a(w_{1\ast }^2+w_{3\ast }^2),}\hfill \\ {\displaystyle d_2={\left(\frac{\frac{\beta }{a}-{\rho }_p}{1-{\rho }_d}\right)}^2-2\left(\frac{\frac{\beta }{a}-{\rho }_p}{1-{\rho }_d}\right)(a+b-a{w}_{1\ast }^2)}\hfill \\ {\displaystyle -2\left(\frac{\frac{\beta }{a}-{\rho }_p}{1-{\rho }_d}\right)(a+b-a{w}_{3\ast }^2)}\hfill \\ {\displaystyle +(a+b-a{w}_{1\ast }^2)(a+b-a{w}_{3\ast }^2)+2,}\hfill \\ {\displaystyle d_3=\frac{2}{a(1-{\rho }_d)}-(a+b-a{w}_{1\ast }^2)}\hfill \\ {\displaystyle -{\left(\frac{\frac{\beta }{a}-{\rho }_p}{1-{\rho }_d}\right)}^2(a+b-a{w}_{1\ast }^2)}\hfill \\ {\displaystyle +\frac{1}{a(1-{\rho }_d)}(a+b-a{w}_{1\ast }^2)(a+b-a{w}_{3\ast }^2)}\hfill \\ {\displaystyle -{\left(\frac{\frac{\beta }{a}-{\rho }_p}{1-{\rho }_d}\right)}^2(a+b-a{w}_{3\ast }^2)}\hfill \\ {\displaystyle +\frac{1}{a(1-{\rho }_d)}(a+b-a{w}_{1\ast }^2)(a+b-a{w}_{3\ast }^2)}\hfill \\ {\displaystyle -(a+b-a{w}_{3\ast }^2),}\hfill \end{array}\right.$$

(40)

and

$$\left\{\begin{array}{c}{\displaystyle d_4=1-\frac{1}{a(1-{\rho }_d)}(a+b-a{w}_{1\ast }^2)+{\left(\frac{\frac{\beta }{a}-{\rho }_p}{1-{\rho }_d}\right)}^2}\hfill \\ {\displaystyle \times (a+b-a{w}_{1\ast }^2)(a+b-a{w}_{3\ast }^2)}\hfill \\ {\displaystyle -\frac{1}{a(1-{\rho }_d)}(a+b-a{w}_{3\ast }^2),}\hfill \\ {\displaystyle d_5=-b^2,}\hfill \\ {\displaystyle d_6=\frac{2b^2}{a(1-{\rho }_d)},}\hfill \\ {\displaystyle d_7={\left(\frac{\frac{\beta }{a}-{\rho }_p}{1-{\rho }_d}\right)}^2.}\hfill \end{array}\right.$$

(41)

Assuming that $s=i\phi =\phi \left(\cos \frac{\pi }{2}+i\sin \frac{\pi }{2}\right)$ is a root of (39), we obtain the following:

$$\begin{array}{ccc}& & {\phi }^{4q}(\cos 2q\pi +i\sin 2q\pi )+d_1{\phi }^{3q}\left(\cos \frac{3q\pi }{2}+i\sin \frac{3q\pi }{2}\right)\hfill \\ & & +d_2{\phi }^{2q}(\cos q\pi +i\sin q\pi )+d_3{\phi }^q\left(\cos \frac{q\pi }{2}+i\sin \frac{q\pi }{2}\right)\hfill \\ & & +d_4+\left[d_5{\phi }^{2q}(\cos q\pi +i\sin q\pi )+d_6{\phi }^q\left(\cos \frac{q\pi }{2}+i\sin \frac{q\pi }{2}\right)+d_7\right]\hfill \\ & & \times (\cos 2\varphi \theta -i\sin 2\varphi \theta )=0,\hfill \end{array}$$

(42)

which leads to

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{B}}_1\left(\phi \right)\cos 2\phi \theta +{\mathcal{B}}_2\left(\phi \right)\cos 2\phi \theta =-{\mathcal{B}}_3\left(\phi \right),}\hfill \\ {\displaystyle {\mathcal{B}}_2\left(\phi \right)\cos 2\phi \theta -{\mathcal{B}}_1\left(\phi \right)\sin 2\phi \theta =-{\mathcal{B}}_4\left(\phi \right),}\hfill \end{array}\right.$$

(43)

where

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{B}}_1\left(\phi \right)=e_1{\phi }^{2q}+e_2{\phi }^q+e_3,}\hfill \\ {\displaystyle {\mathcal{B}}_1\left(\phi \right)=e_4{\phi }^{2q}+e_5{\phi }^q,}\hfill \\ {\displaystyle {\mathcal{B}}_1\left(\phi \right)=e_6{\phi }^{4q}+e_7{\phi }^{3q}+e_8{\phi }^{2q}+e_9{\phi }^q+e_{10},}\hfill \\ {\displaystyle {\mathcal{B}}_1\left(\phi \right)=e_{11}{\phi }^{4q}+e_{12}{\phi }^{3q}+e_{13}{\phi }^{2q}+e_{14}{\phi }^q,}\hfill \end{array}\right.$$

(44)

where

$$\left\{\begin{array}{c}{\displaystyle e_1=d_5\cos q\pi ,}\hfill \\ {\displaystyle e_2=d_6\cos \frac{2q\pi }{2},}\hfill \\ {\displaystyle e_3=d_7,}\hfill \\ {\displaystyle e_4=d_5\sin q\pi ,}\hfill \\ {\displaystyle e_5=d_6\sin \frac{2q\pi }{2},}\hfill \\ {\displaystyle e_6=\cos 2q\pi ,}\hfill \\ {\displaystyle e_7=d_1\cos \frac{3q\pi }{2},}\hfill \\ {\displaystyle e_8=d_2\cos q\pi ,}\hfill \\ {\displaystyle e_9=d_3\cos \frac{2q\pi }{2},}\hfill \\ {\displaystyle e_{10}=d_4,}\hfill \\ {\displaystyle e_{11}=\sin 2q\pi ,}\hfill \\ {\displaystyle e_{12}=d_1\sin \frac{3q\pi }{2},}\hfill \\ {\displaystyle e_{13}=d_2\sin q\pi ,}\hfill \\ {\displaystyle e_{14}=d_3\sin \frac{2q\pi }{2}.}\hfill \end{array}\right.$$

(45)

It follows from (43) that

$${\mathcal{B}}_1^2\left(\phi \right)+{\mathcal{B}}_2^2\left(\phi \right)={\mathcal{B}}_3^2\left(\phi \right)+{\mathcal{B}}_4^2\left(\phi \right),$$

(46)

which results in

$$f_1{\varphi }^{8q}+f_2{\phi }^{7q}+f_3{\phi }^{6q}+f_4{\phi }^{5q}+f_5{\phi }^{4q}+f_6{\phi }^{3q}+f_7{\phi }^{2q}+f_8{\phi }^q+f_9=0,$$

(47)

where

$$\left\{\begin{array}{c}{\displaystyle f_1=e_6^2+e_{11}^2,}\hfill \\ {\displaystyle f_2=2(e_6{e}_7+e_{11}{e}_{12}),}\hfill \\ {\displaystyle f_3=e_7^2+e_{12}^2+2(e_6{e}_8+e_{11}{e}_{13}),}\hfill \\ {\displaystyle f_4=2(e_6{e}_9+e_7{e}_8+e_{11}{e}_{14}+e_{12}{e}_{13}),}\hfill \\ {\displaystyle f_5=e_8^2+e_{13}^2-e_1^2-e_4^2+2(e_6{e}_{10}+e_7{e}_9+e_{12}{e}_{13}),}\hfill \\ {\displaystyle f_6=2(e_7{e}_{10}+e_8{e}_9+e_{13}{e}_{14}-2e_1{e}_2-2e_4{e}_5),}\hfill \\ {\displaystyle f_7=e_9^2+e_{14}^2-e_5^2-e_2^2+2(e_8{e}_{10}-e_1{e}_3),}\hfill \\ {\displaystyle f_8=2(e_9{e}_{10}-e_2{e}_3),}\hfill \\ {\displaystyle f_9=e_{10}^2-e_3^2.}\hfill \end{array}\right.$$

(48)

Denote

$$R\left(\phi \right)=f_1{\varphi }^{8q}+f_2{\phi }^{7q}+f_3{\phi }^{6q}+f_4{\phi }^{5q}+f_5{\phi }^{4q}+f_6{\phi }^{3q}+f_7{\phi }^{2q}+f_8{\phi }^q+f_9$$

(49)

and

$$T\left(\phi \right)=f_1{\phi }^8+f_2{\phi }^7+f_3{\phi }^6+f_4{\phi }^5+f_5{\phi }^4+f_6{\phi }^3+f_7{\phi }^2+f_8\phi +f_9$$

(50)

Lemma 6. 

(i)

Suppose that $d_4+d_7\ne 0$ and $f_k>0(k=1,2,\cdots ,9)$, then Equation (36) does not have any roots with a zero real part.

(ii)

Suppose that $f_9>0$ and ∃${\xi }_0>0$ obey $T\left({\xi }_0\right)<0$, then Equation (36) admits at least two couples of purely imaginary roots.

Proof. 

(i)

In view of (49), we obtain the following:

$$\begin{array}{ccc}& & \frac{dR\left(\phi \right)}{d\phi }=8q{f}_1{\phi }^{8q-1}+7q{f}_2{\phi }^{7q-1}+6q{f}_3{\phi }^{6q-1}+5q{f}_4{\phi }^{5q-1}\hfill \\ & & +4q{f}_5{\phi }^{4q-1}+3q{f}_6{\phi }^{3q-1}+2q{f}_7{\phi }^{2q-1}+q{f}_8{\phi }^{q-1}.\hfill \end{array}$$

(51)

Since $f_l>0(l=,2,\cdots ,8)$, one obtains $\frac{dR\left(\phi \right)}{d\phi }>0$,∀ $\phi >0$. Apart from $R\left(0\right)=f_9>0$, one knows that Equation (47) admits no positive real root. In view of $d_4+d_7\ne 0$, one understands that $s=0$ is not the root of (36). This concludes the proof of (i).

(ii)

Clearly, $T\left(0\right)=f_9>0,T\left({\xi }_0\right)<0({\xi }_0>0)$ and ${lim}_{\phi \to +\infty }\frac{T\left(\phi \right)}{d\phi }=+\infty$, then there exist ${\xi }_1\in (0,{\xi }_0)$ and ${\xi }_2\in ({\xi }_0,+\infty )$, obeying $T\left({\xi }_1\right)=T\left({\xi }_2\right)=0,$ then Equation (47) admits at least both positive real roots. So, (36) admits at least two couples of purely imaginary roots. This concludes the proof of (ii).

□

Assume that Equation (47) admits eight positive real roots ${\phi }_i(i=1,2,\cdots ,8)$. It follows from (43) that

$${\theta }_j^l=\frac{1}{2{\phi }_j}\left[\arccos \left(-\frac{{\mathcal{B}}_1\left({\phi }_j\right){\mathcal{B}}_3\left({\phi }_j\right)+{\mathcal{B}}_2\left({\phi }_j\right){\mathcal{B}}_4\left({\phi }_j\right)}{{\mathcal{B}}_1^2\left({\phi }_j\right)+{\mathcal{B}}_2^2\left({\phi }_j\right)}\right)+2l\pi \right],$$

(52)

where $l=0,1,2,\cdots ,j=1,2,\cdots ,8.$ Let

$${\theta }_{\ast }=\underset{j=1,2,\cdots ,8}{min}\left\{{\theta }_j^0\right\},{\phi }_0{=\phi |}_{\theta ={\theta }_{\ast }}.$$

(53)

Next, the following assumption is needed:

($A_4)$ $V_{11}{V}_{21}+V_{12}{V}_{22}>0,$ where

$$\left\{\begin{array}{c}{\displaystyle V_{11}=4q{\phi }_0^{4q-1}\cos \frac{(4q-1)\pi }{2}+3q{d}_1{\phi }_0^{3q-1}\cos \frac{(3q-1)\pi }{2}}\hfill \\ {\displaystyle +2q{d}_2{\phi }_0^{2q-1}\cos \frac{(2q-1)\pi }{2}+q{d}_3{\phi }_0^{q-1}\cos \frac{(q-1)\pi }{2}}\hfill \\ {\displaystyle +\left(2q{d}_5{\phi }_0^{2q-1}\cos \frac{(2q-1)\pi }{2}+q{d}_6{\phi }_0^{q-1}\cos \frac{(q-1)\pi }{2}\right)\cos 2{\phi }_0{\theta }_{\ast }}\hfill \\ {\displaystyle +\left(2q{d}_5{\phi }_0^{2q-1}\sin \frac{(2q-1)\pi }{2}+q{d}_6{\phi }_0^{q-1}\sin \frac{(q-1)\pi }{2}\right)\sin 2{\phi }_0{\theta }_{\ast },}\hfill \\ {\displaystyle V_{12}=4q{\phi }_0^{4q-1}\sin \frac{(4q-1)\pi }{2}+3q{d}_1{\phi }_0^{3q-1}\sin \frac{(3q-1)\pi }{2}}\hfill \\ {\displaystyle +2q{d}_2{\phi }_0^{2q-1}\sin \frac{(2q-1)\pi }{2}+q{d}_3{\phi }_0^{q-1}\sin \frac{(q-1)\pi }{2}}\hfill \\ {\displaystyle -\left(2q{d}_5{\phi }_0^{2q-1}\cos \frac{(2q-1)\pi }{2}+q{d}_6{\phi }_0^{q-1}\cos \frac{(q-1)\pi }{2}\right)\sin 2{\phi }_0{\theta }_{\ast }}\hfill \\ {\displaystyle +\left(2q{d}_5{\phi }_0^{2q-1}\sin \frac{(2q-1)\pi }{2}+q{d}_6{\phi }_0^{q-1}\sin \frac{(q-1)\pi }{2}\right)\cos 2{\phi }_0{\theta }_{\ast },}\hfill \\ {\displaystyle V_{21}=2{\phi }_0\left(d_5{\phi }_0^{2q}\cos q\pi +d_6{\phi }_0^q\cos \frac{q\pi }{2}+d_7\right)\sin 2{\phi }_0{\theta }_{\ast }}\hfill \\ {\displaystyle +2{\phi }_0\left(d_5{\phi }_0^{2q}\sin q\pi +d_6{\phi }_0^q\sin \frac{q\pi }{2}\right)\cos 2{\phi }_0{\theta }_{\ast },}\hfill \\ {\displaystyle V_{22}=-2{\phi }_0\left(d_5{\phi }_0^{2q}\cos q\pi +d_6{\phi }_0^q\cos \frac{q\pi }{2}+d_7\right)\cos 2{\phi }_0{\theta }_{\ast }}\hfill \\ {\displaystyle +2{\phi }_0\left(d_5{\phi }_0^{2q}\sin q\pi +d_6{\phi }_0^q\sin \frac{q\pi }{2}\right)\sin 2{\phi }_0{\theta }_{\ast }.}\hfill \end{array}\right.$$

(54)

Lemma 7. 

Let $s\left(\theta \right)={\tau }_1\left(\theta \right)+i{\tau }_2\left(\theta \right)$ be the root of Equation (39) near $\theta ={\theta }_{\ast }$, obeying ${\tau }_1\left({\theta }_{\ast }\right)=0,{\tau }_2\left({\theta }_{\ast }\right)={\phi }_0$, then $\mathit{Re}\left[\frac{ds}{d\theta }\right]{|}_{\theta ={\theta }_{\ast },\phi ={\phi }_0}>0.$

Proof. 

Due to (39), we obtain the following:

$$\begin{array}{ccc}& & \left(4q{s}^{4q-1}+3q{d}_1{s}^{3q-1}+2q{d}_2{s}^{2q-1}+q{d}_3{s}^{q-1}\right)\frac{ds}{d\theta }\hfill \\ & & +\left(2q{d}_5{s}^{2q-1}+q{d}_6{s}^{q-1}\right)\frac{ds}{d\theta }{e}^{-2s\theta }\hfill \\ & & -2e^{-2s\theta }\left(\frac{ds}{d\theta }\theta +s\right)(d_5{s}^{2q}+d_6{s}^q+d_7)=0.\hfill \end{array}$$

(55)

Using (55), we obtain the following:

$${\left[\frac{ds}{d\theta }\right]}^{-1}=\frac{V_1\left(s\right)}{V_2\left(s\right)}-\frac{\theta }{s},$$

(56)

where

$$\left\{\begin{array}{c}{\displaystyle V_1\left(s\right)=\left(4q{s}^{4q-1}+3q{d}_1{s}^{3q-1}+2q{d}_2{s}^{2q-1}+q{d}_3{s}^{q-1}\right)}\hfill \\ {\displaystyle +\left(2q{d}_5{s}^{2q-1}+q{d}_6{s}^{q-1}\right)e^{-2s\theta },}\hfill \\ {\displaystyle V_2\left(s\right)=2s{e}^{-2s\theta }\left(d_5{s}^{2q}+d_6{s}^q+d_7\right).}\hfill \end{array}\right.$$

(57)

Then,

$$\mathrm{Re}\left\{\frac{ds}{d\theta }\right\}{|}_{\theta ={\theta }_{\ast },\phi ={\phi }_0}=\mathrm{Re}\left\{\frac{V_1\left(s\right)}{V_2\left(s\right)}\right\}{|}_{\theta ={\theta }_{\ast },\phi ={\phi }_0}=\frac{V_{11}{V}_{21}+V_{12}{V}_{22}}{V_{21}^2+V_{22}^2}.$$

(58)

By $\left(A_4\right)$, we have the following:

$$\mathrm{Re}\left\{{\left[\frac{ds}{d\theta }\right]}^{-1}\right\}{|}_{\theta ={\theta }_{\ast },\phi ={\phi }_0}>0.$$

(59)

□

Next, we provide the following assumption:

(A₅) The following inequalities hold:

$$\left\{\begin{array}{c}{\displaystyle {\Xi }_1=d_4>0,}\hfill \\ {\displaystyle {\Xi }_2=det\left[\begin{array}{cc}{d}_4& 1\\ d_3+d_6& d_2+d_5\end{array}\right]>0,}\hfill \\ {\displaystyle {\Xi }_3=det\left[\begin{array}{ccc}{d}_4& 1& 0\\ d_3+d_6& d_2+d_5& d_4\\ 0& d_4+d_7& d_3+d_6\end{array}\right]>0,}\hfill \\ {\displaystyle {\Xi }_4=(d_3+d_6){\Xi }_3>0.}\hfill \end{array}\right.$$

(60)

Lemma 8. 

If $\theta =0$ and $\left(A_5\right)$ are satisfied, then system (36) remains locally asymptotically stable.

Proof. 

When $\theta =0$, then (36) becomes

$${\lambda }^4+d_1{\lambda }^3+(d_2+d_5){\lambda }^2+(d_3+d_6)\lambda +d_4+d_7=0.$$

(61)

By virtue of $\left(A_5\right)$, one knows that every root ${\lambda }_i$ of (61) obeys $|\arg \left({\lambda }_j\right)|>\frac{q\pi }{2}(j=1,2,3,4).$ Thus, Lemma 8 is correct. □

Relying on the discussion above, the following outcome is lightly acquired:

Theorem 2. 

If $\left(A_1\right),\left(A_4\right),\left(A_5\right)$ are fulfilled, then the positive equilibrium point $W_0(w_{1\ast },w_{2\ast },w_{3\ast }$, $w_{4\ast })$ of system (36) remains locally asymptotically stable if θ falls into the range of $[0,{\theta }_{\ast })$ and system (36) will generate a Hopf bifurcation around the positive equilibrium point of $W_0(w_{1\ast },w_{2\ast },w_{3\ast },w_{4\ast })$, when $\theta ={\theta }_{\ast }.$

5. Bifurcation Control via Hybrid Controller

In this segment, we design a proper hybrid controller, including state feedback and parameter perturbation to control the Hopf bifurcation of system (3). Following the idea of [46,50], we obtain the following fractional-order controlled FitzHugh–Nagumo neural model:

$$\left\{\begin{array}{c}{\displaystyle \frac{d^q{w}_1\left(t\right)}{d{t}^q}=a\left(w_2+w_1-\frac{1}{3}{w}_1^3\right)+b(w_1-w_3(t-\theta )),}\hfill \\ {\displaystyle \frac{d^q{w}_2\left(t\right)}{d{t}^q}={\kappa }_1\left[-\frac{1}{a}(w_1-\alpha +\beta w_2)\right]+{\kappa }_2[w_2(t-\theta )-w_2\left(t\right)],}\hfill \\ {\displaystyle \frac{d^q{w}_3\left(t\right)}{d{t}^q}=a\left(w_4+w_3-\frac{1}{3}{w}_3^3\right)+b(w_3-w_1(t-\theta )),}\hfill \\ {\displaystyle \frac{d^q{w}_4\left(t\right)}{d{t}^q}={\kappa }_1\left[-\frac{1}{a}(w_3-\alpha +\beta w_4)\right]+{\kappa }_2[w_4(t-\theta )-w_4\left(t\right)],}\hfill \end{array}\right.$$

(62)

where ${\kappa }_1,{\kappa }_2$ are feedback gain parameters. Clearly, system (62) owns the same equilibrium point as that of system (3). The linear system of Equation (62) near $W_0(w_{1\ast },w_{2\ast },w_{3\ast },w_{4\ast })$ can be expressed as follows:

$$\left\{\begin{array}{c}{\displaystyle \frac{d^q{w}_1\left(t\right)}{d{t}^q}=(a+b-a{w}_{1\ast }^2)w_1\left(t\right)+a{w}_2\left(t\right)-b{w}_3(t-\theta ),}\hfill \\ {\displaystyle \frac{d^q{w}_2\left(t\right)}{d{t}^q}=-\frac{{\kappa }_1}{a}{w}_1\left(t\right)-\left(\frac{\beta {\kappa }_1}{a}+{\kappa }_2\right)w_2\left(t\right)+{\kappa }_2{w}_2(t-\theta ),}\hfill \\ {\displaystyle \frac{d^q{w}_3\left(t\right)}{d{t}^q}=(a+b-a{w}_{3\ast }^2)w_3\left(t\right)+a{w}_4\left(t\right)-b{w}_1(t-\theta ),}\hfill \\ {\displaystyle \frac{d^q{w}_4\left(t\right)}{d{t}^q}=-\frac{{\kappa }_1}{a}{w}_3\left(t\right)-\left(\frac{\beta {\kappa }_1}{a}+{\kappa }_2\right)w_4\left(t\right)+{\kappa }_2{w}_4(t-\theta ).}\hfill \end{array}\right.$$

(63)

The associated characteristic equation of (63) reads as follows:

$$det\left[\begin{array}{cccc}{s}^q-{\alpha }^{\ast }& -a& b{e}^{-s\theta }& 0\\ {\gamma }^{\ast }& s^q+{\beta }^{\ast }-{\kappa }_2{e}^{-s\theta }& 0& 0\\ b{e}^{-s\theta }& 0& s^q-{\delta }^{\ast }& -a\\ 0& 0& {\gamma }^{\ast }& s^q+{\beta }^{\ast }-{\kappa }_2{e}^{-s\theta }\end{array}\right]=0,$$

(64)

where

$$\left\{\begin{array}{c}{\displaystyle {\alpha }^{\ast }=a+b-a{w}_{1\ast }^2,}\hfill \\ {\displaystyle {\beta }^{\ast }=\frac{\beta {\kappa }_1}{a}+{\kappa }_2}\hfill \\ {\displaystyle {\gamma }^{\ast }=\frac{{\kappa }_1}{a},}\hfill \\ {\displaystyle {\delta }^{\ast }=a+b-a{w}_{3\ast }^2.}\hfill \end{array}\right.$$

(65)

Applying (64), one obtains the following:

$$\begin{array}{ccc}& & s^{4q}+{ϵ}_1{s}^{3q}+{ϵ}_2{s}^{2q}+{ϵ}_3{s}^q+{ϵ}_4+({ϵ}_5{s}^{3q}+{ϵ}_6{s}^{2q}+{ϵ}_7{s}^q+{ϵ}_8)e^{-s\theta }\hfill \\ & & +({ϵ}_9{s}^{2q}+{ϵ}_{10}{s}^q+{ϵ}_{11})e^{-2s\theta }+{ϵ}_{12}{e}^{-3s\theta }=0,\hfill \end{array}$$

(66)

where

$$\left\{\begin{array}{c}{\displaystyle {ϵ}_1=2({\beta }^{\ast }-{\alpha }^{\ast }),}\hfill \\ {\displaystyle {ϵ}_2={\delta }^{\ast }({\alpha }^{\ast }-{\beta }^{\ast })-{\alpha }^{\ast }{\beta }^{\ast }+{\beta }^{\ast }({\beta }^{\ast }-2{\alpha }^{\ast }),}\hfill \\ {\displaystyle {ϵ}_3={\alpha }^{\ast }{\beta }^{\ast }{\delta }^{\ast }+{\beta }^{\ast }[{\delta }^{\ast }({\alpha }^{\ast }-{\beta }^{\ast })-{\alpha }^{\ast }{\beta }^{\ast }]+a{\gamma }^{\ast },}\hfill \\ {\displaystyle {ϵ}_4={\alpha }^{\ast }{\left({\beta }^{\ast }\right)}^2{\delta }^{\ast }-a{\gamma }^{\ast }{\delta }^{\ast },}\hfill \\ {\displaystyle {ϵ}_5=-2{\kappa }_2,}\hfill \\ {\displaystyle {ϵ}_6={\alpha }^{\ast }{\kappa }_2+{\kappa }_2{\delta }^{\ast }-{\beta }^{\ast }{\kappa }_2-{\kappa }_2({\beta }^{\ast }-2{\alpha }^{\ast }),}\hfill \\ {\displaystyle {ϵ}_7={\beta }^{\ast }({\alpha }^{\ast }{\kappa }_2+{\kappa }_2{\delta }^{\ast })-{\alpha }^{\ast }{\kappa }_2{\delta }^{\ast }-{\kappa }_2[{\delta }^{\ast }({\alpha }^{\ast }-{\beta }^{\ast })-{\alpha }^{\ast }{\beta }^{\ast }],}\hfill \\ {\displaystyle {ϵ}_8=-2{\kappa }_2{\beta }^{\ast }{\alpha }^{\ast }{\delta }^{\ast },}\hfill \\ {\displaystyle {ϵ}_9={\kappa }_2,}\hfill \\ {\displaystyle {ϵ}_{10}=-[b^2+{\kappa }_2^2({\alpha }^2{\kappa }_2+{\kappa }_2{\delta }^{\ast })],}\hfill \\ {\displaystyle {ϵ}_{11}={\alpha }^{\ast }{\kappa }_2^2{\delta }^{\ast }-b^2{\beta }^{\ast },}\hfill \\ {\displaystyle {ϵ}_{12}=b^2{\kappa }_1.}\hfill \end{array}\right.$$

(67)

In view of (66), we have the following:

$$\begin{array}{ccc}& & (s^{4q}+{ϵ}_1{s}^{3q}+{ϵ}_2{s}^{2q}+{ϵ}_3{s}^q+{ϵ}_4)e^{s\theta }+({ϵ}_5{s}^{3q}+{ϵ}_6{s}^{2q}+{ϵ}_7{s}^q+{ϵ}_8)\hfill \\ & & +({ϵ}_9{s}^{2q}+{ϵ}_{10}{s}^q+{ϵ}_{11})e^{-s\theta }+{ϵ}_{12}{e}^{-2s\theta }=0,\hfill \end{array}$$

(68)

Assuming that $s=i\vartheta =\vartheta \left(\cos \frac{\pi }{2}+i\sin \frac{\pi }{2}\right)$ is a root of (68), we obtain

$$\begin{array}{ccc}& & [{\vartheta }^{4q}(\cos 2q\pi +i\sin 2q\pi )+{ϵ}_1{\vartheta }^{3q}\left(\cos \frac{3q\pi }{2}+i\sin \frac{3q\pi }{2}\right)\hfill \\ & & +{ϵ}_2{\vartheta }^{2q}(\cos q\pi +i\sin q\pi )+{ϵ}_3{\vartheta }^q\left(\cos \frac{q\pi }{2}+i\sin \frac{q\pi }{2}\right)\hfill \\ & & +{ϵ}_4](\cos \vartheta \theta +i\sin \vartheta \theta )+{ϵ}_5{\vartheta }^{3q}\left(\cos \frac{3q\pi }{2}+i\sin \frac{3q\pi }{2}\right)\hfill \\ & & +{ϵ}_6{\vartheta }^{2q}(\cos q\pi +i\sin q\pi )+{ϵ}_7{\vartheta }^q\left(\cos \frac{q\pi }{2}+i\sin \frac{q\pi }{2}\right)+{ϵ}_8\hfill \\ & & +[{ϵ}_9{\vartheta }^{2q}(\cos q\pi +i\sin q\pi )+{ϵ}_{10}{\vartheta }^q\left(\cos \frac{q\pi }{2}+i\sin \frac{q\pi }{2}\right)+{ϵ}_{11}]\hfill \\ & & \times (\cos \vartheta \theta -i\sin \vartheta \theta )+{ϵ}_{12}(\cos 2\vartheta \theta -i\sin 2\vartheta \theta )=0,\hfill \end{array}$$

(69)

which leads to

$$\left\{\begin{array}{c}{\displaystyle C_1\cos \vartheta \theta +C_2\sin \vartheta \theta +{ϵ}_{12}\cos 2\vartheta \theta =C_3,}\hfill \\ {\displaystyle C_4\cos \vartheta \theta +C_5\sin \vartheta \theta -{ϵ}_{12}\sin 2\vartheta \theta =C_6,}\hfill \end{array}\right.$$

(70)

where

$$\left\{\begin{array}{c}{\displaystyle C_1={\eta }_1{\vartheta }^{4q}+{\eta }_2{\vartheta }^{3q}+{\eta }_3{\vartheta }^{2q}+{\eta }_4{\vartheta }^q+{\eta }_5,}\hfill \\ {\displaystyle C_2={\eta }_6{\vartheta }^{4q}+{\eta }_7{\vartheta }^{3q}+{\eta }_8{\vartheta }^{2q}+{\eta }_9{\vartheta }^q,}\hfill \\ {\displaystyle C_3={\eta }_{10}{\vartheta }^{3q}+{\eta }_{11}{\vartheta }^{2q}+{\eta }_{12}{\vartheta }^q+{\eta }_{13},}\hfill \\ {\displaystyle C_4={\eta }_{14}{\vartheta }^{4q}+{\eta }_{15}{\vartheta }^{3q}+{\eta }_{16}{\vartheta }^{2q}+{\eta }_{17}{\vartheta }^q,}\hfill \\ {\displaystyle C_5={\eta }_{18}{\vartheta }^{4q}+{\eta }_{19}{\vartheta }^{3q}+{\eta }_{20}{\vartheta }^{2q}+{\eta }_{21}{\vartheta }^q+{\eta }_{22},}\hfill \\ {\displaystyle C_6={\eta }_{23}{\vartheta }^{3q}+{\eta }_{24}{\vartheta }^{2q}+{\eta }_{25}{\vartheta }^q,}\hfill \end{array}\right.$$

(71)

where

$$\left\{\begin{array}{c}{\displaystyle {\eta }_1=\cos 2q\pi ,}\hfill \\ {\displaystyle {\eta }_2={ϵ}_1\cos \frac{3q\pi }{2},}\hfill \\ {\displaystyle {\eta }_3=({ϵ}_2+{ϵ}_9)\cos q\pi ,}\hfill \\ {\displaystyle {\eta }_4=({ϵ}_3+{ϵ}_{10})\cos \frac{q\pi }{2},}\hfill \\ {\displaystyle {\eta }_5={ϵ}_4+{ϵ}_{11},}\hfill \\ {\displaystyle {\eta }_6=-\sin 2q\pi }\hfill \\ {\displaystyle {\eta }_7=-{ϵ}_1\sin \frac{3q\pi }{2},}\hfill \\ {\displaystyle {\eta }_8=-({ϵ}_2-{ϵ}_9)\sin q\pi ,}\hfill \\ {\displaystyle {\eta }_9=-({ϵ}_3-{ϵ}_{10})\sin \frac{q\pi }{2},}\hfill \\ {\displaystyle {\eta }_{10}=-{ϵ}_5\cos \frac{3q\pi }{2},}\hfill \\ {\displaystyle {\eta }_{11}=-{ϵ}_6\cos q\pi ,}\hfill \\ {\displaystyle {\eta }_{12}=-{ϵ}_7\cos \frac{q\pi }{2},}\hfill \\ {\displaystyle {\eta }_{13}=-{ϵ}_8,}\hfill \\ {\displaystyle {\eta }_{14}=\sin 2q\pi ,}\hfill \\ {\displaystyle {\eta }_{15}={ϵ}_1\sin \frac{3q\pi }{2},}\hfill \\ {\displaystyle {\eta }_{16}=({ϵ}_2+{ϵ}_9)\sin q\pi ,}\hfill \\ {\displaystyle {\eta }_{17}=({ϵ}_3+{ϵ}_{10})\sin \frac{q\pi }{2},}\hfill \\ {\displaystyle {\eta }_{18}=\cos 2q\pi ,}\hfill \\ {\displaystyle {\eta }_{19}={ϵ}_1\cos \frac{3q\pi }{2},}\hfill \\ {\displaystyle {\eta }_{20}=({ϵ}_2-{ϵ}_9)\cos q\pi ,}\hfill \\ {\displaystyle {\eta }_{21}=({ϵ}_3-{ϵ}_{10})\cos \frac{q\pi }{2},}\hfill \\ {\displaystyle {\eta }_{22}={ϵ}_4-{ϵ}_{11},}\hfill \\ {\displaystyle {\eta }_{23}=-{ϵ}_5\sin \frac{3q\pi }{2},}\hfill \\ {\displaystyle {\eta }_{24}=-{ϵ}_6\sin q\pi ,}\hfill \\ {\displaystyle {\eta }_{25}=-{ϵ}_7\sin \frac{q\pi }{2}.}\hfill \end{array}\right.$$

(72)

By the first equation of (70), we have the following:

$$C_1\cos \vartheta \theta +C_2\left(\pm \sqrt{1-{\cos }^2\vartheta \theta }\right)+{ϵ}_{12}(2{\cos }^2\vartheta \theta -1)=C_3,$$

(73)

which results in the following:

$${\mu }_1{\cos }^4\vartheta \theta +{\mu }_2{\cos }^3\vartheta \theta +{\mu }_3{\cos }^2\vartheta \theta +{\mu }_4\cos \vartheta \theta +{\mu }_5=0,$$

(74)

which results in the following:

$$\left\{\begin{array}{c}{\displaystyle {\mu }_1=4{ϵ}_{12}^2,}\hfill \\ {\displaystyle {\mu }_2=2{ϵ}_{12}({\eta }_1{\vartheta }^{4q}+{\eta }_2{\vartheta }^{3q}+{\eta }_3{\vartheta }^{2q}+{\eta }_4{\vartheta }^q+{\eta }_5),}\hfill \\ {\displaystyle {\mu }_3=-2{ϵ}_{12}({\eta }_{10}{\vartheta }^{3q}+{\eta }_{11}{\vartheta }^{2q}+{\eta }_{12}{\vartheta }^q+{\eta }_{13})}\hfill \\ {\displaystyle +{({\eta }_6{\vartheta }^{4q}+{\eta }_7{\vartheta }^{3q}+{\eta }_8{\vartheta }^{2q}+{\eta }_9{\vartheta }^q)}^2}\hfill \\ {\displaystyle +{({\eta }_1{\vartheta }^{4q}+{\eta }_2{\vartheta }^{3q}+{\eta }_3{\vartheta }^{2q}+{\eta }_4{\vartheta }^q+{\eta }_5)}^2,}\hfill \\ {\displaystyle {\mu }_4=-2({\eta }_{10}{\vartheta }^{3q}+{\eta }_{11}{\vartheta }^{2q}+{\eta }_{12}{\vartheta }^q+{\eta }_{13})}\hfill \\ {\displaystyle \times ({\eta }_1{\vartheta }^{4q}+{\eta }_2{\vartheta }^{3q}+{\eta }_3{\vartheta }^{2q}+{\eta }_4{\vartheta }^q+{\eta }_5),}\hfill \\ {\displaystyle {\mu }_5=({\eta }_{10}{\vartheta }^{3q}+{\eta }_{11}{\vartheta }^{2q}+{\eta }_{12}{\vartheta }^q+{\eta }_{13}+{ϵ}_{12})}\hfill \\ {\displaystyle -{({\eta }_6{\vartheta }^{4q}+{\eta }_7{\vartheta }^{3q}+{\eta }_8{\vartheta }^{2q}+{\eta }_9{\vartheta }^q)}^2.}\hfill \end{array}\right.$$

(75)

By computer, we can solve $\cos \vartheta \theta$ from (74). Here, we assume the following:

$$\cos \vartheta \theta ={\Psi }_1\left(\vartheta \right),$$

(76)

where ${\Psi }_1\left(\vartheta \right)$ stands for a function with respect to $\vartheta$. In the same way, we can also gain the expression of $\sin \vartheta \theta$. Here, we assume the following:

$$\sin \vartheta \theta ={\Psi }_2\left(\vartheta \right),$$

(77)

where ${\Psi }_2\left(\vartheta \right)$ stands for a function with respect to $\vartheta$. Due to (76) and (77), we obtain the following:

$${\Psi }_1^2\left(\vartheta \right)+{\Psi }_2^2\left(\vartheta \right)=1.$$

(78)

By (78), we can solve the expression of $\vartheta$ (say ${\vartheta }_0$). Making use of (76), and obtain the following:

$${\theta }^{\left(l\right)}=\frac{1}{{\vartheta }_0}[\arccos {\Psi }_1\left(\vartheta \right)+2l\pi ],l=0,1,2,\cdots .$$

(79)

Let

$${\theta }_{\diamond }=\underset{\{l=0,1,2,3,\cdots \}}{min}\left\{{\theta }^{\left(l\right)}\right\}.$$

(80)

Therefore, (66) admits a pair of imaginary roots $\pm i{\vartheta }_0$ when $\theta ={\theta }_{\diamond }.$ Now, we prepare the following assumption:

($A_6)$ $Z_{11}{Z}_{21}+Z_{12}{Z}_{22}>0,$ where

$$\left\{\begin{array}{c}{\displaystyle Z_{11}=4q{\vartheta }_0^{4q-1}\cos \frac{(4q-1)\pi }{2}+3q{ϵ}_1{\vartheta }_0^{3q-1}\cos \frac{(3q-1)\pi }{2}}\hfill \\ {\displaystyle +2q{ϵ}_2{\vartheta }_0^{2q-1}\cos \frac{(2q-1)\pi }{2}+q{ϵ}_3{\vartheta }_0^{q-1}\cos \frac{(q-1)\pi }{2}}\hfill \\ {\displaystyle +\left[3q{ϵ}_5{\vartheta }_0^{3q-1}\cos \frac{(3q-1)\pi }{2}+3q{ϵ}_6{\vartheta }_0^{2q-1}\cos \frac{(2q-1)\pi }{2}\right.}\hfill \\ {\displaystyle \left.+q{ϵ}_7{\vartheta }_0^{q-1}\cos \frac{(q-1)\pi }{2}\right]\cos {\vartheta }_0{\theta }_{\diamond }+\left[3q{ϵ}_5{\vartheta }_0^{3q-1}\sin \frac{(3q-1)\pi }{2}\right.}\hfill \\ {\displaystyle \left.+3q{ϵ}_6{\vartheta }_0^{2q-1}\sin \frac{(2q-1)\pi }{2}+q{ϵ}_7{\vartheta }_0^{q-1}\sin \frac{(q-1)\pi }{2}\right]\sin {\vartheta }_0{\theta }_{\diamond }}\hfill \\ {\displaystyle +\left(2q{ϵ}_9{\vartheta }_0^{2q-1}\cos \frac{(2q-1)\pi }{2}+q{ϵ}_{10}{\vartheta }_0^{q-1}\cos \frac{(q-1)\pi }{2}\right)\cos 2{\vartheta }_0{\theta }_{\diamond }}\hfill \\ {\displaystyle +\left(2q{ϵ}_9{\vartheta }_0^{2q-1}\sin \frac{(2q-1)\pi }{2}+q{ϵ}_{10}{\vartheta }_0^{q-1}\sin \frac{(q-1)\pi }{2}\right)\sin 2{\vartheta }_0{\theta }_{\diamond },}\hfill \\ {\displaystyle Z_{12}=4q{\vartheta }_0^{4q-1}\sin \frac{(4q-1)\pi }{2}+3q{ϵ}_1{\vartheta }_0^{3q-1}\sin \frac{(3q-1)\pi }{2}}\hfill \\ {\displaystyle +2q{ϵ}_2{\vartheta }_0^{2q-1}\sin \frac{(2q-1)\pi }{2}+q{ϵ}_3{\vartheta }_0^{q-1}\sin \frac{(q-1)\pi }{2}}\hfill \\ {\displaystyle -\left[3q{ϵ}_5{\vartheta }_0^{3q-1}\cos \frac{(3q-1)\pi }{2}+3q{ϵ}_6{\vartheta }_0^{2q-1}\cos \frac{(2q-1)\pi }{2}\right.}\hfill \end{array}\right.$$

(81)

and

$$\left\{\begin{array}{c}{\displaystyle \left.+q{ϵ}_7{\vartheta }_0^{q-1}\cos \frac{(q-1)\pi }{2}\right]\sin {\vartheta }_0{\theta }_{\diamond }+\left[3q{ϵ}_5{\vartheta }_0^{3q-1}\sin \frac{(3q-1)\pi }{2}\right.}\hfill \\ {\displaystyle \left.+3q{ϵ}_6{\vartheta }_0^{2q-1}\sin \frac{(2q-1)\pi }{2}+q{ϵ}_7{\vartheta }_0^{q-1}\sin \frac{(q-1)\pi }{2}\right]\cos {\vartheta }_0{\theta }_{\diamond }}\hfill \\ {\displaystyle -\left(2q{ϵ}_9{\vartheta }_0^{2q-1}\cos \frac{(2q-1)\pi }{2}+q{ϵ}_{10}{\vartheta }_0^{q-1}\cos \frac{(q-1)\pi }{2}\right)\sin 2{\vartheta }_0{\theta }_{\diamond }}\hfill \\ {\displaystyle +\left(2q{ϵ}_9{\vartheta }_0^{2q-1}\sin \frac{(2q-1)\pi }{2}+q{ϵ}_{10}{\vartheta }_0^{q-1}\sin \frac{(q-1)\pi }{2}\right)\cos 2{\vartheta }_0{\theta }_{\diamond },}\hfill \\ {\displaystyle Z_{21}=\left({ϵ}_5{\vartheta }_0^{3q}\cos \frac{3q\pi }{2}+{ϵ}_6{s}^{2q}\cos q\pi +{ϵ}_7{\vartheta }_0^q\cos \frac{q\pi }{2}+{ϵ}_8\right){\vartheta }_0\sin {\vartheta }_0{\theta }_{\diamond }}\hfill \\ {\displaystyle -\left({ϵ}_5{\vartheta }_0^{3q}\sin \frac{3q\pi }{2}+{ϵ}_6{s}^{2q}\sin q\pi +{ϵ}_7{\vartheta }_0^q\sin \frac{q\pi }{2}\right){\vartheta }_0\cos {\vartheta }_0{\theta }_{\diamond }}\hfill \\ {\displaystyle +2\left({ϵ}_9{\vartheta }_0^{2q}\cos q\pi +{ϵ}_{10}{\vartheta }_0^q\cos \frac{q\pi }{2}+{ϵ}_{11}\right){\vartheta }_0\sin 2{\vartheta }_0{\theta }_{\diamond }}\hfill \\ {\displaystyle -2\left({ϵ}_9{\vartheta }_0^{2q}\sin q\pi +{ϵ}_{10}{\vartheta }_0^q\sin \frac{q\pi }{2}\right){\vartheta }_0\cos 2{\vartheta }_0{\theta }_{\diamond }+3{\vartheta }_0{ϵ}_{12}\sin 3{\vartheta }_0{\theta }_{\diamond },}\hfill \\ {\displaystyle Z_{22}=\left({ϵ}_5{\vartheta }_0^{3q}\cos \frac{3q\pi }{2}+{ϵ}_6{s}^{2q}\cos q\pi +{ϵ}_7{\vartheta }_0^q\cos \frac{q\pi }{2}+{ϵ}_8\right){\vartheta }_0\cos {\vartheta }_0{\theta }_{\diamond }}\hfill \\ {\displaystyle +\left({ϵ}_5{\vartheta }_0^{3q}\sin \frac{3q\pi }{2}+{ϵ}_6{s}^{2q}\sin q\pi +{ϵ}_7{\vartheta }_0^q\sin \frac{q\pi }{2}\right){\vartheta }_0\sin {\vartheta }_0{\theta }_{\diamond }}\hfill \\ {\displaystyle +2\left({ϵ}_9{\vartheta }_0^{2q}\cos q\pi +{ϵ}_{10}{\vartheta }_0^q\cos \frac{q\pi }{2}+{ϵ}_{11}\right){\vartheta }_0\cos 2{\vartheta }_0{\theta }_{\diamond }}\hfill \\ {\displaystyle +2\left({ϵ}_9{\vartheta }_0^{2q}\sin q\pi +{ϵ}_{10}{\vartheta }_0^q\sin \frac{q\pi }{2}\right){\vartheta }_0\sin 2{\vartheta }_0{\theta }_{\diamond }+3{\vartheta }_0{ϵ}_{12}\sin 3{\vartheta }_0{\theta }_{\diamond }.}\hfill \end{array}\right.$$

(82)

Lemma 9. 

Let $s\left(\theta \right)={\varpi }_1\left(\theta \right)+i{\varpi }_2\left(\theta \right)$ be the root of Equation (66) near $\theta ={\theta }_{\diamond }$, obeying ${\varpi }_1\left({\theta }_{\diamond }\right)=0,{\varpi }_2\left({\theta }_{\diamond }\right)={\vartheta }_0$, then $\mathit{Re}\left[\frac{ds}{d\theta }\right]{|}_{\theta ={\theta }_{\diamond },\vartheta ={\vartheta }_0}>0.$

Proof. 

Due to (66), one obtains the following:

$$\begin{array}{ccc}& & \left(4q{s}^{4q-1}+3q{ϵ}_1{s}^{3q-1}+2q{ϵ}_2{s}^{2q-1}+q{ϵ}_3{s}^{q-1}\right)\frac{ds}{d\theta }\hfill \\ & & +\left(3q{ϵ}_5{s}^{3q-1}+3q{ϵ}_6{s}^{q-1}+q{ϵ}_7{s}^{q-1}\right)\frac{ds}{d\theta }{e}^{-s\theta }\hfill \\ & & -2e^{-s\theta }\left(\frac{ds}{d\theta }\theta +s\right)({ϵ}_5{s}^{3q}+{ϵ}_6{s}^{2q}+{ϵ}_7{s}^q+{ϵ}_8)\hfill \\ & & +\left(2q{ϵ}_9{s}^{2q-1}+q{ϵ}_{10}{s}^{q-1}\right)\frac{ds}{d\theta }{e}^{-2s\theta }\hfill \\ & & -2e^{-2s\theta }\left(\frac{ds}{d\theta }\theta +s\right)\left({ϵ}_9{s}^{2q}+{ϵ}_{10}{s}^q+{ϵ}_{11}\right)\hfill \\ & & -3{ϵ}_{12}{e}^{-3s\theta }\left(\frac{ds}{d\theta }\theta +s\right)=0.\hfill \end{array}$$

(83)

Using (83), we obtain the following:

$${\left[\frac{ds}{d\theta }\right]}^{-1}=\frac{Z_1\left(s\right)}{Z_2\left(s\right)}-\frac{\theta }{s},$$

(84)

where

$$\left\{\begin{array}{c}{\displaystyle Z_1\left(s\right)=4q{s}^{4q-1}+3q{ϵ}_1{s}^{3q-1}+2q{ϵ}_2{s}^{2q-1}+q{ϵ}_3{s}^{q-1}}\hfill \\ {\displaystyle +\left(3q{ϵ}_5{s}^{3q-1}+3q{ϵ}_6{s}^{q-1}+q{ϵ}_7{s}^{q-1}\right)e^{-s\theta }}\hfill \\ {\displaystyle +\left(2q{ϵ}_9{s}^{2q-1}+q{ϵ}_{10}{s}^{q-1}\right)e^{-2s\theta },}\hfill \\ {\displaystyle Z_2\left(s\right)=s{e}^{-s\theta }({ϵ}_5{s}^{3q}+{ϵ}_6{s}^{2q}+{ϵ}_7{s}^q+{ϵ}_8)}\hfill \\ {\displaystyle +2s{e}^{-2s\theta }\left({ϵ}_9{s}^{2q}+{ϵ}_{10}{s}^q+{ϵ}_{11}\right)}\hfill \\ {\displaystyle +3s{ϵ}_{12}{e}^{-3s\theta }.}\hfill \end{array}\right.$$

(85)

Then,

$$\mathrm{Re}\left\{\frac{ds}{d\theta }\right\}{|}_{\theta ={\theta }_{\diamond },\vartheta ={\vartheta }_0}=\mathrm{Re}\left\{\frac{Z_1\left(s\right)}{Z_2\left(s\right)}\right\}{|}_{\theta ={\theta }_{\diamond },\vartheta ={\vartheta }_0}=\frac{Z_{11}{Z}_{21}+Z_{12}{Z}_{22}}{Z_{21}^2+Z_{22}^2}.$$

(86)

By $\left(A_6\right)$, we have the following:

$$\mathrm{Re}\left\{{\left[\frac{ds}{d\theta }\right]}^{-1}\right\}{|}_{\theta ={\theta }_{\diamond },\vartheta ={\vartheta }_0}>0.$$

(87)

□

Next, we provide the following assumption:

($A_7)$ The following inequalities hold:

$$\left\{\begin{array}{c}{\displaystyle {\Phi }_1={ϵ}_1>0,}\hfill \\ {\displaystyle {\Phi }_2=det\left[\begin{array}{cc}{ϵ}_1& 1\\ {ϵ}_3+{ϵ}_6+{ϵ}_{10}& {ϵ}_2+{ϵ}_5+{ϵ}_9\end{array}\right]>0,}\hfill \\ {\displaystyle {\Phi }_3=det\left[\begin{array}{ccc}{ϵ}_1& 1& 0\\ {ϵ}_3+{ϵ}_6+{ϵ}_{10}& {ϵ}_2+{ϵ}_5+{ϵ}_9& {ϵ}_1\\ 0& {ϵ}_4+{ϵ}_8+{ϵ}_{11}+{ϵ}_{12}& {ϵ}_3+{ϵ}_6+{ϵ}_{10}\end{array}\right]>0,}\hfill \\ {\displaystyle {\Phi }_4=({ϵ}_3+{ϵ}_6+{ϵ}_{10}){\Phi }_3>0.}\hfill \end{array}\right.$$

(88)

Lemma 10. 

If $\theta =0$ and $\left(A_7\right)$ is satisfied, then system (66) remains locally asymptotically stable.

Proof. 

When $\theta =0$, then (66) becomes

$${\lambda }^4+{ϵ}_1{\lambda }^3+({ϵ}_2+{ϵ}_5+{ϵ}_9){\lambda }^2+({ϵ}_3+{ϵ}_6+{ϵ}_{10})\lambda +{ϵ}_4+{ϵ}_8+{ϵ}_{11}+{ϵ}_{12}=0.$$

(89)

By virtue of $\left(A_7\right)$, one knows that every root ${\lambda }_i$ of (5.28) obeys $|\arg \left({\lambda }_j\right)|>\frac{q\pi }{2}(j=1,2,3,4).$ Thus, Lemma 10 is correct. □

Relying on the discussion above, the following outcome is lightly acquired.

Theorem 3. 

If $\left(A_1\right),\left(A_6\right),\left(A_7\right)$ are fulfilled, then the positive equilibrium point $W_0(w_{1\ast },w_{2\ast }$, $w_{3\ast },w_{4\ast })$ of system (62) remains locally asymptotically stable if θ falls into the range of $[0,{\theta }_{\diamond })$ and system (62) will generate a Hopf bifurcation around the positive equilibrium point $W_0(w_{1\ast },w_{2\ast },w_{3\ast },w_{4\ast })$ when $\theta ={\theta }_{\diamond }.$

Remark 1. 

The choice of the PD controller and the hybrid controller depends on the fractional-order dynamical models. By analyzing the characteristic equation of the controlled model and carrying out computer simulations, then we decide on the choice of controllers. Compared with the previous controllers, the advantages of the PD controller and the hybrid controller lie in their more flexible parameter adjustment. Of course, their disadvantages lie in the computational complexity of the distribution of characteristic roots for the corresponding characteristic equations.

Remark 2. 

Jia et al. [45] investigated the stability and bifurcation of integer-order system (2). In this paper, we explore the stability and bifurcation of fractional-order system (3). Different from model (2), a fractional-order parameter is added in model (3). The characteristic equation of model (3) is more complex than that of model (2). In addition, two different controllers are designed to control the stability and Hopf bifurcation. The analysis of model (3) enriches the bifurcation theory of the fractional-order dynamical system, to some degree.

6. Simulation Outcomes

In this section, we use the standard Euler methods [47] to carry out numerical simulations.

Example 1. 

Consider the following FitzHugh–Nagumo neural model:

$$\left\{\begin{array}{c}{\displaystyle \frac{d{w}_1^{0.91}\left(t\right)}{d{t}^{0.91}}=0.33\left(w_2+w_1-\frac{1}{3}{w}_1^3\right)+1.5(w_1-w_3(t-\theta )),}\hfill \\ {\displaystyle \frac{d{w}_2^{0.91}\left(t\right)}{d{t}^{0.91}}=-\frac{1}{0.33}(w_1-1.88+0.18w_2),}\hfill \\ {\displaystyle \frac{d{w}_3^{0.91}\left(t\right)}{d{t}^{0.91}}=0.33\left(w_4+w_3-\frac{1}{3}{w}_3^3\right)+1.5(w_3-w_1(t-\theta )),}\hfill \\ {\displaystyle \frac{d{w}_4^{0.91}\left(t\right)}{d{t}^{0.91}}=-\frac{1}{0.33}(w_3-1.88+0.18w_4).}\hfill \end{array}\right.$$

(90)

One can lightly see that system (90) admits a unique positive equilibrium point $W_0(w_{1\ast },w_{2\ast }$, $w_{3\ast },w_{4\ast })$=$W_0(1.8382,,1.8382,0.2322)$. Due to the Matlab software, one can lightly obtain ${\varphi }_0=0.7499$ and ${\theta }_0=0.62$. We can check that $\left(A_1\right)$–$\left(A_3\right)$ of Theorem 1 are satisfied. Thus, when $\theta \in [0,0.62)$, then the positive equilibrium point $W_0(1.8382,0.2322,1.8382,0.2322)$ of system (90) remains locally asymptotically stable. Select $\theta =0.6<{\theta }_0=0.62$. The simulation outcomes are presented in Figure 1. One can see from Figure 1 that if θ is less than the delay point ${\theta }_0=0.62$, then $w_1\left(t\right)\to 1.8382,w_2\left(t\right)\to 0.2322,w_3\left(t\right)\to 1.8382,w_4\left(t\right)\to 0.2322$, respectively, as $t\to \infty$. When θ crosses the delay point ${\theta }_0=0.62$, then system (90) loses its stability and generates a cluster of periodic solutions (namely, Hopf bifurcations) at around $W_0(1.8382,0.2322,1.8382,0.2322)$. Select $\theta =0.635$. The simulation outcomes are presented in Figure 2. One can see from Figure 2 that ϑ is greater than the delay point ${\theta }_0=0.62$, then $(w_1\left(t\right),w_2\left(t\right),w_3\left(t\right),w_4\left(t\right))$ is a periodic oscillatory state around $W_0(1.8382,0.2322,1.8382,0.2322)$. To further illustrate the Hopf bifurcation situation, we present the corresponding bifurcation plots (see Figure 3, Figure 4, Figure 5 and Figure 6). Depending on Figure 3, Figure 4, Figure 5 and Figure 6, we can understand the relationship between θ-$w_1$, θ-$w_2$, θ-$w_3$, θ-$w_4$, respectively. Obviously, the bifurcation value of system (90) is $0.62$. Furthermore, the quantitative relation between ${\varphi }_0$ and ${\theta }_0$ is provided via

Table 1. The quantitative relation between ${\varphi }_0$ and ${\theta }_0$ of system (90).

${\mathit{\varphi }}_0$	${\mathit{\theta }}_0$
$0.1805$	$0.1231$
$0.3488$	$0.2519$
$0.4602$	$0.3452$
$0.6169$	$0.4879$
$0.6324$	$0.5027$
$0.7499$	$0.6200$
$0.8753$	$0.7543$
$0.9402$	$0.8276$
$1.0220$	$0.9238$

.

Remark 3. 

In order to achieve local asymptotic stability, we explore the characteristic equation of system (90) according to Section 3. In view of Theorem 1, we know the local asymptotic stability condition of system (90).

Remark 4. 

In the simulation, since $0<q<1$, we can choose the arbitrary q, which falls into the range of $(0,1)$. In this section, we choose $q=0.91$. Of course, we can choose $(0.1,0.5)$. For this case, we can also carry out the numerical simulations for this value for $q\in (0.1,0.5)$. Here, we omit it. But this does not affect the verification of Theorem 1.

Example 2. 

Consider the following fractional-order controlled FitzHugh–Nagumo neural model:

$$\left\{\begin{array}{c}{\displaystyle \frac{d^{0.91}{w}_1\left(t\right)}{d{t}^{0.91}}=0.33\left(w_2+w_1-\frac{1}{3}{w}_1^3\right)+1.5(w_1-w_3(t-\theta )),}\hfill \\ {\displaystyle \frac{d^{0.91}{w}_2\left(t\right)}{d{t}^{0.91}}=-\frac{1}{0.33}(w_1-1.88+0.18w_2)+{\rho }_p(w_2\left(t\right)-w_{2\ast })+{\rho }_d\frac{d^{0.91}(w_2\left(t\right)-w_{2\ast })}{d{t}^{0.91}},}\hfill \\ {\displaystyle \frac{d^{0.91}{w}_3\left(t\right)}{d{t}^{0.91}}=0.33\left(w_4+w_3-\frac{1}{3}{w}_3^3\right)+1.5(w_3-w_1(t-\theta )),}\hfill \\ {\displaystyle \frac{d^{0.91}{w}_4\left(t\right)}{d{t}^{0.91}}=-\frac{1}{0.33}(w_3-1.88+0.18w_4)+{\varrho }_p(w_4\left(t\right)-w_{4\ast })+{\varrho }_d\frac{d^{0.91}(w_4\left(t\right)-w_{4\ast })}{d{t}^{0.91}}.}\hfill \end{array}\right.$$

(91)

One can lightly see that system (91) admits a unique positive equilibrium point $W_0(w_{1\ast },w_{2\ast }$, $w_{3\ast },w_{4\ast })=W_0(1.8382,1.8382,0.2322)$. Set ${\rho }_p=0.2,{\rho }_d=0.8$. Due to the Matlab software, one can lightly obtain ${\phi }_0=0.6877$ and ${\theta }_{\ast }=0.76$. We can check that $\left(A_1\right),\left(A_4\right),\left(A_5\right)$ of Theorem 1 are satisfied. Thus, when $\theta \in [0,0.76)$, then the positive equilibrium point $W_0(1.8382,0.2322,1.8382$, $0.2322)$ of system (91) remains locally asymptotically stable. Select $\theta =0.75<{\theta }_{\ast }=0.76$. The simulation outcomes are presented in Figure 7. One can see from Figure 7 that if θ is less than the delay point, ${\theta }_{\ast }=0.76$, then $w_1\left(t\right)\to 1.8382,w_2\left(t\right)\to 0.2322,w_3\left(t\right)\to 1.8382,w_4\left(t\right)\to 0.2322$, respectively, as $t\to \infty$. When θ crosses the delay point, ${\theta }_{\ast }=0.76$, then system (91) loses its stability and generates a cluster of periodic solutions (namely, Hopf bifurcations) at around $W_0(1.8382,0.2322,1.8382,0.2322)$. Select $\theta =0.8$. The simulation outcomes are presented in Figure 8. One can see in Figure 8 that if ϑ is greater than the delay point, ${\theta }_{\ast }=0.76$, then $(w_1\left(t\right),w_2\left(t\right),w_3\left(t\right),w_4\left(t\right))$ remains a periodic oscillatory state around $W_0(1.8382,0.2322$., $1.8382,0.2322)$. To further illustrate the Hopf bifurcation situation, we present the corresponding bifurcation plots (see Figure 9, Figure 10, Figure 11 and Figure 12). Depending on Figure 9, Figure 10, Figure 11 and Figure 12, we can understand the relationship among θ-$w_1$, θ-$w_2$, θ-$w_3$, and θ-$w_4$, respectively. Obviously, the bifurcation value of system (91) is $0.76$. Furthermore, the quantitative relation between ${\phi }_0$ and ${\theta }_{\ast }$ is presented in

Table 2. The quantitative relation between ${\phi }_0$ and ${\theta }_{\ast }$ of system (91).

${\mathit{\phi }}_0$	${\mathit{\theta }}_{\ast }$
$0.2103$	$0.1789$
$0.2790$	$0.2467$
$0.3383$	$0.3092$
$0.4566$	$0.4456$
$0.0993$	$0.5672$
$0.6089$	$0.6453$
$0.6877$	$0.7600$
$0.7268$	$0.8199$
$0.8112$	$0.9562$

.

Example 3. 

Consider the following fractional-order controlled FitzHugh–Nagumo neural model:

$$\left\{\begin{array}{c}{\displaystyle \frac{d^{0.91}{w}_1\left(t\right)}{d{t}^{0.91}}=0.33\left(w_2+w_1-\frac{1}{3}{w}_1^3\right)+1.5(w_1-w_3(t-\theta )),}\hfill \\ {\displaystyle \frac{d^{0.91}{w}_2\left(t\right)}{d{t}^{0.91}}={\kappa }_1\left[-\frac{1}{0.33}(w_1-1.88+0.18w_2)\right]+{\kappa }_2[w_2(t-\theta )-w_2\left(t\right)],}\hfill \\ {\displaystyle \frac{d^{0.91}{w}_3\left(t\right)}{d{t}^{0.91}}=0.33\left(w_4+w_3-\frac{1}{3}{w}_3^3\right)+1.5(w_3-w_1(t-\theta )),}\hfill \\ {\displaystyle \frac{d^{0.91}{w}_4\left(t\right)}{d{t}^{0.91}}={\kappa }_1\left[-\frac{1}{0.33}(w_3-1.88+0.18w_4)\right]+{\kappa }_2[w_4(t-\theta )-w_4\left(t\right)].}\hfill \end{array}\right.$$

(92)

One can lightly see that system (92) admits a unique positive equilibrium point $W_0(w_{1\ast },w_{2\ast },w_{3\ast },w_{4\ast })=W_0(1.8382,,1.8382,0.2322)$. Set ${\kappa }_1=4.2,{\kappa }_2=0.2$. Due to the Matlab software, one can lightly obtain ${\vartheta }_0=0.2976$ and ${\theta }_{\diamond }=0.46$. We can check that $\left(A_1\right),\left(A_6\right),\left(A_7\right)$ of Theorem 1 are satisfied. Thus, when $\theta \in [0,0.46)$, then the positive equilibrium point $W_0(1.8382,0.2322,1.8382,0.2322)$ of system (92) remains locally asymptotically stable. Select $\theta =0.45<{\theta }_{\diamond }=0.46$. The simulation outcomes are presented in Figure 13. One can find from Figure 13 that if θ is less than the delay point, ${\theta }_{\diamond }=0.46$, then $w_1\left(t\right)\to 1.8382,w_2\left(t\right)\to 0.2322,w_3\left(t\right)\to 1.8382,w_4\left(t\right)\to 0.2322$, respectively, as $t\to \infty$. When θ crosses the delay point, ${\theta }_{\diamond }=0.46$, then system (92) loses its stability and generates a cluster of periodic solutions (namely, Hopf bifurcations) at around $W_0(1.8382,0.2322,1.8382,0.2322)$. Select $\theta =0.473$. The simulation outcomes are presented in Figure 14. One can find from Figure 14 that ϑ is greater than the delay point, ${\theta }_{\diamond }=0.46$, then $(w_1\left(t\right),w_2\left(t\right),w_3\left(t\right),w_4\left(t\right))$ remains in a periodic oscillatory state at around $W_0(1.8382,0.2322,1.8382,0.2322)$. To further illustrate the Hopf bifurcation situation, we present the corresponding bifurcation plots (see Figure 15, Figure 16, Figure 17 and Figure 18). Depending on Figure 15, Figure 16, Figure 17 and Figure 18, we can understand the relationship between θ-$w_1$, θ-$w_2$, θ-$w_3$, and θ-$w_4$, respectively. Obviously, the bifurcation value of system (92) is $0.46$. Furthermore, the quantitative relation between ${\vartheta }_0$ and ${\theta }_{\diamond }$ is presented in

Table 3. The quantitative relation between ${\vartheta }_0$ and ${\theta }_{\diamond }$ of system (92).

${\mathit{\vartheta }}_0$	${\mathit{\theta }}_{\diamond }$
$0.1071$	$0.1542$
$0.1879$	$0.2789$
$0.2358$	$0.3563$
$0.2976$	$0.4600$
$0.3525$	$0.5562$
$0.3934$	$0.6303$
$0.4482$	$0.7328$
$0.5092$	$0.8513$
$0.5496$	$0.9327$

.

Remark 5. 

According to the simulation outcomes of Examples 1–3, we know that the bifurcation values of system (90), system (91), and system (92) are ${\theta }_0=0.62,{\theta }_{\ast }=0.76,and{\theta }_{\diamond }=0.46,$ respectively, which implies that (i) we can enlarge the stability domain and postpone the timing of bifurcation generation of the Hopf bifurcation of system (90) via the $P{D}^p$ controller and (ii) we can narrow the stability domain and advance the timing of bifurcation generation of the Hopf bifurcation of system (90) via the hybrid controller.

7. Conclusions

Depending on earlier studies, we formulate a novel fractional-order coupled FitzHugh–Nagumo neural model with delay. By analyzing the characteristic equation of the formulated fractional-order coupled FitzHugh–Nagumo neural model, we discuss the stability and Hopf bifurcation issue. A delay-independent criterion on stability and Hopf bifurcation for the considered coupled FitzHugh–Nagumo neural model is obtained. Due to the two effective controllers ($P{D}^p$ controller and hybrid controller), the stability domain and the Hopf bifurcation onset time of the formulated fractional-order coupled FitzHugh–Nagumo neural model are successfully dominated. The influence of delay on the stability domain and the Hopf bifurcation onset time is displayed. The acquired theoretical outcomes of this study possess great theoretical value in operating and designing networks. By adjusting the control parameter and time delay, we can effectively control the stability domain and the time of bifurcation onset in neural networks. The control techniques can also be applied to control the stability and bifurcation behavior of other neural networks. In addition, we can control the stability domain and the bifurcation onset time in model (3) via other controllers. This will be explored in future work.