Tsodyks-Markram Neural Mass Model

Replicating BifurcationKit.jl TMModel.jl with all 4 PO discretization methods

Overview

This example replicates the Tsodyks-Markram tutorial from BifurcationKit.jl’s official documentation, demonstrating all four periodic orbit discretization methods:

1. Collocation — Orthogonal collocation at Gauss points
2. Trapezoid — Finite-difference trapezoid rule
3. Shooting — Standard multiple shooting
4. Poincaré — Poincaré shooting on a hyperplane

The model describes short-term synaptic plasticity with three variables: - E: Mean firing rate - x: Synaptic depression (available neurotransmitter) - u: Synaptic facilitation (release probability)

Model Definition (Inline YAML with Dynamics.from_string)

from tvbo import Dynamics

model_yaml = """
name: TsodyksMarkram
description: Tsodyks-Markram mean-field model with synaptic plasticity

parameters:
    J: {name: J, value: 3.07, domain: {lo: 0, hi: 10}}
    alpha: {name: alpha, value: 1.4, domain: {lo: 0.1, hi: 5}}
    E0: {name: E0, value: -2.0, domain: {lo: -10, hi: 2}}
    tau: {name: tau, value: 0.013, domain: {lo: 0.001, hi: 0.1}}
    tauD: {name: tauD, value: 0.20, domain: {lo: 0.01, hi: 1.0}}
    tauF: {name: tauF, value: 1.5, domain: {lo: 0.1, hi: 5.0}}
    U0: {name: U0, value: 0.3, domain: {lo: 0, hi: 1}}

derived_variables:
    SS0: {name: SS0, equation: {lhs: SS0, rhs: "J*u*x*E + E0"}}
    SS1: {name: SS1, equation: {lhs: SS1, rhs: "alpha*log(1 + exp(SS0/alpha))"}}

state_variables:
    E:
        name: E
        initial_value: 0.238616
        equation: {lhs: "Derivative(E, t)", rhs: "(-E + SS1)/tau"}
    x:
        name: x
        initial_value: 0.982747
        equation: {lhs: "Derivative(x, t)", rhs: "(1 - x)/tauD - u*x*E"}
    u:
        name: u
        initial_value: 0.367876
        equation: {lhs: "Derivative(u, t)", rhs: "(U0 - u)/tauF + U0*(1 - u)*E"}

number_of_modes: 1
"""

# Load model from inline YAML string
model = Dynamics.from_string(model_yaml)

Method 1: Collocation (Adaptive Mesh)

Orthogonal collocation with automatic mesh refinement:

from tvbo import SimulationExperiment, Continuation
po_collocation = """
name: tm_po_collocation
dynamics: TsodyksMarkram
free_parameters:
    - name: E0
      domain: {lo: "-2", hi: "-1"}
bothside: true
branches:
    - name: po_from_hopf
      source_point: "hopf:2"
      bothside: false
      continuation:
          name: tm_po_collocation_po_from_hopf
          ds: 0.001
          ds_min: 1.0e-4
          ds_max: 0.1
          max_steps: 110
          tol_stability: 1.0e-5
      discretization:
          method: collocation
          mesh_intervals: 50
          degree: 4
          linear_solver:
              method: COPBLS
          parameters:
              - {name: mesh_adaptation, value: true}
"""

cont = Continuation.from_string(po_collocation)
exp = SimulationExperiment(dynamics=model, continuations=[cont])

print(exp.render_code('bifurcationkit.jl'))

using BifurcationKit
using OrdinaryDiffEq






function TsodyksMarkram!(dx, _x, p, t = 0)

    (;J, alpha, E0, tau, tauD, tauF, U0) = p

    E, x, u = _x



    SS0 = E0 + E .* J .* u .* x
    SS1 = alpha .* log(1 + exp(SS0 ./ alpha))

    dx[1] = (SS1 - E) ./ tau
    dx[2] = (1 - x) ./ tauD - E .* u .* x
    dx[3] = (U0 - u) ./ tauF + E .* U0 .* (1 - u)
    dx
end


# Parameter values
p = (J = 3.07, alpha = 1.4, E0 = -2.0, tau = 0.013, tauD = 0.2, tauF = 1.5, U0 = 0.3)

# Override continuation parameter to start within [p_min, p_max]
p = merge(p, (E0 = -2.0,))

# Initial conditions from model defaults
x0 = [
        0.238616, # Initial value for E
        0.982747, # Initial value for x
        0.367876, # Initial value for u
    ]

# Wrapper: BifurcationKit expects f!(du, x, p) (no explicit time argument)
function TsodyksMarkram_vf!(du, x, p)
    TsodyksMarkram!(du, x, p, 0.0)  # pass dummy time
    return du
end

# Find a steady state via time integration (more robust than raw Newton on x0)
function _find_steady_state(f!, x0, p; T=2000.0)
    function ode_f!(du, u, _p, t)
        f!(du, u, p, t)
    end
    prob_ode = ODEProblem{true, SciMLBase.FullSpecialize}(ode_f!, x0, (0.0, T), p)
    sol = solve(prob_ode; abstol=1e-10, reltol=1e-10, save_everystep=false)
    return sol[:, end]
end

x0_eq = _find_steady_state(TsodyksMarkram!, x0, p)


# Record named state variables for each continuation step
record_from_sol = (x, p; k...) -> (E = x[1], x = x[2], u = x[3],)

# Bifurcation Problem
prob = BifurcationProblem(TsodyksMarkram_vf!, x0_eq, p, (@optic _.E0);
    record_from_solution = record_from_sol)

# ContinuationPar
opts_br = ContinuationPar(p_min = -2.0,
    p_max = -1.0)

using Logging
prev_logger = current_logger()
global_logger(SimpleLogger(devnull, Logging.Error))

br = continuation(prob, PALC(), opts_br; normC = norminf, bothside = true)

global_logger(prev_logger)

bifurcation_result = br




# Record PO envelope (max/min per state variable)
args_po = ( record_from_solution = (x, p; k...) -> begin
        xtt = get_periodic_orbit(p.prob, x, p.p)
        return (
                max_E = maximum(xtt[1,:]),
                min_E = minimum(xtt[1,:]),
                max_x = maximum(xtt[2,:]),
                min_x = minimum(xtt[2,:]),
                max_u = maximum(xtt[3,:]),
                min_u = minimum(xtt[3,:]),
                period = getperiod(p.prob, x, p.p))
    end,
    plot_solution = (x, p; k...) -> begin
        xtt = get_periodic_orbit(p.prob, x, p.p)
        arg = (marker = :d, markersize = 1)
        plot!(xtt.t, xtt[1,:]; label = "E", arg..., k...)
    plot!(xtt.t, xtt[2,:]; label = "x", arg..., k...)
    plot!(xtt.t, xtt[3,:]; label = "u", arg..., k...)
        plot!(br; subplot = 1, putspecialptlegend = false)
        end,
    normC = norminf)

opts_po_cont = ContinuationPar(opts_br, ds = 0.001, dsmin = 0.0001, dsmax = 0.1, max_steps = 110, tol_stability = 1e-05)

hopf_indices = Int[]
for (i, sp) in enumerate(br.specialpoint)
    sp.type == :hopf && push!(hopf_indices, i)
end
if !isempty(hopf_indices)
    hopf_indices = [hopf_indices[2]]
end

po_branches = Any[]
for hopf_idx in hopf_indices
    try
        br_po = continuation(
            br, hopf_idx, opts_po_cont,
            PeriodicOrbitOCollProblem(50, 4, meshadapt = true);
            plot = false,
            args_po...,
            linear_algo = COPBLS(),
            verbosity = 0,
        )
        push!(po_branches, br_po)
    catch e
        @warn "PO continuation from Hopf $hopf_idx failed" exception=(e, catch_backtrace())
    end
end

po_results = (hopf_indices = hopf_indices, branches = po_branches)

result = exp.run('bifurcationkit.jl')
result.plot()

Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython

Click to expand: Generated Julia code for collocation method

print(exp.render_code('bifurcationkit.jl'))

using BifurcationKit
using OrdinaryDiffEq






function TsodyksMarkram!(dx, _x, p, t = 0)

    (;J, alpha, E0, tau, tauD, tauF, U0) = p

    E, x, u = _x



    SS0 = E0 + E .* J .* u .* x
    SS1 = alpha .* log(1 + exp(SS0 ./ alpha))

    dx[1] = (SS1 - E) ./ tau
    dx[2] = (1 - x) ./ tauD - E .* u .* x
    dx[3] = (U0 - u) ./ tauF + E .* U0 .* (1 - u)
    dx
end


# Parameter values
p = (J = 3.07, alpha = 1.4, E0 = -2.0, tau = 0.013, tauD = 0.2, tauF = 1.5, U0 = 0.3)

# Override continuation parameter to start within [p_min, p_max]
p = merge(p, (E0 = -2.0,))

# Initial conditions from model defaults
x0 = [
        0.238616, # Initial value for E
        0.982747, # Initial value for x
        0.367876, # Initial value for u
    ]

# Wrapper: BifurcationKit expects f!(du, x, p) (no explicit time argument)
function TsodyksMarkram_vf!(du, x, p)
    TsodyksMarkram!(du, x, p, 0.0)  # pass dummy time
    return du
end

# Find a steady state via time integration (more robust than raw Newton on x0)
function _find_steady_state(f!, x0, p; T=2000.0)
    function ode_f!(du, u, _p, t)
        f!(du, u, p, t)
    end
    prob_ode = ODEProblem{true, SciMLBase.FullSpecialize}(ode_f!, x0, (0.0, T), p)
    sol = solve(prob_ode; abstol=1e-10, reltol=1e-10, save_everystep=false)
    return sol[:, end]
end

x0_eq = _find_steady_state(TsodyksMarkram!, x0, p)


# Record named state variables for each continuation step
record_from_sol = (x, p; k...) -> (E = x[1], x = x[2], u = x[3],)

# Bifurcation Problem
prob = BifurcationProblem(TsodyksMarkram_vf!, x0_eq, p, (@optic _.E0);
    record_from_solution = record_from_sol)

# ContinuationPar
opts_br = ContinuationPar(p_min = -2.0,
    p_max = -1.0)

using Logging
prev_logger = current_logger()
global_logger(SimpleLogger(devnull, Logging.Error))

br = continuation(prob, PALC(), opts_br; normC = norminf, bothside = true)

global_logger(prev_logger)

bifurcation_result = br




# Record PO envelope (max/min per state variable)
args_po = ( record_from_solution = (x, p; k...) -> begin
        xtt = get_periodic_orbit(p.prob, x, p.p)
        return (
                max_E = maximum(xtt[1,:]),
                min_E = minimum(xtt[1,:]),
                max_x = maximum(xtt[2,:]),
                min_x = minimum(xtt[2,:]),
                max_u = maximum(xtt[3,:]),
                min_u = minimum(xtt[3,:]),
                period = getperiod(p.prob, x, p.p))
    end,
    plot_solution = (x, p; k...) -> begin
        xtt = get_periodic_orbit(p.prob, x, p.p)
        arg = (marker = :d, markersize = 1)
        plot!(xtt.t, xtt[1,:]; label = "E", arg..., k...)
    plot!(xtt.t, xtt[2,:]; label = "x", arg..., k...)
    plot!(xtt.t, xtt[3,:]; label = "u", arg..., k...)
        plot!(br; subplot = 1, putspecialptlegend = false)
        end,
    normC = norminf)

opts_po_cont = ContinuationPar(opts_br, ds = 0.001, dsmin = 0.0001, dsmax = 0.1, max_steps = 110, tol_stability = 1e-05)

hopf_indices = Int[]
for (i, sp) in enumerate(br.specialpoint)
    sp.type == :hopf && push!(hopf_indices, i)
end
if !isempty(hopf_indices)
    hopf_indices = [hopf_indices[2]]
end

po_branches = Any[]
for hopf_idx in hopf_indices
    try
        br_po = continuation(
            br, hopf_idx, opts_po_cont,
            PeriodicOrbitOCollProblem(50, 4, meshadapt = true);
            plot = false,
            args_po...,
            linear_algo = COPBLS(),
            verbosity = 0,
        )
        push!(po_branches, br_po)
    catch e
        @warn "PO continuation from Hopf $hopf_idx failed" exception=(e, catch_backtrace())
    end
end

po_results = (hopf_indices = hopf_indices, branches = po_branches)

Method 2: Trapezoid (Finite Differences)

po_trapezoid = """
name: tm_po_trapezoid
dynamics: TsodyksMarkram
free_parameters:
    - name: E0
      domain: {lo: "-2", hi: "-1"}
bothside: false
branches:
    - name: po_from_hopf
      source_point: "hopf:2"
      bothside: false
      delta_p: 0.001
      continuation:
          name: tm_po_trapezoid_po_from_hopf
          ds: 0.004
          ds_min: 1.0e-4
          ds_max: 0.1
          max_steps: 80
          newton_tol: 1.0e-7
          tol_stability: 1.0e-7
      discretization:
          method: trapezoid
          mesh_intervals: 250
"""

cont_trap = Continuation.from_string(po_trapezoid)
exp_trap = SimulationExperiment(dynamics=model, continuations=[cont_trap])
result_trap = exp_trap.run('bifurcationkit.jl')
result_trap.plot()

Method 3: Standard Shooting

Multiple ODE integrations:

po_shooting = """
name: tm_po_shooting
dynamics: TsodyksMarkram
free_parameters:
    - name: E0
      domain: {lo: "-2", hi: "-1"}
bothside: false
branches:
    - name: po_from_hopf
      source_point: "hopf:2"
      bothside: false
      continuation:
          name: tm_po_shooting_po_from_hopf
          ds: -0.0001
          ds_min: 1.0e-4
          ds_max: 0.1
          max_steps: 110
          tol_stability: 1.0e-4
      discretization:
          method: shooting
          n_sections: 15
          ode_solver:
              method: Rodas5
              abs_tol: 1e-11
              rel_tol: 1e-9
          linear_solver:
              method: MatrixBLS
          parameters:
              - {name: parallel, value: true}
      parameters:
          - {name: max_norm_bound, value: 10}
"""

cont_shoot = Continuation.from_string(po_shooting)
exp_shoot = SimulationExperiment(dynamics=model, continuations=[cont_shoot])
print(exp_shoot.render_code('bifurcationkit.jl'))
result_shoot = exp_shoot.run('bifurcationkit.jl')
result_shoot.plot()

using BifurcationKit
using OrdinaryDiffEq






function TsodyksMarkram!(dx, _x, p, t = 0)

    (;J, alpha, E0, tau, tauD, tauF, U0) = p

    E, x, u = _x



    SS0 = E0 + E .* J .* u .* x
    SS1 = alpha .* log(1 + exp(SS0 ./ alpha))

    dx[1] = (SS1 - E) ./ tau
    dx[2] = (1 - x) ./ tauD - E .* u .* x
    dx[3] = (U0 - u) ./ tauF + E .* U0 .* (1 - u)
    dx
end


# Parameter values
p = (J = 3.07, alpha = 1.4, E0 = -2.0, tau = 0.013, tauD = 0.2, tauF = 1.5, U0 = 0.3)

# Override continuation parameter to start within [p_min, p_max]
p = merge(p, (E0 = -2.0,))

# Initial conditions from model defaults
x0 = [
        0.238616, # Initial value for E
        0.982747, # Initial value for x
        0.367876, # Initial value for u
    ]

# Wrapper: BifurcationKit expects f!(du, x, p) (no explicit time argument)
function TsodyksMarkram_vf!(du, x, p)
    TsodyksMarkram!(du, x, p, 0.0)  # pass dummy time
    return du
end

# Find a steady state via time integration (more robust than raw Newton on x0)
function _find_steady_state(f!, x0, p; T=2000.0)
    function ode_f!(du, u, _p, t)
        f!(du, u, p, t)
    end
    prob_ode = ODEProblem{true, SciMLBase.FullSpecialize}(ode_f!, x0, (0.0, T), p)
    sol = solve(prob_ode; abstol=1e-10, reltol=1e-10, save_everystep=false)
    return sol[:, end]
end

x0_eq = _find_steady_state(TsodyksMarkram!, x0, p)


# Record named state variables for each continuation step
record_from_sol = (x, p; k...) -> (E = x[1], x = x[2], u = x[3],)

# Bifurcation Problem
prob = BifurcationProblem(TsodyksMarkram_vf!, x0_eq, p, (@optic _.E0);
    record_from_solution = record_from_sol)

# ContinuationPar
opts_br = ContinuationPar(p_min = -2.0,
    p_max = -1.0)

using Logging
prev_logger = current_logger()
global_logger(SimpleLogger(devnull, Logging.Error))

br = continuation(prob, PALC(), opts_br; normC = norminf)

global_logger(prev_logger)

bifurcation_result = br




# Record PO envelope (max/min per state variable)
args_po = ( record_from_solution = (x, p; k...) -> begin
        xtt = get_periodic_orbit(p.prob, x, p.p)
        return (
                max_E = maximum(xtt[1,:]),
                min_E = minimum(xtt[1,:]),
                max_x = maximum(xtt[2,:]),
                min_x = minimum(xtt[2,:]),
                max_u = maximum(xtt[3,:]),
                min_u = minimum(xtt[3,:]),
                period = getperiod(p.prob, x, p.p))
    end,
    plot_solution = (x, p; k...) -> begin
        xtt = get_periodic_orbit(p.prob, x, p.p)
        arg = (marker = :d, markersize = 1)
        plot!(xtt.t, xtt[1,:]; label = "E", arg..., k...)
    plot!(xtt.t, xtt[2,:]; label = "x", arg..., k...)
    plot!(xtt.t, xtt[3,:]; label = "u", arg..., k...)
        plot!(br; subplot = 1, putspecialptlegend = false)
        end,
    normC = norminf)

opts_po_cont = ContinuationPar(opts_br, ds = -0.0001, dsmin = 0.0001, dsmax = 0.1, max_steps = 110, tol_stability = 0.0001)

hopf_indices = Int[]
for (i, sp) in enumerate(br.specialpoint)
    sp.type == :hopf && push!(hopf_indices, i)
end
if !isempty(hopf_indices)
    hopf_indices = [hopf_indices[2]]
end

po_branches = Any[]
for hopf_idx in hopf_indices
    try
        prob_ode = ODEProblem(TsodyksMarkram!, copy(x0), (0.0, 1000.0), p; abstol = 1e-11, reltol = 1e-09)
        br_po = continuation(
            br, hopf_idx, opts_po_cont,
            ShootingProblem(15, prob_ode, OrdinaryDiffEq.Rodas5(), parallel = true);
            plot = false,
            args_po...,
            linear_algo = MatrixBLS(),
            verbosity = 0,
            callback_newton = BifurcationKit.cbMaxNorm(10.0),
        )
        push!(po_branches, br_po)
    catch e
        @warn "PO continuation from Hopf $hopf_idx failed" exception=(e, catch_backtrace())
    end
end

po_results = (hopf_indices = hopf_indices, branches = po_branches)

┌ Warning: The precision on the Floquet multipliers is 2.7325737521308788e-8.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 1.6608328073070823e-8.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 6.0082349790657e-8.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 9.880189533591978e-8.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 3.0721707767302713e-7.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 1.3692887806327293e-8.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 2.0132061834608382e-7.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 2.053749993407368e-8.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 4.619286624211238e-8.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 6.629710205322572e-8.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 2.7429833284198744e-7.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 8.543418545033322e-7.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 9.352746726560644e-7.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 1.5273626913157446e-8.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 1.0217241889194383e-8.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 2.4033044820767986e-8.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 8.650176405149776e-8.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 2.462513068246181e-8.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 5.292521616683653e-8.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 2.981830153426688e-7.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 6.7424765105167595e-6.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 7.435102071965049e-5.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.0018416937626528865.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.0006300493108742717.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 1.5178653135453317e-5.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.00029068875963047153.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.0011994743538063198.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.03654867173517135.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 5.173064841587385.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 3.3851362095468835.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 3.508758292012897.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 3.290534010435416.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 3.263381794732664.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 9.994243285299158.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80

┌ Warning: The precision on the Floquet multipliers is 17.314820800964128.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80

Method 4: Poincaré Shooting

Reduced-dimension shooting:

po_poincare = """
name: tm_po_poincare
dynamics: TsodyksMarkram
free_parameters:
    - name: E0
      domain: {lo: "-2", hi: "-1"}
bothside: false
branches:
    - name: po_from_hopf
      source_point: "hopf:2"
      bothside: false
      continuation:
          name: tm_po_poincare_po_from_hopf
          ds: 0.0001
          ds_max: 0.02
          max_steps: 50
          newton_tol: 1.0e-9
          newton_max_iterations: 15
          tol_stability: 1.0e-6
          detect_bifurcation: 2
      discretization:
          method: poincare
          n_sections: 5
          ode_solver:
              method: Rodas5
              abs_tol: 1e-11
              rel_tol: 1e-9
          linear_solver:
              method: MatrixBLS
          parameters:
              - {name: parallel, value: true}
      parameters:
          - {name: max_norm_bound, value: 1.0}
"""

cont_poinc = Continuation.from_string(po_poincare)
exp_poinc = SimulationExperiment(dynamics=model, continuations=[cont_poinc])
print(exp_poinc.render_code('bifurcationkit.jl'))
result_poinc = exp_poinc.run('bifurcationkit.jl')
result_poinc.plot()

using BifurcationKit
using OrdinaryDiffEq






function TsodyksMarkram!(dx, _x, p, t = 0)

    (;J, alpha, E0, tau, tauD, tauF, U0) = p

    E, x, u = _x



    SS0 = E0 + E .* J .* u .* x
    SS1 = alpha .* log(1 + exp(SS0 ./ alpha))

    dx[1] = (SS1 - E) ./ tau
    dx[2] = (1 - x) ./ tauD - E .* u .* x
    dx[3] = (U0 - u) ./ tauF + E .* U0 .* (1 - u)
    dx
end


# Parameter values
p = (J = 3.07, alpha = 1.4, E0 = -2.0, tau = 0.013, tauD = 0.2, tauF = 1.5, U0 = 0.3)

# Override continuation parameter to start within [p_min, p_max]
p = merge(p, (E0 = -2.0,))

# Initial conditions from model defaults
x0 = [
        0.238616, # Initial value for E
        0.982747, # Initial value for x
        0.367876, # Initial value for u
    ]

# Wrapper: BifurcationKit expects f!(du, x, p) (no explicit time argument)
function TsodyksMarkram_vf!(du, x, p)
    TsodyksMarkram!(du, x, p, 0.0)  # pass dummy time
    return du
end

# Find a steady state via time integration (more robust than raw Newton on x0)
function _find_steady_state(f!, x0, p; T=2000.0)
    function ode_f!(du, u, _p, t)
        f!(du, u, p, t)
    end
    prob_ode = ODEProblem{true, SciMLBase.FullSpecialize}(ode_f!, x0, (0.0, T), p)
    sol = solve(prob_ode; abstol=1e-10, reltol=1e-10, save_everystep=false)
    return sol[:, end]
end

x0_eq = _find_steady_state(TsodyksMarkram!, x0, p)


# Record named state variables for each continuation step
record_from_sol = (x, p; k...) -> (E = x[1], x = x[2], u = x[3],)

# Bifurcation Problem
prob = BifurcationProblem(TsodyksMarkram_vf!, x0_eq, p, (@optic _.E0);
    record_from_solution = record_from_sol)

# ContinuationPar
opts_br = ContinuationPar(p_min = -2.0,
    p_max = -1.0)

using Logging
prev_logger = current_logger()
global_logger(SimpleLogger(devnull, Logging.Error))

br = continuation(prob, PALC(), opts_br; normC = norminf)

global_logger(prev_logger)

bifurcation_result = br




# Record PO envelope (max/min per state variable)
args_po = ( record_from_solution = (x, p; k...) -> begin
        xtt = get_periodic_orbit(p.prob, x, p.p)
        return (
                max_E = maximum(xtt[1,:]),
                min_E = minimum(xtt[1,:]),
                max_x = maximum(xtt[2,:]),
                min_x = minimum(xtt[2,:]),
                max_u = maximum(xtt[3,:]),
                min_u = minimum(xtt[3,:]),
                period = getperiod(p.prob, x, p.p))
    end,
    plot_solution = (x, p; k...) -> begin
        xtt = get_periodic_orbit(p.prob, x, p.p)
        arg = (marker = :d, markersize = 1)
        plot!(xtt.t, xtt[1,:]; label = "E", arg..., k...)
    plot!(xtt.t, xtt[2,:]; label = "x", arg..., k...)
    plot!(xtt.t, xtt[3,:]; label = "u", arg..., k...)
        plot!(br; subplot = 1, putspecialptlegend = false)
        end,
    normC = norminf)

opts_po_cont = ContinuationPar(opts_br, ds = 0.0001, dsmax = 0.02, max_steps = 50, tol_stability = 1e-06, detect_bifurcation = 2, newton_options = NewtonPar(tol = 1e-09, max_iterations = 15))

hopf_indices = Int[]
for (i, sp) in enumerate(br.specialpoint)
    sp.type == :hopf && push!(hopf_indices, i)
end
if !isempty(hopf_indices)
    hopf_indices = [hopf_indices[2]]
end

po_branches = Any[]
for hopf_idx in hopf_indices
    try
        prob_ode = ODEProblem(TsodyksMarkram!, copy(x0), (0.0, 1000.0), p; abstol = 1e-11, reltol = 1e-09)
        br_po = continuation(
            br, hopf_idx, opts_po_cont,
            PoincareShootingProblem(5, prob_ode, OrdinaryDiffEq.Rodas5(), parallel = true);
            plot = false,
            args_po...,
            linear_algo = MatrixBLS(),
            verbosity = 0,
            callback_newton = BifurcationKit.cbMaxNorm(1.0),
        )
        push!(po_branches, br_po)
    catch e
        @warn "PO continuation from Hopf $hopf_idx failed" exception=(e, catch_backtrace())
    end
end

po_results = (hopf_indices = hopf_indices, branches = po_branches)

┌ Warning: The precision on the Floquet multipliers is 0.0017909760569191637.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.0018828209805575177.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.0020191811303651376.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.002224173828630299.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.002537209906612921.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.0030244330012906004.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.00379941752373453.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.005003853176015959.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.006903079664666693.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.00992039286214593.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.014688694739140175.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.02202096026559467.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.032548096533339954.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.04522347053831872.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.05160803258031889.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.04213113410655399.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: The precision on the Floquet multipliers is 0.038771208490223696.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: Verbosity toggle: max_iters 
│  Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems).
└ @ SciMLBase ~/.julia/packages/SciMLBase/eTBDr/src/integrator_interface.jl:679
┌ Warning: Verbosity toggle: max_iters 
│  Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems).
└ @ SciMLBase ~/.julia/packages/SciMLBase/eTBDr/src/integrator_interface.jl:679
┌ Warning: The precision on the Floquet multipliers is 0.03747986179125429.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80
┌ Warning: Verbosity toggle: max_iters 
│  Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems).
└ @ SciMLBase ~/.julia/packages/SciMLBase/eTBDr/src/integrator_interface.jl:679
┌ Warning: The precision on the Floquet multipliers is 0.03699846553120567.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetQaD`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/periodicorbit/Floquet.jl:80

┌ Error: Failure to converge with given tolerance = 1.0e-9.
│ Step = 18
│ You can decrease the tolerance or pass a different norm using the argument `normC`.
│ We reached the smallest value [dsmin] valid for ds, namely 0.0001.
│ Stopping continuation at continuation step 18.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/Y3DPO/src/continuation/Contbase.jl:73