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

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.18264813234949256
Period of the periodic orbit ≈ 34.400490310828204
Normal form z⋅(iω + a⋅δp + b⋅|z|²):
┌─ a = 0.07653825213840887 - 0.08333988803527106im
└─ b = 0.07985677803156459 - 0.0018834647622427356im


## 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