Branch switching for a Shil'nikov type homoclinic orbit in the FitzHugh-Nagumo equations.

The FitzHugh-Nagumo (FHN) equations [#!FitzH:61!#,#!NaArYo:62!#] are a simplified version of the Hodgkin-Huxley equations [#!HoHu:52!#]. They model nerve axon dynamics and are given by

$$\begin{split}u_t&=u_{xx}-f_a(u)-w, \\ w_t&=\epsilon(u-\gamma w), \end{split}$$ (27.2)

where

$$f_a(u) = u (u-a)(u-1).$$

Travelling wave solutions of the form $(u,w)(x,t)=(u,w)(\xi)$, where $\xi=x+ct$ are solutions of the following ODE system:

$$\begin{split}\dot u &= v,\\ \dot v &= c v + f_a(u) + w,\\ \dot w &= \frac{\epsilon}{c} (u - \gamma w). \end{split}$$ (27.3)

In particular we consider solitary wave solutions of (27.2). These correspond to orbits homoclinic to $(u,v,w)=0$ in system (27.3). In our numerical example we keep $\gamma=0$.

We aim to find a $2$-homoclinic orbit at a Shil'nikov bifurcation. All the commands given here are in the file fnb.auto. First we obtain a homoclinic orbit using a homotopy technique (see FrDoMo:94), using ISTART=3, for the parameter values $c=0.21, a=0.2, \epsilon=0.0025$.

demo('sib')
ld('fnb')
rn()
sv('1')

Among the output we see:

  BR    PT  TY  LAB     PERIOD       L2-NORM    ...   PAR(17)
   1    20  UZ    3   2.92257E+01   2.37916E-01 ... -1.68000E-09

and a zero of PAR(17) means that a zero of an artificial parameter has been located and the right-hand end point of the corresponding solution belongs to the plane that is tangent to the stable manifold at the saddle. This point still needs to come closer to the equilibrium, which we can achieve by further increasing the period to 300, while keeping PAR(17) at 0:

rn(c='fnb.2',h='fnb.1',s='1')
sv('2')

  BR    PT  TY  LAB     PERIOD       L2-NORM    ...   PAR(2)
   1   190  UZ   10   3.00000E+02   7.37932E-02 ...  1.79286E-01

Next we stop using the homotopy technique and increase the period even further, to 1000.

rn(c='fnb.3',h='fnb.3',s='2')
sv('3')

  BR    PT  TY  LAB     PERIOD       L2-NORM    ...   PAR(2)
   1    80  UZ   13   1.00000E+03   4.04183E-02 ...  1.79286E-01

A continuation in PAR(2)=$a$ and PAR(1)=$c$ needs to be performed to arrive at the place where we wish to find a 2-homoclinic orbit: $a=0$. At the same time we monitor PAR(22) to locate Belyakov points.

rn(c='fnb.4',h='fnb.4',s='3')
sv('4')

  BR    PT  TY  LAB    PAR(2)        L2-NORM    ...   PAR(1)        PAR(22)
   1     6  UZ   15   1.31812E-01   3.28710E-02 ...  2.17166E-01  -6.31243E-06
   1    23  UZ   19  -8.54548E-08   1.56158E-02 ...  2.74218E-01  -9.88772E-02

Hence, there exists a Belyakov point at $(a,c)=(0.131812,0.21766)$. At label 19 we have a lower value of $a$ than at the Belyakov point, and by inspection of the file d.4 we can observe that the equilibrium has one positive eigenvalue and a complex conjugate pair of eigenvalues with negative real part, and conclude that this orbit is of Shil'nikov type. Before starting the homoclinic branch switching, we calculate the adjoint to obtain a `Lin vector':

rn(c='fnb.5',h='fnb.5',s='4')
sv('5')

  BR    PT  TY  LAB    PAR(9)        L2-NORM    ...   PAR(3)     
   1     2  EP   28  -1.00000E+00   1.56158E-02 ...  2.50000E-03

Next, we continue in the time $T_1$ (PAR(21)), the gap $\varepsilon _1$ (PAR(22)) and $c$ (PAR(1)), and by setting ISTART=-2 we try to locate a 2-homoclinic orbit:

rn(c='fnb.6',h='fnb.6',s='5')
sv('6')

In fact we find many of them, exactly as is predicted by the theory:

  BR    PT  TY  LAB    PAR(21)    ...   PAR(1)        PAR(22)
... 
   1   175  UZ   46   1.64818E+02 ...  2.74218E-01   4.59218E-16
   1   180  UZ   47   1.44759E+02 ...  2.74218E-01  -1.43728E-14
   1   184  UZ   48   1.24939E+02 ...  2.74218E-01   1.55506E-13
   1   189  UZ   49   1.04615E+02 ...  2.74218E-01  -2.37665E-11
   1   193  UZ   50   8.53538E+01 ...  2.74218E-01   1.02165E-11
   1   198  UZ   51   6.37899E+01 ...  2.74218E-01  -5.74204E-14

Each of these homoclinic orbits differ by about 20 in the value $T_1$. This is about the time it takes to make one half-turn close to and around the equilibrium, so that orbits differ by the number of half turns around the equilibrium before a big excursion in phase space. Note that the variation of $c$ is so small that it does not appear.

A plot of $T_1$ vs. $\varepsilon _1$ gives insight into how the gap is opened and closed in the continuation process. This is depicted in Figure 27.5.

Figure 27.5: A plot of $\varepsilon _1$ as a function of $T_1$ during our computation of Shil'nikov-type two-homoclinic orbits. Each zero corresponds to a different orbit.

We are now in a position to continue each of these orbits as a normal homoclinic orbit by setting ISTART=1 and ITWIST=0. We leave this as an exercise to the reader.