LabelledArrays Error

[donaldlee3]

Getting the following error with the last line of this example:
ArgumentError: invalid index: :a of type Symbol

using DifferentialEquations
using LabelledArrays

function lorenz!(du,u,p,t)
 du.a = 10.0*(u.b-u.a)
 du.b = u.a*(28.0-u.c) - u.b
 du.c = u.a*u.b - (8/3)*u.c
end
u0 = LVector(a=1.0, b=0.0, c=0.0)
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan)
sol2 = solve(prob)

sol2(0:10, idxs=1)
sol2(0:10, idxs=:a)

Removed all the equations from my actual model and still getting errors:

using DifferentialEquations
using Sundials
using LabelledArrays
using Plots

################################################################################

## For BloodFlow

lag             =   1.           # sec       # introduces a lag in CBF response
tstim           =   60.          # sec

builtin1(x) = 1.5*(-1/5*exp(-x/5)+1/2*exp(-x/2))
builtin2(x,lag,taurise) = exp((x-lag)/taurise)
builtin3(x) = (1.1+1.5*(exp(-x/5)-exp(-x/2)))

const t0        =   lag
const t1        =   tstim
const F0        =   0.012       # sec^-1
const taurise   =   0.1*1/builtin1(0.)
const lastF0    =   F0*builtin3(t1-t0)

## For ExcitatoryConductance

const Ne        =   240.
const finf      =   2.5         # Hz
const f0        =   23.         # Hz
const tau       =   2.5         # sec
const gbar      =   7.8e-06     # mS*cm^-2*sec

################################################################################

function InhibitoryConductance(t)
    return 0.
end

function model(du,u,p,t)
    ############################################################################
    ## BloodFlow

    if p.t_BF_interval == 1        # if t < t0
        #taurise     =   0.1*1/builtin1(0)
        BF          =   F0*(1+0.1*builtin2(t,t0,taurise))
    elseif p.t_BF_interval == 2    # if t0 ≤ t < t1
        BF          =   F0*builtin3(t-t0)
    else                        # if t ≥ t1
        #lastF0      =   F0*builtin3(t1-t0)
        BF          =   F0+(lastF0-F0)*exp(-(t-tstim)/5.)
    end

    ############################################################################
    ## ExcitatoryConductance

    if p.t_EC_interval == 2
        ge          =   Ne*gbar*(finf+(f0-finf)*exp(-t/tau))
    else
        ge          =   0
    end

    ############################################################################

    du .= 0
end

u0  =   LVector(
        # Variables from metabolic model
        Na_n        =   7.9717478717769090e+000,  # 1 ?????????? Need to look into this value
        Na_g        =   1.5105962705116731e+001,  # 2
        GLC_n       =   1.1952401106798076e+000,  # 3
        GLC_g       =   1.1916831724306007e+000,  # 4
        GAP_n       =   4.6031315201704066e-003,  # 5
        GAP_g       =   4.5921537563994577e-003,  # 6
        PEP_n       =   1.6409369231409576e-002,  # 7
        PEP_g       =   1.4203935093932973e-002,  # 8
        PYR_n       =   1.7389794176222959e-001,  # 9
        PYR_g       =   1.6297929416942439e-001,  # 10
        LAC_n       =   5.9812378340393302e-001,  # 11
        LAC_g       =   6.0013834405304178e-001,  # 12
        NADH_n      =   6.2436588774363953e-003,  # 13
        NADH_g      =   1.0386865243459381e-001,  # 14
        NADH_n_mito =   1.2345447794054716e-001,  # 15  # 32
        NADH_g_mito =   1.2483363148236443e-001,  # 16  # 33
        ATP_n       =   2.1973762129517347e+000,  # 17  # 15
        ATP_g       =   2.1951002997965290e+000,  # 18  # 16
        PCr_n       =   4.9460223423181864e+000,  # 19  # 17
        PCr_g       =   4.9242059083825644e+000,  # 20  # 18
        O2_n        =   2.7584951432999300e-002,  # 21  # 19
        O2_g        =   2.7594056537295686e-002,  # 22  # 20
        O2_c        =   6.9808721131857396e+000,  # 23  # 21
        GLC_c       =   4.5003866284951739e+000,  # 24  # 22
        LAC_c       =   5.4885078074249904e-001,  # 25  # 23
        V_v         =   2.0999999999999998e-002,  # 26  # 24
        dHb         =   5.7503371246197574e-002,  # 27  # 25
        GLC_e       =   2.4796908575815539e+000,  # 28  # 26
        LAC_e       =   5.9913109309463841e-001,  # 29  # 27

        # Variables from neuronal model
        V_m         =   -7.3592796321379041e+001,  # 30  # 28
        h           =   9.9371941956423548e-001,  # 31  # 29
        n           =   1.8509072079738374e-002,  # 32  # 30
        Ca_n        =   5.1006943148441310e-005,  # 33  # 31
        )

Ve                      =   0.2
Vcap                    =   0.0055
Vg                      =   0.25
Vn                      =   0.45
MVF                     =   0.07
FourCIN                 =   1.

p   =   LVector(### parameters.jl

        ## Compartment volumes (per unit tissue volume)
        Ve              =   0.2,
        Vcap            =   0.0055,
        Vg              =   0.25,
        Vn              =   0.45,
        MVF             =   0.07,           # mitochondrial volume fraction

        ren             =   Ve/Vn,
        reg             =   Ve/Vg,
        rce             =   Vcap/Ve,
        rcn             =   Vcap/Vn,
        rcg             =   Vcap/Vg,

        ## Surface to Volume ratios and sodium conductances

        SmVn            =   2.5e04,         # cm^-1
        SmVg            =   2.5e04,         # cm^-1

        ## Physical constants and parameters

        R               =   8.314510,       # J mol^-1 K^-1
        F               =   9.64853e04,     # C mol^-1
        RTF             =   26.73,          # mV (corresponds to body temperature 37°C)
        Nae             =   150.,            # mmol/L
        Vm              =   -70.,            # mV

        ## Na/K-ATPase


        KmPump          =   0.5,            # mmol/L


        ## GLC exchange constants

        KtGLCen         =   8.,              # mmol/L
        KtGLCeg         =   8.,              # mmol/L
        KtGLCcg         =   8.,              # mmol/L
        KtGLCce         =   8.,              # mmol/L

        ## Hexokinase-phosphofructokinase system

        KIATP           =   1.,              # mmol/L
        nH              =   4.,
        Kg              =   0.05,           # mmol/L

        ## Phosphoglycerate kinase

        N               =   0.212,          # mmol/L

        KtLACne         =   0.74,           # mmol/L
        KtLACge         =   3.5,            # mmol/L
        KtLACgc         =   1.0,            # mmol/L        # see cremer79.pdf
        KtLACec         =   1.0,            # mmol/L        # and leegsma-vogt01.pdf

        ## Mitochondrial respiration

        KmMito          =   0.04,           # mmol/L
        KO2Mito	        =   0.001,          # mmol/L


        C               =   10.,             # mmol/L

        ## Oxygen exchange constants

        KO2             =   0.0361,         # mmol/L
        HbOP            =   8.6,            # mmol/L
        nh              =   2.73,

        ## Blood flow contributions

        O2a             =   8.35,           # mmol/L
        GLCa            =	4.75,           # mmol/L

        ## Venous flow

        tauv            =   35.,             # secs
        alphav          =   0.5,

        ## Ratio of excitatory conductance

        glia            =   50.,             # mV

        ## Neuronal parameters

        Cm              =   1e-03,          # mF/cm2
        gL              =   0.02,           # mS/cm2
        gNa             =   40.,             # mS/cm2
        gK              =   18.,             # mS/cm2
        gCa             =   0.02,           # mS/cm2
        gmAHP           =   6.5,            # mS/cm2
        KD              =   30e-03,         # mmol/L
        tauCa           =   150e-03,        # sec
        Ca0             =   0.5e-04,        # mmol/L
        EK              =   -80.,            # mV
        ECa             =   120.,            # mV
        Ee              =   0.,              # mV
        Ei              =   -80.,            # mV
        phih            =   4.,
        phin            =   4.,

        ### NewCoeffsStab.jl

        ## This is auto man!

        kCKnps          =	4.33e-002,
        kCKgps	        =	1.35e-003,
        KmNADn          =	4.09e-001,
        KmNADg          =	4.03e+001,
        KmNADHn	        =	4.44e-002,
        KmNADHg	        =	2.69e-002,
        kLDHnps	        =	7.23e+001,
        kLDHgps	        =	1.59,
        MnCyto	        =	4.9e-008,
        MgCyto	        =	2.5e-004,
        MnMito	        =	3.93e+005,
        MgMito	        =	1.06e+004,
        KmADPn	        =	3.41e-003,
        KmADPg	        =	4.83e-004,
        TMaxLACgc       =	2.59e-003,    # Don't see in optimized parameters table, but
                                          # see similar in constrained parameters (DEL)

        ## Constrained parameters

        TMaxGLCen       =   0.041,
        TMaxGLCce       =   0.239,
        TMaxGLCeg       =   0.147,
        TMaxGLCcg       =   0.0016,
        kHKPFKn         =   0.0504,
        kHKPFKg         =   0.185,
        kPGKn           =   3.97,
        kPGKg           =   401.7,
        kPKn            =   36.7,
        kPKg            =   135.2,
        kCKnms          =   0.00028,
        kCKgms          =   0.00001,
        PScapVn         =   1.66,
        PScapVg         =   0.87,
        VMaxMitooutn    =   0.164,
        VMaxMitooutg    =   0.064,
        TnNADH          =   10330,
        TgNADH          =   150.,
        VMaxMitoinn     =   0.1303,
        VMaxMitoing     =   5.7,
        vATPasesn       =   0.1695,
        vATPasesg       =   0.1404,
        kLDHnms         =   0.72,
        kLDHgms         =   0.071,
        TMaxLACne       =   FourCIN*24.3,
        TMaxLACge       =   106.1,
        TMaxLACec       =   0.25,
        LACa            =   0.506,

        ##

        kPumpn          =   2.2e-06,
        gNan            =   0.0136,
        gNag            =   0.0061,
        gKpas           =   0.2035,
        kPumpg          =   4.5e-07,
        vPumpg0         =   0.0687,

        #****************
        #include("ADP.jl")
        qAK             =  0.92,
        A               =  2.212,                                           # mmol/L

        ################for BloodFlow
        t_BF_interval   =   1,

        ################for ExcitatoryConductance
        t_EC_interval   =   1,

        ################for MeanFieldJstim
        x_interval      =   1,

        ################added from u0
        Na0             =   u0.Na_n,
        Vv0             =   u0.V_v
        )

prob = ODEProblem(model, u0, (-60., 800.), p)
sol = solve(prob, CVODE_BDF(), reltol=1e-12, abstol=1e-12, callback=cbs)

sol(-60:0, idxs=1)
sol(-60:0, idxs=:Na_n)

sol(-60:0, idxs=1) gives the error:

MethodError: Cannot `convert` an object of type LArray{Float64,0,Array{Float64,0},(:Na_n, :Na_g, :GLC_n, :GLC_g, :GAP_n, :GAP_g, :PEP_n, :PEP_g, :PYR_n, :PYR_g, :LAC_n, :LAC_g, :NADH_n, :NADH_g, :NADH_n_mito, :NADH_g_mito, :ATP_n, :ATP_g, :PCr_n, :PCr_g, :O2_n, :O2_g, :O2_c, :GLC_c, :LAC_c, :V_v, :dHb, :GLC_e, :LAC_e, :V_m, :h, :n, :Ca_n)} to an object of type Float64

sol(-60:0, idxs=:Na_n) gives the error:

MethodError: Cannot 'convert' an object of type Float64 to an object of type Array{Float64,1}

[donaldlee3]

callback=cbs needs to be removed from the previous code.

I'm reattaching my model code below with equations removed. The lines at the bottom below the hashtags are giving me problems.

using DifferentialEquations
using Sundials
using LabelledArrays
using Plots

################################################################################

## For BloodFlow

lag             =   1.           # sec       # introduces a lag in CBF response
tstim           =   60.          # sec

builtin1(x) = 1.5*(-1/5*exp(-x/5)+1/2*exp(-x/2))
builtin2(x,lag,taurise) = exp((x-lag)/taurise)
builtin3(x) = (1.1+1.5*(exp(-x/5)-exp(-x/2)))

const t0        =   lag
const t1        =   tstim
const F0        =   0.012       # sec^-1
const taurise   =   0.1*1/builtin1(0.)
const lastF0    =   F0*builtin3(t1-t0)

## For ExcitatoryConductance

const Ne        =   240.
const finf      =   2.5         # Hz
const f0        =   23.         # Hz
const tau       =   2.5         # sec
const gbar      =   7.8e-06     # mS*cm^-2*sec

################################################################################

function InhibitoryConductance(t)
    return 0.
end

function model(du,u,p,t)

    ############################################################################
    ## BloodFlow

    if p.t_BF_interval == 1        # if t < t0
        #taurise     =   0.1*1/builtin1(0)
        BF          =   F0*(1+0.1*builtin2(t,t0,taurise))
    elseif p.t_BF_interval == 2    # if t0 ≤ t < t1
        BF          =   F0*builtin3(t-t0)
    else                        # if t ≥ t1
        #lastF0      =   F0*builtin3(t1-t0)
        BF          =   F0+(lastF0-F0)*exp(-(t-tstim)/5.)
    end

    ############################################################################
    ## ExcitatoryConductance

    if p.t_EC_interval == 2
        ge          =   Ne*gbar*(finf+(f0-finf)*exp(-t/tau))
    else
        ge          =   0
    end

    ############################################################################

    du .= 0
end

u0  =   LVector(
        # Variables from metabolic model
        Na_n        =   7.9717478717769090e+000,  # 1 ?????????? Need to look into this value
        Na_g        =   1.5105962705116731e+001,  # 2
        GLC_n       =   1.1952401106798076e+000,  # 3
        GLC_g       =   1.1916831724306007e+000,  # 4
        GAP_n       =   4.6031315201704066e-003,  # 5
        GAP_g       =   4.5921537563994577e-003,  # 6
        PEP_n       =   1.6409369231409576e-002,  # 7
        PEP_g       =   1.4203935093932973e-002,  # 8
        PYR_n       =   1.7389794176222959e-001,  # 9
        PYR_g       =   1.6297929416942439e-001,  # 10
        LAC_n       =   5.9812378340393302e-001,  # 11
        LAC_g       =   6.0013834405304178e-001,  # 12
        NADH_n      =   6.2436588774363953e-003,  # 13
        NADH_g      =   1.0386865243459381e-001,  # 14
        NADH_n_mito =   1.2345447794054716e-001,  # 15  # 32
        NADH_g_mito =   1.2483363148236443e-001,  # 16  # 33
        ATP_n       =   2.1973762129517347e+000,  # 17  # 15
        ATP_g       =   2.1951002997965290e+000,  # 18  # 16
        PCr_n       =   4.9460223423181864e+000,  # 19  # 17
        PCr_g       =   4.9242059083825644e+000,  # 20  # 18
        O2_n        =   2.7584951432999300e-002,  # 21  # 19
        O2_g        =   2.7594056537295686e-002,  # 22  # 20
        O2_c        =   6.9808721131857396e+000,  # 23  # 21
        GLC_c       =   4.5003866284951739e+000,  # 24  # 22
        LAC_c       =   5.4885078074249904e-001,  # 25  # 23
        V_v         =   2.0999999999999998e-002,  # 26  # 24
        dHb         =   5.7503371246197574e-002,  # 27  # 25
        GLC_e       =   2.4796908575815539e+000,  # 28  # 26
        LAC_e       =   5.9913109309463841e-001,  # 29  # 27

        # Variables from neuronal model
        V_m         =   -7.3592796321379041e+001,  # 30  # 28
        h           =   9.9371941956423548e-001,  # 31  # 29
        n           =   1.8509072079738374e-002,  # 32  # 30
        Ca_n        =   5.1006943148441310e-005,  # 33  # 31
        )

Ve                      =   0.2
Vcap                    =   0.0055
Vg                      =   0.25
Vn                      =   0.45
MVF                     =   0.07
FourCIN                 =   1.

p   =   LVector(### parameters.jl

        ## Compartment volumes (per unit tissue volume)
        Ve              =   0.2,
        Vcap            =   0.0055,
        Vg              =   0.25,
        Vn              =   0.45,
        MVF             =   0.07,           # mitochondrial volume fraction

        ren             =   Ve/Vn,
        reg             =   Ve/Vg,
        rce             =   Vcap/Ve,
        rcn             =   Vcap/Vn,
        rcg             =   Vcap/Vg,

        ## Surface to Volume ratios and sodium conductances

        SmVn            =   2.5e04,         # cm^-1
        SmVg            =   2.5e04,         # cm^-1

        ## Physical constants and parameters

        R               =   8.314510,       # J mol^-1 K^-1
        F               =   9.64853e04,     # C mol^-1
        RTF             =   26.73,          # mV (corresponds to body temperature 37°C)
        Nae             =   150.,            # mmol/L
        Vm              =   -70.,            # mV

        ## Na/K-ATPase


        KmPump          =   0.5,            # mmol/L


        ## GLC exchange constants

        KtGLCen         =   8.,              # mmol/L
        KtGLCeg         =   8.,              # mmol/L
        KtGLCcg         =   8.,              # mmol/L
        KtGLCce         =   8.,              # mmol/L

        ## Hexokinase-phosphofructokinase system

        KIATP           =   1.,              # mmol/L
        nH              =   4.,
        Kg              =   0.05,           # mmol/L

        ## Phosphoglycerate kinase

        N               =   0.212,          # mmol/L

        KtLACne         =   0.74,           # mmol/L
        KtLACge         =   3.5,            # mmol/L
        KtLACgc         =   1.0,            # mmol/L        # see cremer79.pdf
        KtLACec         =   1.0,            # mmol/L        # and leegsma-vogt01.pdf

        ## Mitochondrial respiration

        KmMito          =   0.04,           # mmol/L
        KO2Mito	        =   0.001,          # mmol/L


        C               =   10.,             # mmol/L

        ## Oxygen exchange constants

        KO2             =   0.0361,         # mmol/L
        HbOP            =   8.6,            # mmol/L
        nh              =   2.73,

        ## Blood flow contributions

        O2a             =   8.35,           # mmol/L
        GLCa            =	4.75,           # mmol/L

        ## Venous flow

        tauv            =   35.,             # secs
        alphav          =   0.5,

        ## Ratio of excitatory conductance

        glia            =   50.,             # mV

        ## Neuronal parameters

        Cm              =   1e-03,          # mF/cm2
        gL              =   0.02,           # mS/cm2
        gNa             =   40.,             # mS/cm2
        gK              =   18.,             # mS/cm2
        gCa             =   0.02,           # mS/cm2
        gmAHP           =   6.5,            # mS/cm2
        KD              =   30e-03,         # mmol/L
        tauCa           =   150e-03,        # sec
        Ca0             =   0.5e-04,        # mmol/L
        EK              =   -80.,            # mV
        ECa             =   120.,            # mV
        Ee              =   0.,              # mV
        Ei              =   -80.,            # mV
        phih            =   4.,
        phin            =   4.,

        ### NewCoeffsStab.jl

        ## This is auto man!

        kCKnps          =	4.33e-002,
        kCKgps	        =	1.35e-003,
        KmNADn          =	4.09e-001,
        KmNADg          =	4.03e+001,
        KmNADHn	        =	4.44e-002,
        KmNADHg	        =	2.69e-002,
        kLDHnps	        =	7.23e+001,
        kLDHgps	        =	1.59,
        MnCyto	        =	4.9e-008,
        MgCyto	        =	2.5e-004,
        MnMito	        =	3.93e+005,
        MgMito	        =	1.06e+004,
        KmADPn	        =	3.41e-003,
        KmADPg	        =	4.83e-004,
        TMaxLACgc       =	2.59e-003,    # Don't see in optimized parameters table, but
                                          # see similar in constrained parameters (DEL)

        ## Constrained parameters

        TMaxGLCen       =   0.041,
        TMaxGLCce       =   0.239,
        TMaxGLCeg       =   0.147,
        TMaxGLCcg       =   0.0016,
        kHKPFKn         =   0.0504,
        kHKPFKg         =   0.185,
        kPGKn           =   3.97,
        kPGKg           =   401.7,
        kPKn            =   36.7,
        kPKg            =   135.2,
        kCKnms          =   0.00028,
        kCKgms          =   0.00001,
        PScapVn         =   1.66,
        PScapVg         =   0.87,
        VMaxMitooutn    =   0.164,
        VMaxMitooutg    =   0.064,
        TnNADH          =   10330,
        TgNADH          =   150.,
        VMaxMitoinn     =   0.1303,
        VMaxMitoing     =   5.7,
        vATPasesn       =   0.1695,
        vATPasesg       =   0.1404,
        kLDHnms         =   0.72,
        kLDHgms         =   0.071,
        TMaxLACne       =   FourCIN*24.3,
        TMaxLACge       =   106.1,
        TMaxLACec       =   0.25,
        LACa            =   0.506,

        ##

        kPumpn          =   2.2e-06,
        gNan            =   0.0136,
        gNag            =   0.0061,
        gKpas           =   0.2035,
        kPumpg          =   4.5e-07,
        vPumpg0         =   0.0687,

        #****************
        #include("ADP.jl")
        qAK             =  0.92,
        A               =  2.212,                                           # mmol/L

        ################for BloodFlow
        t_BF_interval   =   1,

        ################for ExcitatoryConductance
        t_EC_interval   =   1,

        ################for MeanFieldJstim
        x_interval      =   1,

        ################added from u0
        Na0             =   u0.Na_n,
        Vv0             =   u0.V_v
        )

prob = ODEProblem(model, u0, (-60., 800.), p)
sol = solve(prob, CVODE_BDF(), reltol=1e-12, abstol=1e-12)

################################################################################
plot(sol, idxs=1)
plot(sol, idxs=:Na_n)

sol(-60:0, idxs=1)
sol(-60:0, idxs=:Na_n)

O2_n(t) = sol(t, idxs=findfirst(==(:O2_n), keys(u0)))
O2_n.(-60:-59)
mean(O2_n.(-60:-59))

O2_n(t) = sol(t, idxs=:O2_n)
O2_n.(-60:-59)

Adding on to the problems from the earlier comment:

1. Attempts to plot a single ID are resulting in all the variables being plotted.
2. This works:
   O2_n(t) = sol(t, idxs=findfirst(==(:O2_n), keys(u0)))
   But, this doesn't work:
   O2_n(t) = sol(t, idxs=:O2_n)

ERROR: MethodError: no method matching size(::Symbol)
Closest candidates are:
  size(::CSV.ReversedBuf) at C:\Users\julia\.julia\packages\CSV\MKemC\src\utils.jl:403
  size(::BitArray{1}) at bitarray.jl:104
  size(::BitArray{1}, ::Integer) at bitarray.jl:107

when I try O2_n.(-60:0).
From O2_n.(-60:-59) with the former, I get this:

2-element Array{Any,1}:
 0.0275849514329993
  0-dimensional LArray{Float64,0,Array{Float64,0},(:Na_n, :Na_g, :GLC_n, :GLC_g, :GAP_n, :GAP_g, :PEP_n, :PEP_g, :PYR_n, :PYR_g, :LAC_n, :LAC_g, :NADH_n, :NADH_g, :NADH_n_mito, :NADH_g_mito, :ATP_n, :ATP_g, :PCr_n, :PCr_g, :O2_n, :O2_g, :O2_c, :GLC_c, :LAC_c, :V_v, :dHb, :GLC_e, :LAC_e, :V_m, :h, :n, :Ca_n)}:
:Na_n => 0.0275848602195084

mean(O2_n.(-60:-59)) gives me the following error:

MethodError: no method matching +(::Float64, ::LArray{Float64,0,Array{Float64,0},(:Na_n, :Na_g, :GLC_n, :GLC_g, :GAP_n, :GAP_g, :PEP_n, :PEP_g, :PYR_n, :PYR_g, :LAC_n, :LAC_g, :NADH_n, :NADH_g, :NADH_n_mito, :NADH_g_mito, :ATP_n, :ATP_g, :PCr_n, :PCr_g, :O2_n, :O2_g, :O2_c, :GLC_c, :LAC_c, :V_v, :dHb, :GLC_e, :LAC_e, :V_m, :h, :n, :Ca_n)})
For element-wise addition, use broadcasting with dot syntax: scalar .+ array
Closest candidates are:
  +(::Any, ::Any, !Matched::Any, !Matched::Any...) at operators.jl:538
  +(!Matched::ChainRulesCore.DoesNotExist, ::Any) at C:\Users\julia\.julia\packages\ChainRulesCore\UayCG\src\differential_arithmetic.jl:23
  +(::Float64, !Matched::Float64) at float.jl:401

[ChrisRackauckas]

Transferring to Sundials. It's a Sundials issue since what's happening is that the solution is not actually calculated using LabelledArrays in Sundials

[ChrisRackauckas]

Only standard arrays are supported by the Sundials solver.