🟢 Neural mass equation

The following model is taken from ^[Cortes]:

\[\left\{\begin{array}{l} \tau \dot{E}=-E+g\left(J u x E+E_{0}\right) \\ \dot{x}=\tau_{D}^{-1}(1-x)-u E x \\ \dot{u}=U_0 E(1-u)-\tau_{F}^{-1}(u-U_0) \end{array}\right.\]

with

\[g(y):=\alpha\log(1+exp(y/\alpha)).\]

We use this model as a mean to introduce the basics of BifurcationKit.jl, namely the continuation of equilibria.

It is easy to encode the ODE as follows

using Revise, Plotsimport BifurcationKit as BKimport BifurcationKit: @optic# vector fieldfunction TMvf(z, p)	(;J, α, E0, τ, τD, τF, U0) = p	E, x, u = z	SS0 = J * u * x * E + E0	SS1 = α * log(1 + exp(SS0 / α))	[		(-E + SS1) / τ,		(1.0 - x) / τD - u * x * E,		(U0 - u) / τF +  U0 * (1.0 - u) * E	]end# parameter valuespar_tm = (α = 1.5, τ = 0.013, J = 3.07, E0 = -2.0, τD = 0.200, U0 = 0.3, τF = 1.5, τS = 0.007)# initial conditionz0 = [0.238616, 0.982747, 0.367876]# Bifurcation Problemprob = BK.ODEBifProblem(TMvf, z0, par_tm, (@optic _.E0);	record_from_solution = (x, p; k...) -> (E = x[1], x = x[2], u = x[3]),)

We first compute the branch of equilibria

# continuation options, we limit the parameter range for E0opts_br = BK.ContinuationPar(p_min = -4.0, p_max = -0.9)# continuation of equilibriabr = BK.continuation(prob, BK.PALC(), opts_br;	# we want to compute both sides of the branch of the initial	# value of E0 = -2	bothside = Val(true))scene = plot(br, legend=:topleft)

With detailed information:

br

 ┌─ Curve type: EquilibriumCont
 ├─ Number of points: 61
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter E0 starts at -4.0, ends at -0.9
 ├─ Algo: PALC [Secant]
 └─ Special points:

┌────┬────────────┬─────────────┬────────────────────────────┬───────────┬─────────────┬──────────┬──────┐
│ Nb │    type    │     E0      │          interval          │ precision │   status    │    δ     │ step │
├────┼────────────┼─────────────┼────────────────────────────┼───────────┼─────────────┼──────────┼──────┤
│  1 │   endpoint │ -4.00000000 │                            │           │             │          │    0 │
│  2 │         bp │ -1.46309794 │ (-1.46309794, -1.46302753) │     7e-05 │   converged │ ( 1,  0) │   30 │
│  3 │       hopf │ -1.85075293 │ (-1.85075293, -1.84970242) │     1e-03 │   converged │ ( 2,  2) │   40 │
│  4 │         bp │ -1.86522340 │ (-1.86522340, -1.86494948) │     3e-04 │   converged │ (-1,  0) │   42 │
│  5 │       hopf │ -1.14947379 │ (-1.15202958, -1.14947379) │     3e-03 │   converged │ (-2, -2) │   57 │
│  6 │   endpoint │ -0.90000000 │                            │           │             │          │   60 │
└────┴────────────┴─────────────┴────────────────────────────┴───────────┴─────────────┴──────────┴──────┘

If you only want to compute the branch without the bifurcations (more information is provided Detection of bifurcation points of Equilibria ), change the continuation options to

opts_br = BK.ContinuationPar(p_min = -4.0, p_max = -0.9,	detect_bifurcation = 0)# continuation of equilibriabr = BK.continuation(prob, BK.PALC(), opts_br;	# we want to compute both sides of the branch of the initial	# value of E0 = -2	bothside = Val(true))scene = plot(br, plotfold=true)

References