Hello,
I ran into the following ERROR: type NameState has no field symbolify when trying to create an ODEProblem. The stack trace seemed to me impossible to decipher, any help would be appreciated!
Example below:
using ModelingToolkit, DifferentialEquations
const t = Num(Sym{Real}(:t))
const D = Differential(t)
params = @parameters Cₘ, τCa, Ca∞, calc_multiplier, τav, Na₊g, Na₊ENa, CaS₊g, CaT₊g, Ka₊g, Ka₊EK, KCa₊g, KCa₊EK, Kdr₊g, Kdr₊EK, H₊g, H₊EH, Input₊g, leak₊g, leak₊ELeak
vars = @variables V(t), Ca(t), Avg(t), ECa(t), Na₊mNa(t), Na₊hNa(t), Na₊INa(t), CaS₊mCaS(t), CaS₊hCaS(t), CaS₊ICa(t), CaT₊mCaT(t), CaT₊hCaT(t), CaT₊ICa(t), Ka₊mKa(t), Ka₊hKa(t), Ka₊IK(t), KCa₊mKCa(t), KCa₊IK(t), Kdr₊mKdr(t), Kdr₊IK(t), H₊mH(t), H₊IH(t), Input₊ILeak(t), leak₊ILeak(t)
eqs = [D(V)~(Na₊INa+Ka₊IK+Input₊ILeak+KCa₊IK+Kdr₊IK+leak₊ILeak+CaT₊ICa+H₊IH+CaS₊ICa)/Cₘ,
D(Ca)~(Ca∞-Ca+(calc_multiplier*(CaT₊ICa+CaS₊ICa))/Cₘ)/τCa,
D(Avg)~(Ca-Avg)/τav,
ECa ~ 12.20008405*log(max(3000.0 / Ca, 0.001)),
D(Na₊mNa)~(1.0/(1.0+exp(-0.1890359168241966(25.5 + V)))-Na₊mNa)/(1.32+-1.26/(1+exp(-0.04(120.0 + V)))),
D(Na₊hNa)~((-Na₊hNa+1.0/(1.0+exp(0.19305019305019305(48.9 + V))))*(1.0+exp(-0.1(62.9 + V))))/(0.67(1.5+1.0/(1.0+exp(0.2777777777777778(34.9 + V))))),
D(CaS₊mCaS)~(1.0/(1.0+exp(-0.1234567901234568(33.0 + V)))-CaS₊mCaS)/(1.4+7.0/(exp(0.1(27.0 + V))+exp(-0.07692307692307693(70.0 + V)))),
D(CaS₊hCaS)~(-CaS₊hCaS+1.0/(1.0+exp(0.16129032258064516(60.0 + V))))/(60.0+150.0/(exp(0.1111111111111111(55.0 + V))+exp(-0.0625(65.0 + V)))),
D(CaT₊mCaT)~(-CaT₊mCaT+1.0/(1.0+exp(-0.1388888888888889(27.1 + V))))/(21.7+-21.3/(1.0+exp(-0.04878048780487805(68.1 + V)))),
D(CaT₊hCaT)~(-CaT₊hCaT+1.0/(1.0+exp(0.18181818181818182(32.1 + V))))/(105.0+-89.8/(1.0+exp(-0.0591715976331361(55.0 + V)))),
D(Ka₊mKa)~(-Ka₊mKa+1.0/(1.0+exp(-0.1149425287356322(27.2 + V))))/(11.6+-10.4/(1.0+exp(-0.06578947368421052(32.9 + V)))),
D(Ka₊hKa)~(-Ka₊hKa+1.0/(1.0+exp(0.2040816326530612(56.9 + V))))/(38.6+-29.2/(1.0+exp(-0.03773584905660377(38.9 + V)))),
D(KCa₊mKCa)~(Ca/((1.0+exp(-0.07936507936507936(28.3 + V)))*(3.0+Ca))-KCa₊mKCa)/(90.3+-75.1/(1.0+exp(-0.04405286343612335(46.0 + V)))),
D(Kdr₊mKdr)~(1.0/(1.0+exp(-0.0847457627118644(12.3 + V)))-Kdr₊mKdr)/(7.2+-6.4/(1.0+exp(-0.052083333333333336(28.3 + V)))),
D(H₊mH)~(1.0/(1.0+exp(0.16666666666666666(70.0 + V)))-H₊mH)/(272.0+1499.0/(1.0+exp(-0.1145475372279496(42.2 + V)))),
Na₊INa ~ Na₊g * (Na₊ENa - V) * Na₊hNa * (Na₊mNa^3),
CaS₊ICa ~ CaS₊g * (ECa - V) * CaS₊hCaS * (CaS₊mCaS^3),
CaT₊ICa ~ CaT₊g * (ECa - V) * CaT₊hCaT * (CaT₊mCaT^3),
Ka₊IK ~ Ka₊g * (Ka₊EK - V) * Ka₊hKa * (Ka₊mKa^3),
KCa₊IK ~ KCa₊g * (KCa₊EK - V) * (KCa₊mKCa^4),
Kdr₊IK ~ Kdr₊g * (Kdr₊EK - V) * (Kdr₊mKdr^4),
H₊IH ~ H₊g * (H₊EH - V) * H₊mH,
Input₊ILeak ~ Input₊g, #think of input current as a tunable parameter, just like other conductances
leak₊ILeak ~ (leak₊ELeak - V) * leak₊g];
ics = Dict(
Kdr₊g => 50.0,
Na₊ENa => 50.0,
Na₊g => 1000.0,
leak₊g => 0.1,
CaS₊hCaS => 0.5,
leak₊ELeak => -50.0,
Ca∞ => 0.05,
KCa₊g => 186.0,
Kdr₊mKdr => 0.0172529,
CaT₊hCaT => 0.993774,
Input₊ILeak => 0.0,
V => -60.0,
KCa₊EK => -80.0,
H₊g => 0.2,
Ka₊g => 160.0,
Ka₊hKa => 0.653091,
CaS₊g => 12.0,
CaT₊g => 6.0,
τCa => 20.0,
CaT₊mCaT => 0.0102574,
Avg => 1.91964,
Ca => 0.05,
KCa₊mKCa => 0.00122546,
Input₊g => 0.0,
τav => 3000.0,
Na₊mNa => 0.00146895,
H₊EH => -20.0,
Ka₊mKa => 0.0225301,
Kdr₊EK => -80.0,
Cₘ => 10.0,
Ka₊EK => -80.0,
H₊mH => 0.158869,
Na₊hNa => 0.894999,
calc_multiplier => 0.939488,
CaS₊mCaS => 0.0344452,
ECa => 85.0,
Na₊INa => 0.0,
CaS₊ICa => 0.0,
CaT₊ICa => 0.0,
Ka₊IK => 0.0,
KCa₊IK => 0.0,
Kdr₊IK => 0.0,
H₊IH => 0.0,
leak₊ILeak => 0.0
);
sys_basic = structural_simplify(ODESystem(eqs, t, vars, params; defaults=ics, name=:neuron));
prob_basic = ODEProblem{true}(sys_basic, [], (0.0,3000.0), [])
sol_basic = OrdinaryDiffEq.solve(prob_basic, Tsit5(), dt=0.01; save_idxs=[1,2,3])
"""
Recursive gradient of average calcium.
Alternative dynamics need
- gradient of v wrt conductances
∇_g v = [u₁, u₂, ..., uₙ]'
where
u̇ = - (g̃/c) u + Φ(x(t))'
and
Φ(x(t)) = [E_Na * m_Na^3(t) * h_Na(t), E_Cas * m_CaS^3(t) * h_CaS(t), ...]
- gradient of <Ca> wrt conductances
∇_g <Ca> = [z₁, z₂, ..., zₙ]'
where ż = - (1/τₐ)z + w
and ẇ = - (1/τcₐ)w + [∇_g Φ_Ca(u(t),x(t),v(t)) * g_Ca + Φ_Ca(x(t))']
"""
function parametrizedby_g()
sum([ics[Na₊g] * ics[Na₊mNa]^3 * ics[Na₊hNa],
ics[CaS₊g] * ics[CaS₊mCaS]^3 * ics[CaS₊hCaS],
ics[CaT₊g] * ics[CaT₊mCaT]^3 * ics[CaT₊hCaT],
ics[Ka₊g] * ics[Ka₊mKa]^3 * ics[Ka₊hKa] ,
ics[KCa₊g] * ics[KCa₊mKCa]^4,
ics[Kdr₊g] * ics[Kdr₊mKdr]^4,
ics[H₊g] * ics[H₊mH]])
#returns a constant, need to choose where gates are evaluated. Atm iniitial conditions are at -60
end
function Φᵥ()
[Na₊ENa * Na₊mNa^3 * Na₊hNa,
ECa * CaS₊mCaS^3 * CaS₊hCaS, #ECa is function of t
ECa * CaT₊mCaT^3 * CaT₊hCaT,
Ka₊EK * Ka₊mKa^3 * Ka₊hKa,
KCa₊EK * KCa₊mKCa^4,
Kdr₊EK * Kdr₊mKdr^4,
H₊EH * H₊mH]
end
function sigmoid(V, vhalf, slope)
return 1.0 / (1.0 + exp((V + vhalf) / slope))
end
function dsigmoid_dV(V, vhalf, slope)
(-exp((V + vhalf) / slope)) / (slope * (sigmoid(V, vhalf, slope))^2)
end
function dΦ_Ca_dv()
[3 * CaS₊mCaS^2 * CaS₊dmCaS_dV * CaS₊hCaS * (ECa - V) + CaS₊mCaS^3 * CaS₊dhCaS_dV * (ECa - V) - CaS₊mCaS^3 * CaS₊hCaS,
3 * CaT₊mCaT^2 * CaT₊dmCaT_dV * CaT₊hCaT * (ECa - V) + CaT₊mCaT^3 * CaT₊dhCaT_dV * (ECa - V) - CaT₊mCaT^3 * CaT₊hCaT]
end
function ∇_gΦ_Ca()
[dΦ_Ca_dv()[1] .* u dΦ_Ca_dv()[2] .* u]
#return a 7x2 matrix
end
function Φ_Ca()
[0.0,
CaS₊mCaS^3 * CaS₊hCaS * (ECa - V),
CaT₊mCaT^3 * CaT₊hCaT * (ECa - V),
0.0,
0.0,
0.0,
0.0]
end
extra_params = @parameters g̃
extra_vars = @variables (w(t))[1:7], (u(t))[1:7], (z(t))[1:7], CaS₊dmCaS_dV(t), CaS₊dhCaS_dV(t), CaT₊dmCaT_dV(t), CaT₊dhCaT_dV(t)
extra_eqs = [
D(CaS₊dmCaS_dV) ~ (-CaS₊dmCaS_dV + dsigmoid_dV(V, 33.0, -8.1)) / (1.4 + 7.0 / (exp(0.1(27.0 + V)) + exp(-0.07692307692307693(70.0 + V)))),
D(CaS₊dhCaS_dV) ~ (-CaS₊dhCaS_dV + dsigmoid_dV(V, 60.0, 6.2)) / (60.0 + 150.0 / (exp(0.1111111111111111(55.0 + V)) + exp(-0.0625(65.0 + V)))),
D(CaT₊dmCaT_dV) ~ (-CaT₊dmCaT_dV + dsigmoid_dV(V, 27.1, -7.2)) / (21.7 + -21.3 / (1.0 + exp(-0.04878048780487805(68.1 + V)))),
D(CaT₊dhCaT_dV) ~ (-CaT₊dhCaT_dV + dsigmoid_dV(V, 32.1, 5.5)) / (105.0 + -89.8 / (1.0 + exp(-0.0591715976331361(55.0 + V)))),
[D(u[i]) ~ -(g̃ / Cₘ) * (u[i] + Φᵥ()[i]) for i in 1:7]...,
[D(w[i]) ~ -(1 / τCa) * (w[i] + (∇_gΦ_Ca()*[CaS₊g, CaT₊g] + Φ_Ca())[i]) for i in 1:7]...,
[D(z[i]) ~ -(1 / τav) * (z[i] + w[i]) for i in 1:7]...,
]
extra_ics = merge(Dict(
g̃ => parametrizedby_g(),
CaS₊dmCaS_dV => 0.0,
CaS₊dhCaS_dV => 0.0,
CaT₊dmCaT_dV => 0.0,
CaT₊dhCaT_dV => 0.0),
Dict([u...] .=> [0.0 for i = 1:7]),
Dict([w...] .=> [0.0 for i=1:7]),
Dict([z...] .=> [0.0 for i = 1:7]))
Eqs = vcat(eqs, extra_eqs)
Params = vcat(params, extra_params)
Vars = vcat(vars, extra_vars)
Ics = merge(ics, extra_ics)
sys_recursive = ODESystem(Eqs, t, Vars, Params; defaults=Ics, name=:neuron);
sys_recursive = structural_simplify(sys_recursive)
prob_recursive = ODEProblem{true}(sys_recursive, [], (0.0, 5000.0), []) # error here!
Version Info:
Julia 1.10.3
ModelingToolkit v9.14.0
DifferentialEquations v7.13.0
There’s an issue which points to the same error so I wondered if it’s a bug that will be likely resolved soon or if there’s a workaround atm?
Thanks!