Example: the fast Morris-Lecar equations

We consider the following system
$\displaystyle \left\{\begin{array}{rcl}
\dot v&=&y-0.5(v+0.5)-2w(v+0.7)-m_{\infty}(v-1)\\
\dot w&=&1.15(w_{\infty}-w)\tau
\end{array}\right.$     (74)

where $m_{\infty}(v)=(1+\tanh((v+0.01)/0.15))/2$, $w_{\infty}(v)=(1+\tanh((v-z)/0.145))/2$ and $\tau=\cosh((v-0.1)/0.29).$ Here $v$ and $w$ are the state variables and $y$ and $z$ are the parameters. This is a modification of the fast subsystem of the Morris-Lecar equations studied in [37],[38]; the Morris-Lecar equations were introduced in [33] as a model for the electrical activity in the barnacle giant muscle fiber. In our model y corresponds to the slow variable in the fast Morris-Lecar equations; z is the potential for which $w_{\infty}(z)=\frac{1}{2}$.

We find a stable equilibrium (EP) for $y=0.110472$ and $z=0.1$ at $(0.04722,0.32564)$ by time integration. We continue this equilibrium with free parameter $y$ for decreasing values of $y$ by running testEquilMLfast.m:

p=[0.11047;0.1];ap1=[1];
[x0,v0]=init_EP_EP(@MLfast,[0.047222;0.32564],p,ap1);
opt=contset;
opt=contset(opt,'Singularities',1);
opt=contset(opt,'MaxNumPoints',65);
opt=contset(opt,'MinStepsize',0.00001);
opt=contset(opt,'MaxStepsize',0.01);
opt=contset(opt,'Backward',1);

[x,v,s,h,f]=cont(@equilibrium,x0,[],opt);
cpl(x,v,s,[3 1]);

The output is as follows:

>> testEquilMLfast
first point found
tangent vector to first point found
label = H , x = ( 0.036756 0.294770 0.075658 )
First Lyapunov coefficient = 8.238955e+00
label = LP, x = ( -0.033738 0.136501 -0.020727 )
a=-1.036700e+01
Neutral saddle
label = H , x = ( -0.119894 0.045956 0.033207 )
label = LP, x = ( -0.244914 0.008514 0.083257 )
a=-2.697425e+00

elapsed time  = 0.1 secs
npoints curve = 65

Figure 20: Computed equilibrium curve in the MLfast model
\includegraphics[scale=0.6]{ex/LPLC1.eps}
We find a Hopf (H) bifurcation point at $y=0.075659$, two limit points (LP) at $y=-0.020727$ and $y=0.083257$ and a neutral saddle (H) at $y=0.033207$. There are stable equilibria before the first H point and after the second LP point and unstable equilibria between the first H point and the second LP point. The Lyapunov coefficient in the first Hopf point $l_1=8.234573$ is positive which means that the periodic orbits are born unstable. The results are plotted using the plot function cpl, cf. §3.9 where the fourth argument is used to select the third and first component of the solution which are the parameter $y$ and the coordinate $v$. The results can be seen in Figure 20.

The detected Hopf point is used to start a limit cycle continuation. We choose $N=30$ test intervals and $m=4$ collocation points for the discretization.

testEquilMLfast
x1=x(1:2,s(2).index);p=[x(end,s(2).index);0.1];
[x0,v0]=init_H_LC(@MLfast,x1,p,ap1,0.0001,30,4);
opt=contset;
opt=contset(opt,'IgnoreSingularity',1);
opt=contset(opt,'Singularities',1);
opt=contset(opt,'MaxNumPoints',50);
[x2,v2,s2,h2,f2]=cont(@limitcycle,x0,v0,opt);

plotcycle(x2,v2,s2,[1 2]);
The output is as follows:
>> testLCMLfast
first point found
tangent vector to first point found
label = H , x = ( 0.036756 0.294770 0.075658 )
First Lyapunov coefficient = 8.238955e+00
label = LP, x = ( -0.033738 0.136501 -0.020727 )
a=-1.036700e+01
label = H , x = ( -0.119894 0.045956 0.033207 )
label = LP, x = ( -0.244914 0.008514 0.083257 )
a=-2.697425e+00

elapsed time  = 0.1 secs
npoints curve = 65
first point found
tangent vector to first point found
Limit point cycle (period = 4.222012e+00, parameter = 8.456948e-02)
Normal form coefficient = -2.334578e-01

elapsed time  = 4.0 secs
npoints curve = 50

The periodic orbit is initially unstable. We detect a limit point of cycles LPC at $y=0.084569$. At this point the stability is gained. Afterwards the stability is preserved but the period tends to infinity and the periodic orbits end in a homoclinic orbit. The results can be seen in Figure 21.

Figure 21: Computed limit cycle curve started from a Hopf bifurcation
\includegraphics[scale=0.7]{ex/LPLC2.eps}

We now compute a curve of fold bifurcations of limit cycles. The starting vector x0 is calculated from the LPC on the previously computed curve of limit cycles, using init_LPC_LPC. Continuation is done using a call to the standard continuer with limitpointcycle as curve definition file. We free both $y$ and $z$ to continue the LPC curve through this LPC point. The computations are done by executing the script testLPCMLfast.m:

testEquilMLfast
x1=x(1:2,s(2).index);p=[x(end,s(2).index);0.1];
[x0,v0]=init_H_LC(@MLfast,x1,p,ap1,0.0001,30,4);
opt=contset;
opt=contset(opt,'IgnoreSingularity',1);
opt=contset(opt,'Singularities',1);
opt=contset(opt,'MaxNumPoints',50);
opt=contset(opt,'FunTolerance',0.0000001);
opt=contset(opt,'VarTolerance',0.0000001);
[x2,v2,s2,h2,f2]=cont(@limitcycle,x0,v0,opt);
[x0,v0]=init_LPC_LPC(@MLfast,x2,s2(2),[1 2],30,4);
opt=contset;
opt=contset(opt,'FunTolerance',0.0001);
opt=contset(opt,'VarTolerance',0.0001);
opt=contset(opt,'MaxNumPoints',30);
%opt=contset(opt,'Backward',1);
opt=contset(opt,'Singularities',1);
[x3,v3,s3,h3,f3]=cont(@limitpointcycle,x0,v0,opt);
plotcycle(x3,v3,s3,[1 2]);

The output is as follows:

>> testLPCMLfast
first point found
tangent vector to first point found
label = H , x = ( 0.036756 0.294770 0.075658 )
First Lyapunov coefficient = 8.238955e+00
label = LP, x = ( -0.033738 0.136501 -0.020727 )
a=-1.036700e+01
Neutral saddle
label = H , x = ( -0.119894 0.045956 0.033207 )
label = LP, x = ( -0.244914 0.008514 0.083257 )
a=-2.697425e+00

elapsed time  = 0.4 secs
npoints curve = 65
first point found
tangent vector to first point found
Limit point cycle (period = 4.222012e+00, parameter = 8.456948e-02)
Normal form coefficient = -2.334578e-01

elapsed time  = 3.9 secs
npoints curve = 50
first point found
tangent vector to first point found

elapsed time  = 7.1 secs
npoints curve = 30

The results are plotted using the standard plot function plotcycle where the fourth argument is used to select the coordinates. The results can be seen in Figure 22. We note that it shrinks to a single point. The labels of the plot are added manually .

Figure 22: Computed fold of limit cycles curve started from a fold bifurcation of limitcycles
\includegraphics[scale=0.8]{ex/Figure22.eps}