🟢 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, Plots
using BifurcationKit
const BK = BifurcationKit

# vector field
function 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 values
par_tm = (α = 1.5, τ = 0.013, J = 3.07, E0 = -2.0, τD = 0.200, U0 = 0.3, τF = 1.5, τS = 0.007)

# initial condition
z0 = [0.238616, 0.982747, 0.367876]

# Bifurcation Problem
prob = BifurcationProblem(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 E0
opts_br = ContinuationPar(p_min = -4.0, p_max = -0.9)

# continuation of equilibria
br = continuation(prob, PALC(), opts_br;
	# we want to compute both sides of the branch of the initial
	# value of E0 = -2
	bothside = 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
 └─ Special points:

- #  1, endpoint at E0 ≈ -4.00000000,                                                                     step =   0
- #  2,       bp at E0 ≈ -1.46309794 ∈ (-1.46309794, -1.46302753), |δp|=7e-05, [converged], δ = ( 1,  0), step =  30
- #  3,     hopf at E0 ≈ -1.85075293 ∈ (-1.85075293, -1.84970242), |δp|=1e-03, [converged], δ = ( 2,  2), step =  40
- #  4,       bp at E0 ≈ -1.86522340 ∈ (-1.86522340, -1.86494948), |δp|=3e-04, [converged], δ = (-1,  0), step =  42
- #  5,     hopf at E0 ≈ -1.14947379 ∈ (-1.15202958, -1.14947379), |δp|=3e-03, [converged], δ = (-2, -2), step =  57
- #  6, endpoint at E0 ≈ -0.90000000,                                                                     step =  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 = ContinuationPar(p_min = -4.0, p_max = -0.9,
	detect_bifurcation = 0)

# continuation of equilibria
br = continuation(prob, PALC(), opts_br;
	# we want to compute both sides of the branch of the initial
	# value of E0 = -2
	bothside = true)

scene = plot(br, plotfold=true)

References