Stability and Bifurcation in a Hopfield Neuron Model with Delays

Suqi Ma¹, S. J. Hogan²

¹Department of Mathematics, China Agricultural University, Beijing, China.
²Department of Mathematics, Bristol University, Bristol, UK.

DOI: 10.4236/ijmnta.2025.141001                

Abstract

The Hopfield system is an artificial neuron model that can be applied to neuron memory and information processing. Like synaptic connections of neurons, the inhibition or excitable feedback often incorporate delay effects, which are either discrete or distributed time delays. With delays varying, Hopf bifurcation of distributed time delay is investigated and stability regime is partitioned by Hopf curves on the parameter plane. Lyapunov-Schmidt reduction skills combined with center manifold theory are applied to discuss the stability of bifurcating periodical solutions arising from Hopf points. In addition, DDE-Biftool software significantly provides the numerical computation of stability analysis of periodical solutions appearing in discrete time delay Hopfield system. The period-doubling bifurcation of periodical solutions, which form P2 circles and P4 circles, respectively, by continuous periodical solutions with varying free parameters, is discussed.

Keywords

Hopfield Network Model, Distributed Time Delay, Center Manifold, Period Doubling Bifurcation

1. Introduction

The Hopfield network model is an artificial neuron biological model with fundamental inspiration to deepen people’s understanding of associative memory of the human brain. The Hopfield network utilize connections to store and retrieve different patterns, mimicking natural efficiency in information processing, like synaptic connections of neurons [1]-[4]. The Hopfield network model explores the intricate connections among neurons, which often expresses the delayed excitable or inhibition feedback due to the neuron memory with discrete or distributed time delays [5] [6].

As is well known, people investigate the Hopfield network model with its oscillation dynamical behavior undertaken Hopf bifurcation occurring. The authors in paper [1] [7]-[9] have reported Hopf bifurcation with Hopfield network model underlying its symmetrical character. Some authors also explore second Hopf bifurcation of Hopfield network model and both Chaos and hyperchaos atractors are reported [7] [10]-[12]. However, with the known results of spatial symmetry of neuron models, which manifests significant importance in model efficiency, periodical transition phenomena have gradually become an interesting topic in the relevant field [13] [14]. Inspired by the cells prosperous dynamical phenomena observed in numerical simulation results, we develop the period doubling bifurcation of periodical solutions, which has bifurcating branches with doubly periods and usually leads the routes to chaos or quasi-periodical attractors. We know scheme in this paper is mainly dependent on an artificial handbook named DDE-Biftool software [15] [16]. Our referred work devastating and attracting the experiences and skills in conducting the continuous periodical solutions with bifurcating branches of period doubling bifurcation too [15]-[17].

The often neuron Hopfield model with distributed delay is put forward as the following

$\begin{array}{l}{u^{\prime }}_1\left(t\right)=-b_1{u}_1\left(t\right)+c_1\tanh \left(u_1\left(t-{\tau }_1\right)\right)+a_{12}{c}_2{\displaystyle {\int }_0^{{\tau }_2}\tanh }\left(u_2\left(t-{\tau }_2+s\right)\right)\text{d}s,\\ {u^{\prime }}_2\left(t\right)=-b_2{u}_2\left(t\right)+a_{21}{c}_1{\displaystyle {\int }_0^{{\tau }_2}\tanh }\left(u_1\left(t-{\tau }_2+s\right)\right)\text{d}s+c_2\tanh \left(u_2\left(t-{\tau }_1\right)\right)\end{array}$ (1.1)

Wherein, three distributed time delays period with average delay ${\tau }_i\left(i=1,2,3\right)$ are referred. The coefficients $a_{12},a_{21}$ represent the two mode of information transmitting between neurons, and $c_1,c_2$ describes which neuron is invoked to link up.

The simple version of discrete time delays model of the above Hopfield network is also described as

$\begin{array}{l}{u^{\prime }}_1\left(t\right)=-b_1{u}_1\left(t\right)+c_1\tanh \left(u_1\left(t-{\tau }_1\right)\right)+a_{12}{c}_2\tanh \left(u_2\left(t-{\tau }_2\right)\right),\\ {u^{\prime }}_2\left(t\right)=-b_2{u}_2\left(t\right)+a_{21}{c}_1\tanh \left(u_1\left(t-{\tau }_2\right)\right)+c_2\tanh \left(u_2\left(t-{\tau }_1\right)\right)\end{array}$ (1.2)

The interesting oscillating phenomena of system (1.2) are explored as shown in Figure 1(a) and Figure 1(b). Both the chaos and quasi-periodical attractors are found with free parameters $c_2$ and time delays ${\tau }_i\left(i=1,2\right)$ chosen. With fixed parameter $a_{12}=1$ , $a_{21}=1.2$ , system (1.2) manifests periodical oscillation with nearly spatial symmetry. With the progress to be remarkable, we expand doubling period bifurcation of system with the helpful job of Floquet multiplier computed by DDE-Biftool. As for system (1.1), we expand the stability analysis of trivial solution and further reduction system combined with the center manifold theory [11] [18] [19]. Hopf bifurcation occurs if a pair of imaginary roots cross over the imaginary axis and simultaneously, the transversal condition is satisfied. We apply Lyapunov-Schmidt reduction method to derive the normal formal by the trunction system of (1.2) by projecting solution operator onto the center subspace. The stability of bifurcating periodical solution is determined by the normal form and the bifurcation direction is computed. The example of continuing the bifurcating periodical solution by DDE-Biftool is done by continuously varying two-time delays.

Figure 1. The nearly symmetry attractors of system (1.2). (a) with $c_1=-2$ , $c_2=0.9$ , ${\tau }_1=1.2$ , ${\tau }_2=19$ ; (b) $c_1=-2$ , $c_2=0.6$ , ${\tau }_1=1.2$ , ${\tau }_2=26$ .

Since the distributed time delay Equation (1.1) has a zero characteristic root undertaken Hopf bifurcation happens, the bifurcation of periodical oscillation phenomena is difficult. We transform system (1.1) into DDEs as follows

$\begin{array}{l}{u^{\prime }}_1\left(t\right)=-b_1{u}_1\left(t\right)+c_1\tanh \left(u_1\left(t-{\tau }_1\right)\right)+a_{12}{c}_2{u}_3,\\ {u^{\prime }}_2\left(t\right)=-b_2{u}_2\left(t\right)+a_{21}{c}_1{u}_4+c_2\tanh \left(u_2\left(t-{\tau }_1\right)\right),\\ {u^{\prime }}_3\left(t\right)=\tanh \left(u_2\left(t\right)\right)-\tanh \left(u_2\left(t-{\tau }_2\right)\right),\\ {u^{\prime }}_4\left(t\right)=\tanh \left(u_1\left(t\right)\right)-\tanh \left(u_1\left(t-{\tau }_2\right)\right)\end{array}$ (1.3)

The bifurcating periodical solution is continuously continued by varying time delays. For system (1.2), we replace distributed terms in the Equation by its discrete time delays parts, which is more easily to result in limit cycle bifurcation by DDE-Biftool software. The period doubling bifurcation of periodical solution is found by Floquent multiplier attains at −1. The doubly period solutions of P2 and P4 periodical oscillation are also continued by varying free parameters. The P2 and P4 islands of periodical solutions are simulated.

The whole paper is organized as listed. In Section 2, the stability analysis of system (1.1) is done and Hopf bifurcation arises as the system loss its stability. In Section 3, the normal form is computed by applying the Lyapunov-Schmidt reduction skills combined with the center manifold theory. In Section 4, the period-doubling bifurcation branches the P2 and P4 periodical solutions in system (1.2) is simulated by software. The conclusion is given finally.

2. Stability Analysis

The stability property of the distributed time delay system (1.1) is investigated. System (1.1) loss stability as the rightmost characteristic root with positive real part. The linear version of system (1.1) is written as the following

$\begin{array}{l}{u^{\prime }}_1\left(t\right)=-b_1{u}_1\left(t\right)+c_1{u}_1\left(t-{\tau }_1\right)+a_{12}{c}_2{\displaystyle {\int }_0^{{\tau }_2}{u}_2}\left(t-{\tau }_2+s\right)\text{d}s,\\ {u^{\prime }}_2\left(t\right)=-b_2{u}_2\left(t\right)+a_{21}{c}_1{\displaystyle {\int }_0^{{\tau }_2}{u}_1}\left(t-{\tau }_2+s\right)\text{d}s+c_2{u}_2\left(t-{\tau }_1\right)\end{array}$ (2.1)

The characteristic Equation of Equation (2.1) can be written as

$\text{Δ}\left(\lambda \right)=\left|\begin{array}{cc}-b_1-\lambda -\frac{c_1}{{\tau }_1}\frac{1}{\lambda }\left({\text{e}}^{-\lambda {\tau }_1}-1\right)& -a_{12}\frac{c_2}{{\tau }_2}\frac{1}{\lambda }\left({\text{e}}^{-\lambda {\tau }_2}-1\right)\\ -a_{21}\frac{c_1}{{\tau }_1}\frac{1}{\lambda }\left({\text{e}}^{-\lambda {\tau }_1}-1\right)& -b_2-\lambda -\frac{c_2}{{\tau }_2}\frac{1}{\lambda }\left({\text{e}}^{-\lambda {\tau }_2}-1\right)\end{array}\right|$ (2.2)

For simplicity, Equation (2.2) is rewritten as

$\begin{array}{l}\Delta \left(\lambda ,{\tau }_1,{\tau }_2\right)=H\left(\lambda \right){\text{e}}^{-\lambda {\tau }_1}+P\left(\lambda \right)+Q\left(\lambda \right){\text{e}}^{-\lambda {\tau }_2}+W\left(\lambda \right){\text{e}}^{-2\lambda {\tau }_1}+G\left(\lambda \right){\text{e}}^{-2\lambda {\tau }_2}=0\hfill \end{array}$ (2.3)

with

$\begin{array}{l}P\left(\lambda \right)=-b_1^2{\lambda }^2-\left(b_1+b_2\right){\lambda }^3-{\lambda }^4+a_{12}{c}_2{a}_{21}{c}_1,\\ Q\left(\lambda \right)=-2a_{12}{c}_2{a}_{21}{c}_1,\\ H\left(\lambda \right)=b_1{\lambda }^2{c}_2+{\lambda }^3{c}_2+b_2{\lambda }^2{c}_1+{\lambda }^3{c}_1,\\ W\left(\lambda \right)=-{\lambda }^2{c}_1{c}_2,G\left(\lambda \right)=a_{12}{c}_2{a}_{21}{c}_1\end{array}$

By setting $\lambda =i\omega ,\left(\omega >0\right)$ , we calculate the imaginary roots and seek for the critical value $\left({\tau }_1,{\tau }_2\right)$ of Hopf bifurcation. Setting the real part from the imaginary part of the coefficients as mentioned above,

$\begin{array}{l}{P}_1=\Re \left(P\left(i\omega \right)\right),P_2=\Im \left(P\left(i\omega \right)\right),Q_1=\Re \left(Q\left(i\omega \right)\right),Q_2=\Im \left(Q\left(i\omega \right)\right),\\ H_1=\Re \left(H\left(i\omega \right)\right),H_2=\Im \left(H\left(i\omega \right)\right),W_1=\Re \left(W\left(i\omega \right)\right),W_2=\Im \left(W\left(i\omega \right)\right),\\ G_1=\Re \left(G\left(i\omega \right)\right),G_2=\Im \left(G\left(i\omega \right)\right).\end{array}$

We also set

$\varphi =\mathrm{mod}\left(\omega {\tau }_1,2\pi \right)+2n\pi ,\psi =\mathrm{mod}\left(\omega {\tau }_2,2\pi \right)+2m\pi$ (2.4)

for $m,n=0,1,2,\cdots$ . Furthermore, one substitute $\varphi ,\psi$ into Equation (2.3), by the related triangle equality to solve

$\begin{array}{l}\cos \left(\varphi \right)=-\frac{g1\cos \left(2\psi \right)+g_2\sin \left(2\psi \right)+g_3\cos \left(\psi \right)+g_4\sin \left(\psi \right)+g_5}{HH},\\ \sin \left(\varphi \right)=-\frac{g_1\sin \left(2\psi \right)-g_2\cos \left(2\psi \right)-g_4\cos \left(\psi \right)+g_3\sin \left(\psi \right)+g_6}{HH},\\ HH=d_1\sin \left(\psi \right)\cos \left(2\psi \right)+d_2\sin \left(\psi \right)\sin \left(2\psi \right)+d_3\cos \left(\psi \right)\cos \left(2\psi \right)\\ +d_4\cos \left(\psi \right)\sin \left(2\psi \right)+d_5\cos \left(2\psi \right)+d_7\cos \left(\psi \right)+d_6\sin \left(2\psi \right)\\ +d_8\sin \left(\psi \right)+d_9.\end{array}$

with

$\begin{array}{l}{g}_1=G_1{H}_1+G_2{H}_2,g_2=-G_1{H}_2+G_2{H}_1,g_3=H_1{Q}_1+H_2{Q}_2,\\ g_4=H_1{Q}_2-H_2{Q}_1,g_5=H_1{P}_1-H_1{W}_1+H_2{P}_2-H_2{W}_2,\\ g_6=-H_1{P}_2-H_1{W}_2+H_2{P}_1+H_2{W}_1,\\ d_1=-2G_2{Q}_1+2G_1{Q}_2,d_2=2G_2{Q}_2+2G_1{Q}_1,d_3=2G_1{Q}_1+2G_2{Q}_2,\\ d_4=-2G_1{Q}_2+2G_2{Q}_1,d_5=2G_1{P}_1+2G_2{P}_2,d_6=2G_2{P}_1-2G_1{P}_2,\\ d_7=2P_2{Q}_2+2P_1{Q}_1,d_8=2P_1{Q}_2-2P_2{Q}_1,\\ d_9=-W_2^2-W_1^2+P_2^2+P_1^2+Q_1^2+Q_2^2+G_1^2+G_2^2.\end{array}$

Set the functions

$\begin{array}{l}F\left(\psi \right)=\cos {\left(\varphi \right)}^2+\sin {\left(\varphi \right)}^2-1\equiv 0,\\ G\left(\varphi ,\psi \right)=\tan \left(\varphi \right)-\frac{g_1\sin \left(2\psi \right)-g_2\cos \left(2\psi \right)-g_4\cos \left(\psi \right)+g_3\sin \left(\psi \right)+g_6}{g1\cos \left(2\psi \right)+g_2\sin \left(2\psi \right)+g_3\cos \left(\psi \right)+g_4\sin \left(\psi \right)+g_5}\end{array}$ (2.5)

By Equation (2.3), we also get

$\begin{array}{c}D\left(\varphi ,\psi \right)=(H_1\cos \left(\varphi \right)+H_2\sin \left(\varphi \right)+P_1+Q_1\cos \left(\psi \right)+Q_2\sin \left(\psi \right)+G_1\cos \left(2\psi \right)\\ +{G_2\sin \left(2\psi \right))}^2+(-H_1\sin \left(\varphi \right)+P_2+Q_2\cos \left(\psi \right)-Q_1\sin \left(\psi \right)\\ +H_2\cos \left(\varphi \right)+G_2\cos \left(2\psi \right)-{G_1\sin \left(2\psi \right))}^2-W_1^2-W_2^2\end{array}$ (2.6)

By Equations (2.5) and Equation (2.6), we solve $\left(\varphi ,\psi ,\omega \right)$ which satisfy the characteristic Equation, hence after Hopf bifurcation value $\left({\tau }_1,{\tau }_2\right)$ are given as

$\{\begin{array}{ll}{\tau }_1=\frac{\varphi +2n\pi }{\omega },\hfill & \text{for}n=0,1,2,\cdots \hfill \\ {\tau }_2=\frac{\psi +2m\pi }{\omega },\hfill & \text{for}m=0,1,2,\cdots \hfill \end{array}$ (2.7)

For example, fixed parameter $a_{12}=1$ , $a_{21}=1.2$ , $c_1=-2$ , $c_2=0.6023012058$ , we derive $\left(\varphi ,\psi ,\omega \right)$ from Equations (2.5) and Equation (2.6) which is

$\begin{array}{l}\omega =1.7197,\cos \left(\varphi \right)=0.2172524861,\sin \left(\varphi \right)=0.9761154426,\\ \cos \left(\psi \right)=-0.3828472279,\sin \left(\psi \right)=0.9238116692\end{array}$

and the threshold values for Hopf bifurcation is listed as

$\begin{array}{l}{\tau }_{1,1}=0.7850148049,{\tau }_{2,1}=1.141869341,\\ {\tau }_{1,2}=4.438666784,{\tau }_{2,2}=4.795521320\end{array}$

We draw a picture of Hopf bifurcation lines on $\left({\tau }_1,{\tau }_2\right)$ -plane, as shown in Figure 2. The stability property of the trivial solution is plotted by DDE-Biftool, as shown in Figure 2(a). Hopf curves are also drawn on $\left({\tau }_1,{\tau }_2\right)$ -plane, wherein the blue Hopf lines separate the stable regimes from the unstable regimes, whilst the red Hopf lines denote Hopf lines in the unstable regime, as shown in Figure 2(b). To determine the transversal condition for Hopf bifurcation, we compute the differential $\text{Δ}\left(\lambda ,{\tau }_1,{\tau }_2\right)$ with respect to its delay arguments to get

Figure 2. Hopf bifurcation on $\left({\tau }_1,{\tau }_2\right)$ parameter plane. (a) The imaginary roots at threshold value of Hopf bifurcation, which is sub-plotted with $\left({\tau }_1,{\tau }_2\right)=\left(0.7850148049,1.141869341\right)$ , $\left(0.7850148049,4.795521320\right)$ , $\left(4.438666784,1.141869341\right)$ , and $\left(4.438666784,4.795521320\right)$ ; (b) Hopf bifurcation curves on $\left({\tau }_1,{\tau }_2\right)$ parameter plane.

$\begin{array}{l}\frac{\text{d}\lambda }{\text{d}{\tau }_2}\left(\frac{\text{d}H}{\text{d}\lambda }{\text{e}}^{-\lambda {\tau }_1}+\frac{\text{d}P}{\text{d}\lambda }+\frac{\text{d}Q}{\text{d}\lambda }{\text{e}}^{-\lambda {\tau }_2}+\frac{\text{d}W}{\text{d}\lambda }{\text{e}}^{-2\lambda {\tau }_1}+\frac{\text{d}G}{\text{d}\lambda }{\text{e}}^{-2\lambda {\tau }_2}\right)\\ +\frac{\text{d}\lambda }{\text{d}{\tau }_2}\left(H{\text{e}}^{-\lambda {\tau }_1}\left(-{\tau }_1\right)+Q{\text{e}}^{-\lambda {\tau }_2}\left(-{\tau }_2\right)+W{\text{e}}^{-2\lambda {\tau }_1}\left(-2{\tau }_1\right)+G{\text{e}}^{-2\lambda {\tau }_2}\left(-2{\tau }_2\right)\right)\\ +Q{\text{e}}^{-\lambda {\tau }_2}\left(-\lambda \right)+G{\text{e}}^{-2\lambda {\tau }_2}\left(-2\lambda \right)=0\end{array}$ (2.8)

Therefore, we have

$\begin{array}{c}\frac{1}{\frac{\text{d}\lambda }{\text{d}{\tau }_2}}=\frac{\frac{\text{d}H}{\text{d}\lambda }{\text{e}}^{-\lambda {\tau }_1}+\frac{\text{d}P}{\text{d}\lambda }+\frac{\text{d}Q}{\text{d}\lambda }{\text{e}}^{-\lambda {\tau }_2}+\frac{\text{d}W}{\text{d}\lambda }{\text{e}}^{-2\lambda {\tau }_1}+\frac{\text{d}G}{\text{d}\lambda }{\text{e}}^{-2\lambda {\tau }_2}}{Q{\text{e}}^{-\lambda {\tau }_2}\left(-\lambda \right)+G{\text{e}}^{-2\lambda {\tau }_2}\left(-2\lambda \right)}\\ +\frac{H{\text{e}}^{-\lambda {\tau }_1}\left(-{\tau }_1\right)+Q{\text{e}}^{-\lambda {\tau }_2}\left(-{\tau }_2\right)+W{\text{e}}^{-2\lambda {\tau }_1}\left(-2{\tau }_1\right)+G{\text{e}}^{-2\lambda {\tau }_2}\left(-2{\tau }_2\right)}{Q{\text{e}}^{-\lambda {\tau }_2}\left(-\lambda \right)+G{\text{e}}^{-2\lambda {\tau }_2}\left(-2\lambda \right)}\end{array}$ (2.9)

Noticed if stretch along Hopf line, with given $\omega$ , one gets

$F\left(\psi \right)=0,G\left(\varphi ,\psi \right)=0$ (2.10)

that is, it easily computes that

$\{\begin{array}{ll}{\psi }_m={\psi }_0\left(\omega \right)+2m\pi \hfill & \Rightarrow {\tau }_{2,m}\left(\omega \right)=\frac{\psi }{\omega }\hfill \\ {\varphi }_n={\varphi }_0\left(\omega \right)+2n\pi \hfill & \Rightarrow {\tau }_{1,n}\left(\omega \right)=\frac{\varphi }{\omega }\hfill \end{array}$ (2.11)

wherein, from Equation (2.10), one calculates $\psi ,\varphi$ to get ${\psi }_0=\mathrm{mod}\left(\psi ,2\pi \right)$ , ${\varphi }_0=\mathrm{mod}\left(\varphi ,2\pi \right)$ . Hence, not mazed by Hopf line, we have

$\text{Δ}\left(i\omega ,{\tau }_{1,n},{\tau }_{2,m}\right)=0$ (2.12)

If differential Equation (2.12) with respect to $\omega$ , we also have

$\begin{array}{l}i\left(\frac{\text{d}H}{\text{d}\lambda }{\text{e}}^{-i\varphi }+\frac{\text{d}P}{\text{d}\lambda }+\frac{\text{d}Q}{\text{d}\lambda }{\text{e}}^{-i\psi }+\frac{\text{d}W}{\text{d}\lambda }{\text{e}}^{-i2\varphi }+\frac{\text{d}G}{\text{d}\lambda }{\text{e}}^{-i2\psi }\right)\\ +i\left(H{\text{e}}^{-i\varphi }\left(-{\tau }_1\right)+Q{\text{e}}^{-i\psi }\left(-{\tau }_2\right)+W{\text{e}}^{-i2\varphi }\left(-2{\tau }_1\right)+G{\text{e}}^{-i2\psi }\left(-2{\tau }_2\right)\right)\\ +{{\tau }^{\prime }}_{2,m}\left(\omega \right)\left(Q{\text{e}}^{-i\psi }\left(-i\omega \right)+G{\text{e}}^{-i2\psi }\left(-2i\omega \right)\right)\\ +{{\tau }^{\prime }}_{1,n}\left(\omega \right)\left(H{\text{e}}^{-i\varphi }\left(-i\omega \right)+W{\text{e}}^{-i2\varphi }\left(-2i\omega \right)\right)=0\end{array}$ (2.13)

Furthermore, one calculates the transversal condition from Equation (2.9) to get

$\begin{array}{l}\delta \left(\omega ,{\tau }_{1,n},{\tau }_{2,m}\right)={\frac{1}{\frac{\text{d}\lambda }{\text{d}{\tau }_2}}|}_{\lambda =i\omega }\\ =\frac{\frac{\text{d}H}{\text{d}\lambda }{\text{e}}^{-i\varphi }+\frac{\text{d}P}{\text{d}\lambda }+\frac{\text{d}Q}{\text{d}\lambda }{\text{e}}^{-i\psi }+\frac{\text{d}W}{\text{d}\lambda }{\text{e}}^{-i2\varphi }+\frac{\text{d}G}{\text{d}\lambda }{\text{e}}^{-i2\psi }}{Q{\text{e}}^{-i\psi }\left(-i\omega \right)+G{\text{e}}^{-i2\psi }\left(-2i\omega \right)}\\ +\frac{H{\text{e}}^{-\lambda {\tau }_1}\left(-{\tau }_1\right)+Q{\text{e}}^{-\lambda {\tau }_2}\left(-{\tau }_2\right)+W{\text{e}}^{-2\lambda {\tau }_1}\left(-2{\tau }_1\right)+G{\text{e}}^{-2\lambda {\tau }_2}\left(-2{\tau }_2\right)}{Q{\text{e}}^{-i\psi }\left(-i\omega \right)+G{\text{e}}^{-i2\psi }\left(-2i\omega \right)}\\ =-\frac{{{\tau }^{\prime }}_{2,m}\left(\omega \right)\left(Q{\text{e}}^{-i\psi }\left(-i\omega \right)+G{\text{e}}^{-i2\psi }\left(-2i\omega \right)\right)}{iQ{\text{e}}^{-i\psi }\left(-i\omega \right)+G{\text{e}}^{-i2\psi }\left(-2i\omega \right)}\\ -\frac{{{\tau }^{\prime }}_{1,n}\left(\omega \right)\left(H{\text{e}}^{-i\varphi }\left(-i\omega \right)+W{\text{e}}^{-i2\varphi }\left(-2i\omega \right)\right)}{iQ{\text{e}}^{-i\psi }\left(-i\omega \right)+G{\text{e}}^{-i2\psi }\left(-2i\omega \right)}\end{array}$

Therefore, one gets

$\begin{array}{c}\Re \left(\delta \left(\omega ,{\tau }_{1,n},{\tau }_{2,m}\right)\right)=-\Re \frac{{{\tau }^{\prime }}_{1,n}\left(\omega \right)\left(H{\text{e}}^{-i\varphi }-2W{\text{e}}^{-i2\varphi }\right)}{i\left(Q{\text{e}}^{-i\psi }-2G{\text{e}}^{-i2\psi }\right)}\\ =-{{\tau }^{\prime }}_{1,n}\left(\omega \right)\Re \frac{\left(H{\text{e}}^{-i\varphi }-2W{\text{e}}^{-i2\varphi }\right)\left(-i\left(\overline{Q}{\text{e}}^{i\psi }-2\overline{G}{\text{e}}^{i2\psi }\right)\right)}{{\left|i\left(Q{\text{e}}^{-i\psi }-2G{\text{e}}^{-i2\psi }\right)\right|}^2}\end{array}$ (2.14)

Based on the above discussion, one concludes that Hopf bifurcation occurs at $\left({\tau }_{1,n},{\tau }_{2,m}\right)$ if and only if the transversal condition is satisfied with$\Re \left(\delta \left(\omega ,{\tau }_{1,n},{\tau }_{2,m}\right)\right)\ne 0$ . The periodical solution bifurcates from Hopf point is simulated and stability analysis of periodical solution is carried out by normal form computation method in the next section.

3. Norm Form Analysis

As discussed in Section 2, Hopf bifurcation occurs at $\left({\tau }_1^{\text{*}},{\tau }_2^{\text{*}}\right)$ point while a pair of imaginary roots cross the imaginary axis with the transversal condition being satisfied. We apply Lyapunov-Schmidt reduction method combined with center manifold theory to compute the normal form near Hopf point. The perturbation method is explored to investigate the bifurcating direction of Hopf bifurcation and analyze the stability of the bifurcating periodical solution.

Set ${\tau }_1={\tau }_1^{\text{*}}+ϵ{\tau }_m,{\tau }_2={\tau }_2^{\text{*}}+ϵ{\tau }_n$ , and rewrite system (1.1) to its third trunction form as

$\begin{array}{c}{u^{\prime }}_1\left(t\right)=-b_1{u}_1\left(t\right)+c_1\left(u_1\left(t-{\tau }_1\right)-\frac{1}{3}{u}_1{\left(t-{\tau }_1\right)}^3\right)+a_{12}{c}_2{\displaystyle {\int }_0^{{\tau }_2}{u}_2}\left(t-{\tau }_2+s\right)\text{d}s\\ -\frac{1}{3}{a}_{12}{c}_2{\displaystyle {\int }_0^{{\tau }_2}{u}_2{\left(t-{\tau }_2+s\right)}^3\text{d}s}+o\left(\Vert u_1^3\Vert ,\Vert u_2^3\Vert ,\Vert u_1{\left(t-{\tau }_1\right)}^3\Vert ,\Vert u_2{\left(t-{\tau }_2\right)}^3\Vert \right)\\ {u^{\prime }}_2\left(t\right)=-b_2{u}_2\left(t\right)+a_{21}{\displaystyle {\int }_0^{{\tau }_2}{c}_1}\left(u_1\left(t-{\tau }_2+s\right)-\frac{1}{3}{u}_1{\left(t-{\tau }_2+s\right)}^3\right)\text{d}s\\ +c_2\left(u_2\left(t-{\tau }_1\right)-\frac{1}{3}{u}_2{\left(t-{\tau }_1\right)}^3\right)+o\left(\Vert u_1^3\Vert ,\Vert u_2^3\Vert ,\Vert u_1{\left(t-{\tau }_1\right)}^3\Vert ,\Vert u_2{\left(t-{\tau }_2\right)}^3\Vert \right)\end{array}$ (3.1)

Based on the fundamental theory of DDEs, Equations (3.1) is defined on its phase space $\varphi \in C\left(\left[-\tau ,0\right]\to R^2\right)$ , which is a Banach space with the super norm $\Vert \varphi \Vert ={\max }_{-\tau \le \theta \le 0}\left|\varphi \left(\theta \right)\right|$ with $\tau =\max \left({\tau }_1,{\tau }_2\right)$ . Furthermore, we write Equations (3.1) into the following nonlinear system

$\left(\begin{array}{c}{{\varphi }^{\prime }}_1\left(0\right)\\ {{\varphi }^{\prime }}_2\left(0\right)\end{array}\right)=ℒ\left(ϵ\right)\left(\begin{array}{c}{\varphi }_1\\ {\varphi }_2\end{array}\right)+F\left({\varphi }_1,{\varphi }_2\right)$ (3.2)

with the linearized version

$ℒ\left(0\right)\left(\begin{array}{c}{\varphi }_1\\ {\varphi }_2\end{array}\right)=\left(\begin{array}{c}-b_1{\varphi }_1\left(0\right)+c_1{\varphi }_1\left(-{\tau }_1^{\text{*}}\right)+a_{12}{c}_2{\displaystyle {\int }_0^{{\tau }_2^{\text{*}}}{\varphi }_2\left(-{\tau }_2^{\text{*}}+s\right)\text{d}s}\\ -b_2{\varphi }_2\left(0\right)+a_{21}{c}_1{\displaystyle {\int }_0^{{\tau }_2^{\text{*}}}{\varphi }_1\left(-{\tau }_2^{\text{*}}+s\right)\text{d}s}+c_2{\varphi }_2\left(-{\tau }_1^{\text{*}}\right)\end{array}\right)$ (3.3)

with its perturbation part

$\begin{array}{c}ℒ\left(ϵ\right)\left(\begin{array}{c}{\varphi }_1\\ {\varphi }_2\end{array}\right)=ℒ\left(0\right)\left(\begin{array}{c}{\varphi }_1\\ {\varphi }_2\end{array}\right)+c_1\left(\begin{array}{c}{\varphi }_1\left(-{\tau }_1\right)-{\varphi }_1\left(-{\tau }_1^{\text{*}}\right)\\ {\varphi }_2\left(-{\tau }_1\right)-{\varphi }_2\left(-{\tau }_1^{\text{*}}\right)\end{array}\right)\\ +\left(\begin{array}{c}{a}_{12}{c}_2\left({\displaystyle {\int }_0^{{\tau }_2}{\varphi }_2\left(-{\tau }_2+s\right)\text{d}s}-{\displaystyle {\int }_0^{{\tau }_2^{\text{*}}}{\varphi }_2\left(-{\tau }_2^{\text{*}}+s\right)\text{d}s}\right)\\ a_{21}{c}_1\left({\displaystyle {\int }_0^{{\tau }_2}{\varphi }_1\left(-{\tau }_2+s\right)\text{d}s}-{\displaystyle {\int }_0^{{\tau }_1^{\text{*}}}{\varphi }_1\left(-{\tau }_2^{\text{*}}+s\right)\text{d}s}\right)\end{array}\right)\end{array}$ (3.4)

We also express the nonlinear part as the following

$\begin{array}{l}F\left({\varphi }_1,{\varphi }_2\right)=\left(\begin{array}{c}-\frac{1}{3}{c}_1{\varphi }_1{\left(-{\tau }_1^{\text{*}}\right)}^3-\frac{1}{3}{a}_{12}{c}_2{\displaystyle {\int }_0^{{\tau }_2^{\text{*}}}{\varphi }_2{\left(-{\tau }_2^{\text{*}}+s\right)}^3\text{d}s}\\ -\frac{1}{3}{a}_{21}{\displaystyle {\int }_0^{{\tau }_2^{\text{*}}}{\varphi }_1{\left(-{\tau }_2^{\text{*}}+s\right)}^3\text{d}s}-\frac{1}{3}{c}_2{\varphi }_2{\left(-{\tau }_1^{\text{*}}\right)}^3\end{array}\right)\hfill \end{array}$ (3.5)

According to the property of DDEs, solution of Equation (3.1) is continuous on Banach space $C=C\left(\left[-\tau ,0\right],R^2\right)$ . With few discontinuity jumps, the solution operator is differentiable on the extended space $BC=C\left(\left[-\tau ,0\right],R^2\right)$ . We define solution with its domain as

$\left\{\varphi \in C\left(\left[-\tau ,0\right],R^2\right),{\varphi }^{\prime }\left(\theta \right)\in BC,\mathcal{A}\varphi ={\varphi }^{\prime }\left(\theta \right)+X_0\phi \right\}$

with $X_0$ is the fundamental solution matrix.

Consider the linear operator (3.2), it’s an infinitesimal generator of the strong continuous semigroup in phase space $BC$ , and we define the new linear opeartor $\mathcal{A}\left(ϵ\right):BC\to BC$ and its adjoint operator $A^{\text{*}}\left(ϵ\right):B{C}^{\text{*}}\to B{C}^{\text{*}}$ wherein $B{C}^{\text{*}}=C\left(\left[0,\tau \right],R^2\right)$ , that is

$\mathcal{A}\left(ϵ\right)\varphi =\{\begin{array}{ll}{\varphi }^{\prime }\left(\theta \right),\hfill & \text{for}-\tau \le \theta <0,\hfill \\ ℒ\left(0\right)\varphi +X_0\left({\varphi }^{\prime }\left(0\right)-ℒ\left(0\right)\varphi \right),\hfill & \text{for}\theta =0\hfill \end{array}$ (3.6)

and

${\mathcal{A}}^{\text{*}}\left(0\right)\psi =\{\begin{array}{ll}-{\psi }^{\prime }\left(-s\right),\hfill & \text{for}0\le s\le \tau ,\hfill \\ {ℒ}^{\text{*}}\left(0\right)\psi ,\hfill & \text{for}s=0\hfill \end{array}$ (3.7)

Based on Reize theorem, there exists the bounded variation matrix $\eta \left(\theta \right)$ to represent

$ℒ\left(0\right)\varphi ={\displaystyle {\int }_{-\tau }^0\text{d}\eta \left(\theta \right)\varphi \left(\theta \right)}$ (3.8)

and

${ℒ}^{\text{*}}\left(0\right)\psi ={\displaystyle {\int }_{-\tau }^0\text{d}{\eta }^{\text{T}}\left(-s\right)\psi \left(-s\right)}$ (3.9)

For $\varphi \in C\left(\left[-\tau ,0\right],R^2\right)$ , $\psi \in C\left(\left[0,\tau \right],R^2\right)$ , we define the inner product by its bilinear form

$\langle \psi ,\varphi \rangle ={\overline{\psi }}^{\text{T}}\left(0\right)\varphi \left(0\right)-{\displaystyle {\int }_{-\tau }^0{\overline{\psi }}^{\text{T}}\left(\xi -\theta \right)\text{d}\eta \left(\theta \right)\varphi \left(xi\right)\text{d}\xi }$ (3.10)

Set $U={\left(x,y\right)}^{\text{T}}$ , for $-\tau \le \theta \le 0$ , with definition $U\left(t+\theta \right)=U_t\left(\theta \right)$ , Equation (3.1) can be written into its differential operator form

${U^{\prime }}_t\left(\theta \right)=\mathcal{A}{U}_t\left(\theta \right)+ℛU_t\left(\theta \right),\text{for}-\tau \le \theta \le 0$ (3.11)

The reduction technique by using Schimdt-Lyapunov method is to project the solution onto the center manifold. Considering the linear version of differential operator (3.3), Hopf bifurcation occurs at $ϵ=0$ , and the associated characteristic roots set is finite which is denoted as $\Lambda =\left\{i\omega ,-i\omega \right\}$ as verified in Section 2. With the assumption of other eigenvalues having negative real parts, the phase space can be appended onto its center manifold. Hence the eigenspace is decomposed into the center subspace $P$ associated with eigenvalues of $\text{Λ}$ and its complementary space is represented by $Q$ . We suppose the eigenspace $P$ is spanned by $P=span\left\{q\left(\theta \right),\overline{q}\left(\theta \right)\right\}$ , $-\tau \le \theta \le 0$ ,where eigen vector $q\left(\theta \right)$ and its conjugate vector $\overline{q}\left(\theta \right)$ are respectively satisfied

$\mathcal{A}q\left(\theta \right)=i\omega q\left(\theta \right),\mathcal{A}\overline{q}\left(\theta \right)=-i\omega \overline{q}\left(\theta \right)$ (3.12)

We represent the eigenbasis $\Phi \left(\theta \right)=\left(q\left(\theta \right),\overline{q}\left(\theta \right)\right)$ , and correspondingly, the eigenbasis $\Psi \left(s\right)$ of the conjugate linear operator ${\mathcal{A}}^{\text{*}}\left(0\right)$ is denoted as $\Psi \left(s\right)=\left(p\left(s\right),\overline{p}\left(s\right)\right)$ given that

${\mathcal{A}}^{*}\left(0\right)p\left(s\right)=i\omega p\left(s\right),{\mathcal{A}}^{*}\left(0\right)\overline{p}\left(s\right)=-i\omega \overline{p}\left(s\right)$ (3.13)

We also make the equality $\langle \Phi \left(\theta \right),\Psi \left(s\right)\rangle =I$ satisfied.

Based on the Lyapunov-Schmidt reduction technique, the solution of the differential operator equation $U_t$ is decomposed into the direct summation of center space and its complementary space, that is

$U_t=zq\left(\theta \right)+\overline{z}\overline{q}\left(\theta \right)+y_t$ (3.14)

since $C=P\oplus Q$ . Therefore, by defining the projection operator $\Pi :C\to P$ , we have $\Pi U_t=\langle \Psi ,U_t\rangle \Phi$ , to get

$Z^{\prime }=BZ+{\overline{\Psi }}^{\text{T}}\left(0\right)X_0\left({\varphi }^{\prime }\left(0\right)-L\left(0\right)\varphi \right)+{\overline{\Psi }}^{\text{T}}\left(0\right)X_0ℛ\left(\varphi \right)$ (3.15)

with $B=\left(\begin{array}{cc}i\omega & 0\\ 0& -i\omega \end{array}\right)$ and $Z={\left(z_1,z_2\right)}^{\text{T}}$ .

We also rewrite the linear system of (3.13) into

$\begin{array}{c}{Z}^{\prime }=BZ+c_1{\overline{\Psi }}^{\text{T}}\left(0\right)\left(\begin{array}{c}{\varphi }_1\left(-{\tau }_1\right)-{\varphi }_1\left(-{\tau }_1^{*}\right)\\ {\varphi }_2\left(-{\tau }_1\right)-{\varphi }_2\left(-{\tau }_1^{*}\right)\end{array}\right)\\ +{\overline{\Psi }}^{\text{T}}\left(0\right)\left(\begin{array}{c}{a}_{12}{c}_2\left({\displaystyle {\int }_0^{{\tau }_2}{\varphi }_2\left(-{\tau }_2+s\right)\text{d}s}-{\displaystyle {\int }_0^{{\tau }_2^{\text{*}}}{\varphi }_2\left(-{\tau }_2^{\text{*}}+s\right)\text{d}s}\right)\\ a_{21}{c}_1\left({\displaystyle {\int }_0^{{\tau }_2}{\varphi }_1\left(-{\tau }_2+s\right)\text{d}s}-{\displaystyle {\int }_0^{{\tau }_1^{\text{*}}}{\varphi }_1\left(-{\tau }_2^{\text{*}}+s\right)\text{d}s}\right)\end{array}\right)\end{array}$ (3.16)

For example, we choose the eigenvectors

$q\left(\theta \right)=\left(\begin{array}{c}\frac{i}{\omega }{a}_{12}{c}_2\left({\text{e}}^{-i\omega {\tau }_2}-1\right)\\ b_1+i\omega -c_1{\text{e}}^{-i\omega {\tau }_1}\end{array}\right){\text{e}}^{i\omega \theta },p\left(s\right)=\left(\begin{array}{c}-\frac{i}{\omega }{a}_{21}{c}_1\left({\text{e}}^{i\omega {\tau }_2}-1\right)\\ b_1-i\omega -c_1{\text{e}}^{i\omega {\tau }_1}\end{array}\right){\text{e}}^{-i\omega s}$ (3.17)

for $-\tau \le \theta \le 0$ , $0\le s\le \tau$ . The linear system (3.14) is written into

$z^{\prime }\left(t\right)=i\omega z-zϵ\left(i{\tau }_m\left({\overline{p}}_1{q}_1+{\overline{p}}_2{q}_2\right)c_1\omega {\text{e}}^{-i\omega {\tau }_1}+{\tau }_n\left({\text{e}}^{i\omega {\tau }_2}-1\right)\left(a_{12}{c}_2{\overline{p}}_1{q}_2+a_{21}{c}_1{\overline{p}}_2{q}_1\right)\right)$(3.18)

and the nonlinear system is written as

$\begin{array}{c}{z}^{\prime }\left(t\right)=i\omega z+{\overline{p}}_1(-\frac{1}{3}{c}_1{\left(z{q}_1{\text{e}}^{-i\omega {\tau }_1^{\text{*}}}+\overline{z}{\overline{q}}_1{\text{e}}^{i\omega {\tau }_1^{\text{*}}}\right)}^3\\ -\frac{1}{3}{a}_{12}{c}_2{\displaystyle {\int }_0^{{\tau }_2^{\text{*}}}{\left(z{q}_2{\text{e}}^{-i\omega \left({\tau }_2^{\text{*}}-s\right)}+\overline{z}{\overline{q}}_2{\text{e}}^{i\omega \left({\tau }_2^{\text{*}}-s\right)}\right)}^3\text{d}s})\\ +{\overline{p}}_2(-\frac{1}{3}{a}_{21}{\displaystyle {\int }_0^{{\tau }_2^{\text{*}}}{\left(z{q}_1{\text{e}}^{-i\omega \left({\tau }_2^{\text{*}}-s\right)}-\overline{z}{\overline{q}}_1{\text{e}}^{i\omega \left({\tau }_2^{\text{*}}-s\right)}\right)}^3\text{d}s}\\ -\frac{1}{3}{c}_2{\left(z{q}_2{\text{e}}^{-i\omega {\tau }_1^{\text{*}}}+\overline{z}{\overline{q}}_2{\text{e}}^{i\omega {\tau }_1^{\text{*}}}\right)}^3)\\ =i\omega z-\left({\overline{p}}_1{c}_1{q}_1^2{\overline{q}}_1+{\overline{p}}_2{c}_2{q}_2^2{\overline{q}}_2\right){\text{e}}^{-i\omega {\tau }_1}+\frac{i}{\omega }\left({\overline{p}}_1{a}_{12}{c}_2{\overline{q}}_2{q}_2^2+{\overline{p}}_2{a}_{21}{c}_1{\overline{q}}_1{q}_1^2\right)\\ -\frac{i}{\omega }\left({\overline{p}}_1{a}_{12}{c}_2{q}_2^2{\overline{q}}_2+a_{21}{\overline{p}}_2{\overline{q}}_1{c}_1{q}_1^2\right){\text{e}}^{-i\omega {\tau }_2}\\ =i\omega z+f_{21}{z}^2\overline{z}\end{array}$ (3.19)

We also have

${y^{\prime }}_t\left(\theta \right)={U^{\prime }}_t\left(\theta \right)-<\Psi ,{U^{\prime }}_t\left(\theta \right)\Phi \left(\theta \right)=\mathcal{A}{y}_t+ℛ\left(\varphi \right)-{\overline{\Psi }}^{\text{T}}\left(0\right)X_0ℛ\left(\varphi \right)$ (3.20)

With the above discussion, we have the dimensional reduction system

$z^{\prime }\left(t\right)=i\omega z+ϵ\alpha z+f_{21}{z}^2\overline{z}$ (3.21)

with

$\begin{array}{l}\alpha =-\left(i{\tau }_m\left({\overline{p}}_1{q}_1+{\overline{p}}_2{q}_2\right)c_1\omega {\text{e}}^{-i\omega {\tau }_1}+{\tau }_n\left({\text{e}}^{i\omega {\tau }_2}-1\right)\left(a_{12}{c}_2{\overline{p}}_1{q}_2+a_{21}{c}_1{\overline{p}}_2{q}_1\right)\right)\\ f_{21}=-\left({\overline{p}}_1{c}_1{q}_1^2{\overline{q}}_1+{\overline{p}}_2{c}_2{q}_2^2{\overline{q}}_2\right){\text{e}}^{-i\omega {\tau }_1}+\frac{i}{\omega }\left({\overline{p}}_1{a}_{12}{c}_2{\overline{q}}_2{q}_2^2+{\overline{p}}_2{a}_{21}{c}_1{\overline{q}}_1{q}_1^2\right)\\ -\frac{i}{\omega }\left({\overline{p}}_1{a}_{12}{c}_2{q}_2^2{\overline{q}}_2+a_{21}{\overline{p}}_2{\overline{q}}_1{c}_1{q}_1^2\right){\text{e}}^{-i\omega {\tau }_2}\end{array}$

with $a_{12}=1$ , $a_{21}=1.2$ , $c_1=-2$ , $c_2=0.6023012058$ , $\omega =1.7197$ , ${\tau }_1=0.7850148049$ , ${\tau }_2=1.141869341$ , $b_1=0.1$ , $b_2=0.1$ , we compute $\alpha =\left(1.775180730+3.921586787i\right)ϵ{\tau }_m+\left(-1.465601884+1.222555424i\right){\tau }_n$ , and $l_1\left(0\right)=-0.1949329832+0.4889169011\text{e}-1i$ . Henceafter, the stable periodical solution bifurcates from Hopf point. The continued periodical solutions by varying two-time delays is shown in Figure 3.

Figure 3. The continuous of periodical solution as varying parameter ${\tau }_1$ and ${\tau }_2$ . (a) time delay steps ${\tau }_m=0.02$ , ${\tau }_n=-0.008$ ; (b) time delay steps ${\tau }_m=0.02$ , ${\tau }_n=0.008$ .

For system (1.3), we also simulate torus solution by increasing timed delay ${\tau }_2$ , as shown in Figure 4. The asymmetry solution of system (1.3) is computed by DDE23 software. The chosen parameter are listed as ${\tau }_1=1.9$ and ${\tau }_2=13.9,14.9,15.9,16.9,17.9,18.9$ respectively.

4. Numerical Simulation

As discussed in Section 2, Hopf bifurcation occurs with the corresponding transversal condition being satisfied, and the periodical oscillation solution arise. We simulate the periodical solutions both in time series solution and phase portraits as system (1.2) loss stability underlying Hopf bifurcation. Some novel periodical solutions with spatial symmetry are simulated which inspire us the enthusiasm to pay attention to the continuation of periodical solutions further. Whether the spatial symmetry property can be continued with doubly period or quadruple period? DDE-Biftool is artificial mathematical software that analyzes system bifurcation behavior. Herein we apply DDE-Biftool to continue periodical solutions of system (1.2) with multiple period. Firstly, we use DDE23 to simulate the spatial symmetrical solution with doubly and quadruple period, which can often be progressed with Runga-Kutta algorithm too. Both the time series solution and phase portraits of symmetrical solutions are shown in Figure 3. With fixed parameter $a_{12}=1$ ,

Figure 4. The phase portraits of strange attractors in system (1.1) with time delays varying. (a) ${\tau }_1=1.9$ , ${\tau }_2=13.9$ ; (b) ${\tau }_1=1.9$ , ${\tau }_2=14.9$ ; (c) ${\tau }_1=1.9$ , ${\tau }_2=15.9$ ; (d) ${\tau }_1=1.9$ , ${\tau }_2=16.9$ ; (e) ${\tau }_1=1.9$ , ${\tau }_2=17.9$ ; (f) ${\tau }_1=1.9$ , ${\tau }_2=18.9$ .

$a_{21}=1.2$ , $c_1=-2$ , $c_2=1.2$ and time delays ${\tau }_1=1.9$ , ${\tau }_2=2.2$ , the doubly period solutions are observed, however the quadruple solution is found with ${\tau }_1=1.3$ , ${\tau }_2=2.8$ .

Figure 5. The P2 and P4 solution with spatial symmetry are observed with time delays ${\tau }_1=1.9$ , ${\tau }_2=2.2$ and ${\tau }_1=1.3$ , ${\tau }_2=2.8$ respectively. (a) Time series solution of P2 solution is simulated; (b) Time series solution of P4 solution; (c) The phase portraits of P2 solution in (a); (d) The phase portraits of P4 solution in (b).

Using DDE-Biftool, the continuation work of period solutions as varying free parameter continuously is completed by br_contn command, which manifests the solution mind branches for ever. The question is answered, that the doubly period solution and quadruple period solution observed in Figure 4 can be continued with spatial symmetry while varying free parameter $c_1$ or $c_2$ . The track of period solution branches expressed by maximal magnitude forms a circle as varying $c_1$ and $c_2$ respectively, as shown in Figure 5. The blue circles shown in Figure 4(a) and Figure 4(b) represent the continued solution branches with doubly period, whilst the red circles are the quadruple period solution continued

Figure 6. The continuation of P2 solutions and P4 solutions as varying free parameter. (a) The blue color circle of continued P2 solutions as varying $c_1$ , with $c_2=1.2$ , ${\tau }_1=1.9$ , ${\tau }_2=2.2$ ; The red color circle of continued P4 solutions as varying $c_1$ , with $c_2=1.2$ , ${\tau }_1=1.3$ , ${\tau }_2=2.8$ ; (b) The blue color circle of continued P2 solutions as varying $c_2$ , with $c_1=-2$ , ${\tau }_1=1.9$ , ${\tau }_2=2.2$ ; The red color circle of continued P4 solutions as varying $c_2$ , with $c_1=-2$ , ${\tau }_1=1.3$ , ${\tau }_2=2.8$ ; (c) The corresponding period of P2 solutions and P4 solutions (as shown in (a)) versus $c_1$ ; (d) The corresponding period of P2 solutions and P4 solutions (as shown in (b)) versus $c_2$ .

branches. The corresponding period versus $c_1$ and $c_2$ respectively are shown in Figure 4(c) and Figure 4(d). The oscillation rhythm keeps its symmetry property, which is reversal by its time series solutions and named as P2 solutions and P4 solutions respectively. The continued P2 solutions and P4 solutions manifest spatial symmetry, as shown in Figure 5. The blue color pictures are periodical solutions continued by varying parameter $c_1$ , as shown in Figure 6(a) and Figure 6(b), whilst the green color pictures are produced by continued periodical solutions with $c_2$ free parameter, as shown in Figure 5(c) and Figure 5(d). It seems that some islands of P2 solutions continued circles peered in eyesight if varying $c_1$ or $c_2$ parameter. We endeavor to simulate P2 solutions and it is feasible as varying $c_1$ and $c_2$ , which formes three circles with maximal magnitudes versus free parameter. It is noted that with fixed parameter $b_1=0.1$ , $b_2=0.1$ , $a_{12}=1$ and $a_{21}=1.2$ , and time delay being ${\tau }_1=1.9$ , ${\tau }_2=2.2$ , free parameter is listed as three arrays, alike $c_1=-2$ , $c_2=1.2$ and $c_1=-4$ , $c_2=2.4$ and $c_1=-6$ , $c_2=3.6$ , and the continuation of P2 solution is a noval work since the spatial symmetry property of periodical solutions preserved further. As shown in Figure 6, the island circles are drawn respectively with $c_1=-2$ , $c_1=-4$ and $c_1=-6$ respectively, as varying $c_2$ , the P2 solutions with spatial symmetry is observed and continued by using brcontn programm in DDE-Biftool. The island of period circles versus $c_2$ free parameter are usually observed as shown in Figure 7(b) too.

Figure 7. The continuous P2 and P4 periodical solutions with varying free parameter. (a) P2 solutions continued as varying $c_1$ ; (b) P4 solutions as varying $c_2$ ; (c) P2 solutions obtained as varying $c_2$ ; (d) P4 solutions simulated as varying $c_2$ . The numerical simulation method are using DDE-Biftool software.

For system (1.1), the periodical oscillation happens undertaken Hopf bifurcation at $\left({\tau }_1,{\tau }_2\right)$ point as discussed in Section 2. However, the strange attractor is focused in system as equilibrium loss its stability. We simulate strange attractor with fixed parameter $b_1=1$ , $b_2=1$ , $a_{12}=0.1$ , $a_{21}=0.12$ , $c_1=-2$ , $c_2=1$ . Then varying time delays ${\tau }_1$ and $ta{u}_2$ with a bigger ${\tau }_2$ assumed, the strange attractors are observed with the phase portraits are shown in Figures 8(a)-(e).

Figure 8. The continuous circls with P2 solutions of spatial symmetry found. (a) Three P2 circles produced as varying magnitudes using continuous method in DDE-Biftool; (b) Three P2 circles formed with period v.s. $c_2$ .

Figure 9. BP bifurcation with zero equilibrium solution as varying parameter $c_1$ , with other parameters $c_2=5$ , ${\tau }_1=1.6$ , ${\tau }_2=0.8$ . (a) BP bifurcation happens at $c_1=-0.4371$ ; the red color curves denote the unstable equilibrium solutions, whilst the blue color curves represent the stable equilibria; (b) The phase picture of continuous of periodical solutions arising from Hopf bifurcation point at zero solution with $c_1=-4.435$ ; (c) The phase picture of continuous of periodical solutions as varying parameter $c_1$ , which arise from equilibrium $E_0=\left(5.9663,0.7584\right)$ at $c_1=-2.592$ ; (d) The continuous of periodical solutions arise from equilibrium $E_0=\left(-5.9663,-0.7584\right)$ at $c_1=-2.592$ .

Branch point bifurcation phenomena occurs at $c_1=-0.4371$ . As shown in Figure 9, there exists a pair of equilibria which is symmetry about zero point in the plane at $c_1=-0.4371$ . Henceafter, the branches of equilibria arising from BP points are shown in Figure 9(a), which is unstable at red line whilst stable in green color. Notice Hopf bifurcation occurs $c_1=-2.592$ , and the bifurcating periodical solutions are arise from Hopf points. As shown in Figure 9(c) and Figure 9(d), the bifurcating solutions are symmetry. With trivial solution, Hopf bifurcation occurs at $c_1=-4.443$ , and the periodical solutions is continued as shown in Figure 9(b).

5. Conclusion

The Hopfield network model is the artificial neuron biological model which relates to the peoples associative memory alike with synaptic connections of neurons. By distributed time delay or discrete time delay, the information feedback between neurons was significantly effective with its past time history. For the distributed delay system, Hopf bifurcation happened as the characteristic roots crossed over the imaginary axis from Hopf left plane to right half plane. Applying center manifold theorem, the normal form was usually computed by writing ODEs system on extended phase space, hence the solution operator was projected onto the center subspace by dimension reduction technique. The stability of bifurcating periodical solutions was calculated and the bifurcating direction resulted from the linearized system. The bifurcating periodical solutions were simulated numerically by their equivalent DDEs, which were composed of four equations. For the discrete time delay system, DDE-Biftool software is applied to investigate the period-doubling bifurcation of the continued limit cycles further. The continuation of P2 and P4 periodical solutions was done. The spatial symmetry P2 and P4 solutions were continued by circles by varying free parameters.

Availability of Data and Materials

All data generated or analyzed during this study are included in this published article.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Hopfield, J.J. (1982) Neural Networks and Physical Systems with Emergent Collective Computational Abilities. Proceedings of the National Academy of Sciences, 79, 2554-2558. https://doi.org/10.1073/pnas.79.8.2554
[2]	Hopfield, J.J. (1984) Neurons with Graded Response Have Collective Computational Properties Like Those of Two-State Neurons. Proceedings of the National Academy of Sciences, 81, 3088-3092. https://doi.org/10.1073/pnas.81.10.3088
[3]	Korner, E., Kupper, R., Rahman, M.K.M. and Shkuro, Y. (2007) Neurocomputing Research Developments. Nova Science Publishers.
[4]	Skarda, C.A. and Freeman, W.J. (1987) How Brains Make Chaos in Order to Make Sense of the World. Behavioral and Brain Sciences, 10, 161-173. https://doi.org/10.1017/s0140525x00047336
[5]	Skarda, C.A. and Freeman, W.J. (1990) Chaos and the New Science of the Brain. Neuroscience, 1, 275-285.
[6]	Kundu, A., Das, P. and Roy, A.B. (2013) Complex Dynamics of a Four Neuron Network Model Having a Pair of Short-Cut Connections with Multiple Delays. Nonlinear Dynamics, 72, 643-662. https://doi.org/10.1007/s11071-012-0742-2
[7]	Das, A., Das, P. and Roy, A.B. (2002) Chaos in a Three-Dimensional General Model of Neural Network. International Journal of Bifurcation and Chaos, 12, 2271-2281. https://doi.org/10.1142/s0218127402005820
[8]	Lewis, J.E. and Glass, L. (1991) Steady States, Limit Cycles, and Chaos in Models of Complex Biological Networks. International Journal of Bifurcation and Chaos, 1, 477-483. https://doi.org/10.1142/s0218127491000373
[9]	Huang, Y. and Yang, X. (2006) Hyperchaos and Bifurcation in a New Class of Four-Dimensional Hopfield Neural Networks. Neurocomputing, 69, 1787-1795. https://doi.org/10.1016/j.neucom.2005.11.001
[10]	Yang, X. and Huang, Y. (2007) Chaos and Two-Tori in a New Family of 4-CNNS. International Journal of Bifurcation and Chaos, 17, 953-963. https://doi.org/10.1142/s0218127407017677
[11]	Guckenheimer, J. and Holmes, P. (1997) Nonlinear Oscillation, Dynamical Systems and Bifurcations of Vector Fields. Springer.
[12]	Rech, P.C. (2011) Dynamics of a Neuron Model in Different Two-Dimensional Parameter-Spaces. Physics Letters A, 375, 1461-1464. https://doi.org/10.1016/j.physleta.2011.02.037
[13]	Luonan Chen, and Aihara, K. (1999) Global Searching Ability of Chaotic Neural Networks. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, 46, 974-993. https://doi.org/10.1109/81.780378
[14]	Yang, X. and Yuan, Q. (2005) Chaos and Transient Chaos in Simple Hopfield Neural Networks. Neurocomputing, 69, 232-241. https://doi.org/10.1016/j.neucom.2005.06.005
[15]	Engelborghs, K., Luzyanina, T. and Samae, G. (2001) DDEBIFTOOL, a Matlab Package for Bifurcation Analysis of Delay Differential Equations, Technical Report TW330.
[16]	Sieber, J., Engelborghs, K., Samaey, G. and Roose, D. (2015) DDEBIFTOOL Manual Bifurcation Analysis of Delay Differential Equations. https://arxiv.org/abs/1406.7144
[17]	Ma, S. (2022) Bifurcation Analysis of Periodic Oscillation in a Hematopoietic Stem Cells Model with Time Delay Control. Mathematical Problems in Engineering, 2022, 1-10. https://doi.org/10.1155/2022/7304280
[18]	Hale, J. (2003) Theory of Functional Differential Equations. World Publishing Corporation.
[19]	Ma, S.Q, Lu, Q. and Hogan, S. (2007) The Normal Form for Double Hopf bifurcation for Stuart-Landau System with Nonlinear Delay Feedback and Delay-Dependent Parameters. Advances in Complex Systems, 10, 423-448. https://doi.org/10.1142/s0219525907001227