(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, Plots
using BifurcationKit
const BK = BifurcationKit
g(z) = (tanh(z − 1) + tanh(1)) * cosh(1)^2
function neuron2VF(x, xd, p)
(; a,b,c,d) = p
[
-x[1] - a * g(b * xd.u[1][1]) + c * g(d * xd.u[2][2]),
-x[2] - a * g(b * xd.u[1][2]) + c * g(d * xd.u[2][1])
]
end
delaysF(par) = [par.τ1, par.τ2]
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, (@optic _.a))
optn = NewtonPar(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.1, n_inversion = 4)
br = continuation(prob, PALC(tangent=Bordered()), opts) ┌─ Curve type: EquilibriumCont
├─ Number of points: 11
├─ Type of vectors: Vector{Float64}
├─ Parameter a starts at 0.25, ends at 1.0
├─ Algo: PALC [Bordered]
└─ Special points:
- # 1, hopf at a ≈ +0.29890362 ∈ (+0.29853073, +0.29890362), |δp|=4e-04, [converged], δ = ( 2, 2), step = 3
- # 2, hopf at a ≈ +0.34537483 ∈ (+0.34523499, +0.34537483), |δp|=1e-04, [converged], δ = ( 2, 2), step = 4
- # 3, hopf at a ≈ +0.36830746 ∈ (+0.36662946, +0.36830746), |δp|=2e-03, [converged], δ = ( 2, 2), step = 5
- # 4, hopf at a ≈ +0.57734590 ∈ (+0.57292648, +0.57734590), |δp|=4e-03, [converged], δ = ( 2, 2), step = 7
- # 5, hopf at a ≈ +1.00000000 ∈ (+0.86018861, +1.00000000), |δp|=1e-01, [converged], δ = ( 2, 2), step = 10
- # 6, endpoint at a ≈ +1.00000000, step = 10
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)SuperCritical - Hopf bifurcation point at a ≈ 0.34537482504934847.
Frequency ω ≈ 0.18264813234949276
Period of the periodic orbit ≈ 34.40049031082817
Normal form z⋅(iω + a⋅δa + b⋅|z|²):
┌─ a = 0.07653825210387889 + 0.08333988810342174im
└─ b = -0.042960448860069445 - 0.04757891550603535im
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, (@optic _.c),
ContinuationPar(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, (@optic _.c),
ContinuationPar(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))
plot!(scene, brhopf2, vars = (:a, :c), xlims = (0, 0.7), ylims = (-0.1, 1))
scene
Continuation of branch points
We follow the branch points in the parameter plane $(a, c)$.
prob2 = ConstantDDEBifProblem(neuron2VF, delaysF, x0, (pars..., a = 0.12), (@optic _.c))
br2 = continuation(prob2, PALC(), ContinuationPar(opts, p_max = 1.22);)
# change tolerance for avoiding error computation of the eigenvalues
opts_bp = br.contparams
@reset opts_bp.newton_options.eigsolver.σ = 1e-7
# index of the first branch point
index_bp = findfirst(x -> x.type == :bp, br2.specialpoint)
brbp = continuation(br2, index_bp, (@optic _.a),
ContinuationPar(opts_bp; 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 = 0,
plot = true,
detect_codim2_bifurcation = 2,
bothside = false,
start_with_eigen = true)
scene = plot(brbp, vars = (:a, :c), branchlabel = "Branch points")
plot!(scene, brhopf, vars = (:a, :c), branchlabel = "Hopf")
plot!(scene, brhopf2, vars = (:a, :c), branchlabel = "Hopf")
scene
Branch of periodic orbits
We compute the branch of periodic orbits from a Hopf bifurcation.
pars = (a = 0.069, b = 2., c = 0.6, d = 1.2, τ1 = 11.6, τ2 = 20.3)
delaysF(par) = [par.τ1, par.τ2]
prob3 = ConstantDDEBifProblem(neuron2VF, delaysF, zeros(2), pars, (@optic _.c))
br3 = continuation(prob3, PALC(), ContinuationPar(opts, p_max = 1.0);) ┌─ Curve type: EquilibriumCont
├─ Number of points: 9
├─ Type of vectors: Vector{Float64}
├─ Parameter c starts at 0.6, ends at 1.0
├─ Algo: PALC [Secant]
└─ Special points:
- # 1, hopf at c ≈ +0.77139744 ∈ (+0.76971944, +0.77139744), |δp|=2e-03, [converged], δ = ( 2, 2), step = 5
- # 2, hopf at c ≈ +0.80914893 ∈ (+0.80911440, +0.80914893), |δp|=3e-05, [converged], δ = ( 2, 2), step = 6
- # 3, hopf at c ≈ +0.92847320 ∈ (+0.92405378, +0.92847320), |δp|=4e-03, [converged], δ = ( 2, 2), step = 7
- # 4, nd at c ≈ +1.00000000 ∈ (+0.92847320, +1.00000000), |δp|=7e-02, [converged], δ = ( 3, 2), step = 8
- # 5, endpoint at c ≈ +1.00000000, step = 8
opts_po_cont = ContinuationPar(dsmax = 0.02, ds = -0.001, max_steps = 185, nev = 10, tol_stability = 1e-4)
probpo = PeriodicOrbitOCollProblem(20, 5; N = 2,
jacobian = BK.AutoDiffDense(),)
br_po = continuation(
br3, 1, ContinuationPar(opts_po_cont; detect_bifurcation = 3),
probpo;
normC = norminf,
eigsolver = BK.FloquetGEV(DDE_DefaultEig(maxit=100, tol = 1e-12, σ = 1e-3)),
)
plot(br_po, vars = (:param, :amplitude))
plot(br3, br_po)