MWE to illustrate the requirement for fix_unassigned_nonconstant_par_to_init (without it I get: ERROR: ExtraVariablesSystemException: The system is unbalanced. There are 4 highest order derivative variables and 3 equations.)

mwe_fix_unassigned_nonconstant_par_to_init.jl

using SBMLToolkit
using SBML
using Catalyst
using Downloads

case = "00170"

function read_case(case)
    base_url = "https://raw.githubusercontent.com/sbmlteam/sbml-test-suite/master/cases/semantic/$case/$case"
    sbml_url = base_url*"-sbml-l3v2.xml"
    sbml = String(take!(Downloads.download(sbml_url, IOBuffer())))
end

# Read case
sbml = read_case(case)

# Read SBML
ml = readSBMLFromString(sbml, doc -> begin
            set_level_and_version(3, 2)(doc)
            convert_simplify_math(doc)
        end)

## Works
function mwe_ReactionSystem1(model::SBML.Model; kwargs...)  # Todo: requires unique parameters (i.e. SBML must have been imported with localParameter promotion in libSBML)
    # length(model.events) > 0 ? error("Model contains events. Please import with `ODESystem(model)`") : nothing  @Anand: how to suppress this when called from ODESystem
    rxs = SBMLToolkit.mtk_reactions(model)
    u0map = SBMLToolkit.get_u0map(model)
    parammap = SBMLToolkit.get_paramap(model)
    defs = ModelingToolkit._merge(Dict(u0map), Dict(parammap))

    algrules, obsrules, raterules = SBMLToolkit.get_rules(model)
    # constant_eqs = fix_constants_at_init(model)  # Todo PL: take this out
    unassigned_pars = SBMLToolkit.fix_unassigned_nonconstant_par_to_init(model)  # @Sam: if I do not do that, there are too few equations. Do you have better suggestions?
    for o in obsrules
        defs[o.lhs] = substitute(o.rhs, defs)
    end
    constraints_sys = ODESystem(vcat(algrules, raterules, obsrules, unassigned_pars),
                                Catalyst.DEFAULT_IV; name = gensym(:CONSTRAINTS))

    # Hacky way to keep zero_ode species with contant=boundaryCondition=false in system
    rs = ReactionSystem(rxs, Catalyst.DEFAULT_IV, first.(u0map), first.(parammap); defaults = defs, name = gensym(:SBML),
        constraints = constraints_sys, kwargs...)
    odessys = convert(ODESystem, rs; include_zero_odes=true, combinatoric_ratelaws=false)
    u0map, fixed_to_init = SBMLToolkit.get_u0map(model, equations(odessys))  # Todo PL: Use input=true instead
    constraints_sys = ODESystem(vcat(algrules, raterules, obsrules, unassigned_pars, fixed_to_init),
                                Catalyst.DEFAULT_IV; name = gensym(:CONSTRAINTS))
    ReactionSystem(rxs, Catalyst.DEFAULT_IV, first.(u0map), first.(parammap); defaults = defs, name = gensym(:SBML),
        constraints = constraints_sys, kwargs...)
end

rs1 = mwe_ReactionSystem1(ml)  # 4 states, 4 equations
sys1 = convert(ODESystem, rs1; include_zero_odes = false, combinatoric_ratelaws=false)
ssys = structural_simplify(sys1)

## Fails
function mwe_ReactionSystem2(model::SBML.Model; kwargs...)  # Todo: requires unique parameters (i.e. SBML must have been imported with localParameter promotion in libSBML)
    # length(model.events) > 0 ? error("Model contains events. Please import with `ODESystem(model)`") : nothing  @Anand: how to suppress this when called from ODESystem
    rxs = SBMLToolkit.mtk_reactions(model)
    u0map = SBMLToolkit.get_u0map(model)
    parammap = SBMLToolkit.get_paramap(model)
    defs = ModelingToolkit._merge(Dict(u0map), Dict(parammap))

    algrules, obsrules, raterules = SBMLToolkit.get_rules(model)
    # constant_eqs = fix_constants_at_init(model)  # Todo PL: take this out
    # unassigned_pars = SBMLToolkit.fix_unassigned_nonconstant_par_to_init(model)  # @Sam: if I do not do that, there are too few equations. Do you have better suggestions?
    for o in obsrules
        defs[o.lhs] = substitute(o.rhs, defs)
    end
    constraints_sys = ODESystem(vcat(algrules, raterules, obsrules),
                                Catalyst.DEFAULT_IV; name = gensym(:CONSTRAINTS))

    # Hacky way to keep zero_ode species with contant=boundaryCondition=false in system
    rs = ReactionSystem(rxs, Catalyst.DEFAULT_IV, first.(u0map), first.(parammap); defaults = defs, name = gensym(:SBML),
        constraints = constraints_sys, kwargs...)
    odessys = convert(ODESystem, rs; include_zero_odes=true, combinatoric_ratelaws=false)
    u0map, fixed_to_init = SBMLToolkit.get_u0map(model, equations(odessys))  # Todo PL: Use input=true instead
    constraints_sys = ODESystem(vcat(algrules, raterules, obsrules, fixed_to_init),
                                Catalyst.DEFAULT_IV; name = gensym(:CONSTRAINTS))
    ReactionSystem(rxs, Catalyst.DEFAULT_IV, first.(u0map), first.(parammap); defaults = defs, name = gensym(:SBML),
        constraints = constraints_sys, kwargs...)
end

rs2 = mwe_ReactionSystem2(ml)  # 4 states, 3 equations
sys2 = convert(ODESystem, rs2; include_zero_odes = false, combinatoric_ratelaws=false)
ssys = structural_simplify(sys2)  # ERROR: ExtraVariablesSystemException: The system is unbalanced. There are 4 highest order derivative variables and 3 equations.