Equating ModelingToolkit parameters and variables?

[AndreaRH]

Hi!

I have a first ODESystem with a parameter gNa and a second ODESystem with a variable g(t). How to create a third system that links the two together so that every time gNa appears in the first system we actually take g(t) ?
Note: I want to preserve the first system which has gNa as a parameter, don’t want to define it as a variable from the start. Just want to extend it when it has dynamics.

Here’s a minimally working example that reproduces the following error.

using ModelingToolkit
const t = Num(Sym{Real}(:t))
const D = Differential(t)

struct Na{T}
    gNa::T
    RNa::T
    mNa::T
    hNa::T
end
m∞(::Na, V) = 1.0 / (1.0 + exp((V + 25.5) / -5.29))
h∞(::Na, V) = 1.0 / (1.0 + exp((V + 48.9) / 5.18))

function channel_dynamics(ch::Na)
    states = @variables mNa(t) hNa(t) INa(t) V(t)
    parameters = @parameters gNa ENa
    eqs = [D(mNa) ~ (m∞(ch, V) - mNa),
        D(hNa) ~ (h∞(ch, V) - hNa),
        INa ~ gNa * mNa^3 * hNa * (ENa - V),
        D(V) ~ INa]
    defaultmap = [mNa => ch.mNa, hNa => ch.hNa, V => -50., gNa => ch.gNa]
    return eqs, states, parameters, defaultmap
end

function regul_indep(ch::Na)
    states = @variables g(t) R(t) Ca(t)
    parameters = @parameters τch Ca_tgt τg
    eqs = [D(R) ~ (Ca_tgt - Ca) / τch,
        D(g) ~ (R - g) / τg,
        D(Ca) ~ (-Ca + 0.5)]
    defaultmap = [g => ch.gNa, R => ch.RNa, Ca =>2.0]
    return eqs, states, parameters, defaultmap
end

ch = Na(70.0, 70.0, 0.0, 0.0)

eqs, states, parameters, defaultmap = channel_dynamics(ch)
@named sys1 = ODESystem(eqs, t, states, parameters; defaults = defaultmap)

eqs, states, parameters, defaultmap = regul_indep(ch)
@named sys2 = ODESystem(eqs, t, states, parameters; defaults = defaultmap)

sys = extend(sys1, sys2)
append!(equations(sys), [sys1.ps[1] ~ sys2.states[1]]) #gNa ~ g(t)
final = alias_elimination(sys)
final = structural_simplify(final)

When I add the extra equation [gNa ~ g(t)] to sys I can see the equations are all OK with alias_elimination. The problem is with structural_simplify, which gives the error
ERROR: LoadError: ExtraEquationsSystemException: The system is unbalanced. There are 7 highest order derivative variables and 8 equations.

I found this related issue but doesn’t seem to be solved: https://github.com/SciML/ModelingToolkit.jl/issues/936

Over-specifying the system with more equations shouldn’t be a problem. What’s the correct way to do this?

[ChrisRackauckas]

Did you open an issue on this?

@YingboMa

[AndreaRH]

No, because I found this issue already existed:

Solving chemical systems of equations with both equilibrium and kinetic species

This is building on a different discourse post I made previously: https://discourse.julialang.org/t/constraints-in-catalyst-jl-for-weak-acids-and-ph-equilibrium/68037/2 And on Github: https://github.com/SciML/Catalyst.jl/issues/398 I am trying to solve a chemical system of equations that combines both equilibrium and kinetic species. Below is the sample code I am using: using OrdinaryDiffEq using Catalyst using DifferentialEquations using Plots using Latexify using ModelingToolkit @paramet…

It’s the same problem I encountered, except the message I got about More equations than variables, here are the potential extra equation(s) did not say which equation was the extra one.

So in my example I decided to re-write the equation I needed to take the g(t) variable.

gstate = sys2.states[1]
gparam = sys1.ps[1]
my_idx = 3 # for the INa eq. 
LHS = equations(sys1)[my_idx].lhs
RHS = equations(sys1)[my_idx].rhs
equations(sys1)[my_idx] = LHS ~ gstate / gparam * RHS 
sys = alias_elimination(extend(sys1, sys2))

[YingboMa]

The system isn’t valid because it’s over-determined. If you want to extend the dynamics later, you shouldn’t set it in the base model.

[AndreaRH]

Yes, I suppose I could have my gNa parameter as a constant variable and then add the dynamics.