Neuron model (codim 2, periodic orbits)
Consider the neuron model
\[\left\{\begin{array}{l} \dot{x}_1(t)=-\kappa x_1(t)+\beta \tanh \left(x_1\left(t-\tau_s\right)\right)+a_{12} \tanh \left(x_2\left(t-\tau_2\right)\right) \\ \dot{x_2}(t)=-\kappa x_2(t)+\beta \tanh \left(x_2\left(t-\tau_s\right)\right)+a_{21} \tanh \left(x_1\left(t-\tau_1\right)\right) \end{array}\right.\]
Continuation and codim 1 bifurcations
We first instantiate the model
using Revise, DDEBifurcationKit, Parameters, LinearAlgebra, Plots
using BifurcationKit
const BK = BifurcationKit
function neuronVF(x, xd, p)
@unpack κ, β, a12, a21, τs, τ1, τ2 = p
[
-κ * x[1] + β * tanh(xd[3][1]) + a12 * tanh(xd[2][2]),
-κ * x[2] + β * tanh(xd[3][2]) + a21 * tanh(xd[1][1])
]
end
delaysF(par) = [par.τ1, par.τ2, par.τs]
pars = (κ = 0.5, β = -1, a12 = 1, a21 = 0.5, τ1 = 0.2, τ2 = 0.2, τs = 1.5)
x0 = [0.01, 0.001]
prob = ConstantDDEBifProblem(neuronVF, delaysF, x0, pars, (@lens _.τs))
optn = NewtonPar(verbose = true, eigsolver = DDE_DefaultEig())
opts = ContinuationPar(p_max = 13., p_min = 0., newton_options = optn, ds = -0.01, detect_bifurcation = 3, nev = 5, dsmax = 0.2, n_inversion = 4)
br = continuation(prob, PALC(), opts; verbosity = 1, plot = true, bothside = true, normC = norminf) ┌─ Curve type: EquilibriumCont
├─ Number of points: 60
├─ Type of vectors: Vector{Float64}
├─ Parameter τs starts at 13.0, ends at 0.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, endpoint at τs ≈ +13.00000000, step = 0
- # 2, hopf at τs ≈ +8.90969445 ∈ (+8.90527503, +8.90969445), |δp|=4e-03, [converged], δ = (-2, -2), step = 15, eigenelements in eig[ 16], ind_ev = 4
- # 3, hopf at τs ≈ +1.59721596 ∈ (+1.59714604, +1.59721596), |δp|=7e-05, [converged], δ = (-2, -2), step = 43, eigenelements in eig[ 44], ind_ev = 2
- # 4, endpoint at τs ≈ +0.00000000, step = 59
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 τs ≈ 8.909694447616149.
Frequency ω ≈ 0.8592867628281184
Period of the periodic orbit ≈ 7.312093679297607
Normal form z⋅(iω + a⋅δp + b⋅|z|²):
┌─ a = 0.01253261906342409 - 0.09759405211685644im
└─ b = -0.04708036126526864 + 0.0320112794492103im
Continuation of Hopf points
We follow the Hopf points in the parameter plane $(a_{21},\tau_s)$. We tell the solver to consider br.specialpoint[3] and continue it.
# continuation of the first Hopf point
brhopf = continuation(br, 3, (@lens _.a21),
setproperties(br.contparams, detect_bifurcation = 1, dsmax = 0.04, max_steps = 230, p_max = 15., p_min = -1.,ds = -0.02);
detect_codim2_bifurcation = 2,
# bothside = true,
start_with_eigen = true)
# continuation of the second Hopf point
brhopf2 = continuation(br, 2, (@lens _.a21),
setproperties(br.contparams, detect_bifurcation = 1, dsmax = 0.1, max_steps = 56, p_max = 15., p_min = -1.,ds = -0.01, n_inversion = 4);
detect_codim2_bifurcation = 2,
start_with_eigen = true,
bothside=true)
scene = plot(brhopf, brhopf2, legend = :top)Branch of periodic orbits
We change the continuation parameter and study the bifurcations as function of $a_{21}$.
prob2 = ConstantDDEBifProblem(neuronVF, delaysF, x0, pars, (@lens _.a21))
br2 = BK.continuation(prob2, PALC(), setproperties(opts, ds = 0.1, p_max = 3., n_inversion=8); verbosity = 0, plot = false, normC = norminf) ┌─ Curve type: EquilibriumCont
├─ Number of points: 12
├─ Type of vectors: Vector{Float64}
├─ Parameter a21 starts at 0.5, ends at 3.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 a21 ≈ +0.80713224 ∈ (+0.80711498, +0.80713224), |δp|=2e-05, [converged], δ = ( 2, 2), step = 2, eigenelements in eig[ 3], ind_ev = 2
- # 2, bp at a21 ≈ +2.25000297 ∈ (+2.24999434, +2.25000297), |δp|=9e-06, [ guessL], δ = (-1, 0), step = 8, eigenelements in eig[ 9], ind_ev = 2
- # 3, endpoint at a21 ≈ +3.00000000, step = 11
We then compute the branch of periodic orbits from the Hopf bifurcation points using orthogonal collocation.
# continuation parameters
opts_po_cont = ContinuationPar(dsmax = 0.1, ds= -0.0001, dsmin = 1e-4, p_max = 10., p_min=-0., max_steps = 120, detect_bifurcation = 0, save_sol_every_step=1)
@set! opts_po_cont.newton_options.tol = 1e-8
@set! opts_po_cont.newton_options.verbose = false
# arguments for periodic orbits
args_po = ( record_from_solution = (x, p) -> begin
xtt = BK.get_periodic_orbit(p.prob, x, nothing)
return (max = maximum(xtt[1,:]),
min = minimum(xtt[1,:]),
period = getperiod(p.prob, x, nothing))
end,
plot_solution = (x, p; k...) -> begin
xtt = BK.get_periodic_orbit(p.prob, x, nothing)
plot!(xtt.t, xtt[1,:]; label = "V1", k...)
plot!(xtt.t, xtt[2,:]; label = "V2", k...)
plot!(br2; subplot = 1, putspecialptlegend = false)
end,
normC = norminf)
probpo = PeriodicOrbitOCollProblem(60, 4; N = 2, jacobian = BK.AutoDiffDense())
br_pocoll = @time continuation(
br2, 1, opts_po_cont,
probpo;
verbosity = 0, plot = false,
args_po...,
δp = 0.001,
normC = norminf,
)
scene = plot(br2, br_pocoll)We can plot the periodic orbit as they approach the homoclinic point.
scene = plot(layout = 2)
for ii = 1:10:110
solpo = BK.get_periodic_orbit(br_pocoll.γ.prob.prob, br_pocoll.sol[ii].x, nothing)
plot!(scene, solpo.t ./ solpo.t[end], solpo.u[1,:], label = "", subplot = 1)
end
xlabel!(scene, "t / period", subplot = 1)
plot!(scene, br_pocoll, vars = (:param, :period), subplot = 2, xlims=(2.2,2.4))
scene