(Another) Neuron model (codim 2, periodic orbits)

• (Another) Neuron model (codim 2, periodic orbits)
• • Continuation and codim 1 bifurcations
  • Normal forms computation
  • Continuation of Hopf points
  • Continuation of Fold points

Consider the neuron model

\[\begin{aligned} & \dot{x}_1(t)=-x_1(t)-a g\left(b x_1\left(t-\tau_1\right)\right)+c g\left(d x_2\left(t-\tau_2\right)\right) \\ & \dot{x}_2(t)=-x_2(t)-a g\left(b x_2\left(t-\tau_1\right)\right)+c g\left(d x_1\left(t-\tau_2\right)\right) \end{aligned}\]

where $g(z)=[\tanh (z-1)+\tanh (1)] \cosh (1)^2$.

Continuation and codim 1 bifurcations

We first instantiate the model

using Revise, DDEBifurcationKit, Parameters, Plots
using BifurcationKit
const BK = BifurcationKit

g(z) = (tanh(z − 1) + tanh(1))*cosh(1)^2
function neuron2VF(x, xd, p)
   @unpack a,b,c,d = p
   [
      -x[1] - a * g(b*xd[1][1]) + c * g(d*xd[2][2]),
      -x[2] - a * g(b*xd[1][2]) + c * g(d*xd[2][1])
   ]
end

function delaysF(par)
   [par.τ1, par.τ2]
end

pars = (a = 0.25, b = 2., c = 15/29, d = 1.2, τ1 = 12.7, τ2 = 20.2)
x0 = [0.01, 0.001]

prob = ConstantDDEBifProblem(neuron2VF, delaysF, x0, pars, (@lens _.a))

optn = NewtonPar(verbose = false, eigsolver = DDE_DefaultEig(maxit=100))
opts = ContinuationPar(p_max = 1., p_min = 0., newton_options = optn, ds = 0.01, detect_bifurcation = 3, nev = 9, dsmax = 0.2, n_inversion = 4)
br = continuation(prob, PALC(tangent=Bordered()), opts)

 ┌─ Curve type: EquilibriumCont
 ├─ Number of points: 10
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter a starts at 0.25, ends at 1.0
 ├─ Algo: PALC
 └─ Special points:

If `br` is the name of the branch,
ind_ev = index of the bifurcating eigenvalue e.g. `br.eig[idx].eigenvals[ind_ev]`

- #  1,     hopf at a ≈ +0.29890362 ∈ (+0.29853073, +0.29890362), |δp|=4e-04, [converged], δ = ( 2,  2), step =   3, eigenelements in eig[  4], ind_ev =   2
- #  2,     hopf at a ≈ +0.34537483 ∈ (+0.34523499, +0.34537483), |δp|=1e-04, [converged], δ = ( 2,  2), step =   4, eigenelements in eig[  5], ind_ev =   4
- #  3,     hopf at a ≈ +0.36830746 ∈ (+0.36662946, +0.36830746), |δp|=2e-03, [converged], δ = ( 2,  2), step =   5, eigenelements in eig[  6], ind_ev =   6
- #  4,     hopf at a ≈ +0.57449141 ∈ (+0.57260366, +0.57449141), |δp|=2e-03, [converged], δ = ( 2,  2), step =   7, eigenelements in eig[  8], ind_ev =   8
- #  5,     hopf at a ≈ +1.00000000 ∈ (+1.00000000, +1.00000000), |δp|=0e+00, [    guess], δ = ( 2,  2), step =   9, eigenelements in eig[ 10], ind_ev =  10
- #  6, endpoint at a ≈ +1.00000000,                                                                      step =   9

We then plot the branch

scene = plot(br)

Normal forms computation

As in BifurcationKit.jl, it is straightforward to compute the normal forms.

hopfpt = BK.get_normal_form(br, 2)

SubCritical - Hopf bifurcation point at a ≈ 0.34537482504934847.
Frequency ω ≈ 0.1826481323494928
Period of the periodic orbit ≈ 34.40049031082816
Normal form z⋅(iω + a⋅δp + b⋅|z|²):
┌─ a = 0.07653825207812184 + 0.08333988803350242im
└─ b = 0.07985677803156456 + 0.001883464762243053im

Continuation of Hopf points

We follow the Hopf points in the parameter plane $(a,c)$. We tell the solver to consider br.specialpoint[1] and continue it.

# continuation of the first Hopf point
brhopf = continuation(br, 1, (@lens _.c),
         setproperties(br.contparams, detect_bifurcation = 1, dsmax = 0.01, max_steps = 100, p_max = 1.1, p_min = -0.1,ds = 0.01, n_inversion = 2);
         verbosity = 0,
         detect_codim2_bifurcation = 2,
         bothside = true,
         start_with_eigen = true)

brhopf2 = continuation(br, 2, (@lens _.c),
         setproperties(br.contparams, detect_bifurcation = 1, dsmax = 0.01, max_steps = 100, p_max = 1.1, p_min = -0.1,ds = -0.01);
         verbosity = 0,
         detect_codim2_bifurcation = 2,
         bothside = true,
         start_with_eigen = true)

scene = plot(brhopf, vars = (:a, :c), xlims = (0,0.7), ylims = (0,1))
plot!(scene, brhopf2, vars = (:a, :c), xlims = (-0,0.7), ylims = (-0.1,1))
scene

Continuation of Fold points

We follow the Fold points in the parameter plane $(a, c)$. We tell the solver to consider br2.specialpoint[3] and continue it.

prob2 = ConstantDDEBifProblem(neuron2VF, delaysF, x0, (@set pars.a = 0.12), (@lens _.c))
br2 = continuation(prob2, PALC(), setproperties(opts, p_max = 1.22);)


# change tolerance for avoiding error computation of the EV
opts_fold = br.contparams
@set! opts_fold.newton_options.eigsolver.σ = 1e-7

brfold = continuation(br2, 3, (@lens _.a),
         setproperties(opts_fold; detect_bifurcation = 1, dsmax = 0.01, max_steps = 70, p_max = 0.6, p_min = -0.6,ds = -0.01, n_inversion = 2, tol_stability = 1e-6);
         verbosity = 1, plot = true,
         detect_codim2_bifurcation = 2,
         update_minaug_every_step = 1,
         bothside = false,
         start_with_eigen = true)

scene = plot(brfold, vars = (:a, :c), branchlabel = "Fold")
plot!(scene, brhopf, vars = (:a, :c), branchlabel = "Hopf")
plot!(scene, brhopf2, vars = (:a, :c), branchlabel = "Hopf")
scene