Use equations directly (#59)

Author: sosiristseng

.vscode/settings.json

@@ -461,9 +461,9 @@ version = "1.15.1"
 
 [[deps.DifferentiationInterface]]
 deps = ["ADTypes", "LinearAlgebra"]
-git-tree-sha1 = "f620da805b82bec64ab4d5f881c7592c82dbc08a"
+git-tree-sha1 = "54d7b8c74408048aea9c055ac8573b2b5c5ec11f"
 uuid = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
-version = "0.7.3"
+version = "0.7.4"
 
     [deps.DifferentiationInterface.extensions]
     DifferentiationInterfaceChainRulesCoreExt = "ChainRulesCore"
@@ -996,9 +996,9 @@ version = "3.1.1+0"
 
 [[deps.JuliaFormatter]]
 deps = ["CommonMark", "Glob", "JuliaSyntax", "PrecompileTools", "TOML"]
-git-tree-sha1 = "7b3a1a49526f1549e24032520afd812790c37e5b"
+git-tree-sha1 = "f512fefd5fdc7dd1ca05778f08f91e9e4c9fdc37"
 uuid = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
-version = "2.1.5"
+version = "2.1.6"
 
 [[deps.JuliaSyntax]]
 git-tree-sha1 = "937da4713526b96ac9a178e2035019d3b78ead4a"

Manifest.toml

@@ -88,33 +88,18 @@ function build_model(; name=:ecmesys, bcl=0second, istim=-80μAcm⁻², tstart=0
     end
 
     # Subsystems
-    lccsys = get_lcc_sys(; ca_ss, ca_o, k_i, k_o, vm)
-    ryrsys = get_ryr_sys(; ca_jsr, ca_ss)
-    cksys = get_ck_sys(; atp_i, adp_i)
-    forcesys = get_force_sys(; atp_i, adp_i, ca_i)
-    jcasys = get_jca_sys(; atp_i, adp_i, ca_i, ca_nsr, ca_jsr, ca_ss, ca_o, na_i, na_o, vm)
-    iksys = get_ik_sys(; na_i, na_o, k_i, k_o, mg_i, vm, atp_i, adp_i)
-    inasys = get_ina_sys(; na_i, na_o, ca_i, ca_o, vm)
-    inaksys = get_inak_sys(; atp_i, adp_i, vm, na_i, na_o, k_o)
-    tcassys = get_tca_sys(; atp_m, adp_m, nad_m, nadh_m, h_m, ca_m, pi_m, mg_m, use_mg)
-    @unpack suc, fum, oaa, vSL, vIDH, vKGDH, vMDH = tcassys
-    etcsys = get_etc_sys(; h_i, h_m, DOX, MT_PROT, O2, nad_m, nadh_m, suc, fum, oaa)
-    c5sys = get_c5_sys(; dpsi, h_i, h_m, atp_i, adp_i, atp_m, adp_m, pi_m, MT_PROT, mg_i, mg_m, use_mg)
-    @unpack vROS, vNADHC1, vHres, vO2 = etcsys
-    rossys = get_ros_sys(; dpsi, sox_m, nadph_i, V_MITO_V_MYO)
-    mitocasys = get_mitoca_sys(; na_i, ca_m, ca_i, dpsi)
-
-    @unpack vTrROS, vIMAC, vSOD_i = rossys
-    @unpack INa, INsNa, ICaB, INaB = inasys
-    @unpack IPMCA, Jup, INaCa, Jtr, Jxfer = jcasys
-    @unpack IK, IK1, IKp, IKatp = iksys
-    @unpack INaK = inaksys
-    @unpack ICaK, ICaL = lccsys
-    @unpack vCK_mito = cksys
-    @unpack vANT, vC5, vHu, vHleak = c5sys
-    @unpack Jrel = ryrsys
-    @unpack vAm, Jtrpn = forcesys
-    @unpack vUni, vNaCa = mitocasys
+    @unpack eqs_cicr16, ICaL, ICaK, Jrel, Jtr, Jxfer = get_cicr16_eqs(; vm, ca_ss, ca_i, ca_jsr, ca_nsr, k_i, ca_o, k_o)
+    @unpack eqs_ck, vCK_mito, vCK_cyto, vTR_crp = get_ck_eqs(; atp_i, adp_i)
+    @unpack eqs_force, vAm, Jtrpn = get_force_eqs(; atp_i, adp_i, ca_i)
+    @unpack eqs_jca, IPMCA, Jup, INaCa = get_jca_eqs(; atp_i, adp_i, ca_i, ca_nsr, ca_o, na_i, na_o, vm)
+    @unpack eqs_ik, IK1, IK, IKp, IKatp = get_ik_eqs(; na_i, na_o, k_i, k_o, vm, atp_i, adp_i, mg_i)
+    @unpack eqs_ina, INa, INsNa, INaB, ICaB = get_ina_eqs(; na_i, na_o, ca_i, ca_o, vm)
+    @unpack eqs_nak, INaK = get_inak_eqs(; atp_i, adp_i, vm, na_i, na_o, k_o)
+    @unpack eqs_tca, vIDH, vKGDH, vMDH, vSL, suc, fum, oaa = get_tca_eqs(; atp_m, adp_m, nad_m, nadh_m, ca_m, h_m, pi_m, mg_m, use_mg)
+    @unpack eqs_etc, vHres, vROS, vSDH, vNADHC1, vO2 = get_etc_eqs(; nad_m, nadh_m, dpsi, sox_m, suc, fum, oaa, DOX, MT_PROT, O2, h_i, h_m)
+    @unpack eqs_c5, vANT, vC5, vHu, vHleak = get_c5_eqs(; dpsi, h_i, h_m, atp_i, adp_i, atp_m, adp_m, pi_m, MT_PROT, mg_i, mg_m)
+    @unpack eqs_ros, vTrROS, vIMAC, vCAT, vGPX_i, vGR_i, vSOD_i = get_ros_eqs(; dpsi, sox_m, nadph_i, V_MITO_V_MYO)
+    @unpack eqs_mitoca, vUni, vNaCa = get_mitoca_eqs(; na_i, ca_m, ca_i, dpsi)
 
     eqs = [
         ΣA_i ~ atp_i + adp_i,
@@ -136,21 +121,19 @@ function build_model(; name=:ecmesys, bcl=0second, istim=-80μAcm⁻², tstart=0
         D(O2) ~ V_MT / (V_MYO + V_MT) * (-vROS - vO2) + kdiffO2 * (O2_o - O2),
     ]
 
+    eqs = [eqs; eqs_cicr16; eqs_ck; eqs_force; eqs_jca; eqs_ik; eqs_ina; eqs_nak; eqs_tca; eqs_etc; eqs_c5; eqs_ros; eqs_mitoca]
+
     if bcl > 0 && duty > 0 && istim != 0
         ts0 = collect(tstart:bcl:tend)
         ts1 = ts0 .+ duty
 
-        sstart = ModelingToolkit.SymbolicDiscreteCallback( ts0 => [iStim ~ istim], discrete_parameters = iStim, iv = t)
-        send = ModelingToolkit.SymbolicDiscreteCallback( ts1 => [iStim ~ 0], discrete_parameters = iStim, iv = t)
+        sstart = ModelingToolkit.SymbolicDiscreteCallback( ts0 => [iStim ~ istim]; discrete_parameters = [iStim], iv = t)
+        send = ModelingToolkit.SymbolicDiscreteCallback( ts1 => [iStim ~ 0]; discrete_parameters = [iStim], iv = t)
         sys = System(eqs, t; name, discrete_events = [sstart, send])
     else
         sys = System(eqs, t; name)
     end
 
-    for s in (lccsys, ryrsys, cksys, forcesys, jcasys, iksys, inasys, inaksys, tcassys, etcsys, c5sys, rossys, mitocasys)
-        sys = extend(s, sys; name)
-    end
-
     if simplify
         sys = mtkcompile(sys)
     end

src/ECMEDox.jl

@@ -1,7 +1,9 @@
-# Good old 16-state CICR (Cortassa 2006), including L-type calcium channel (LCC) and Ryanodine receptor (RyR)
-"LCC ODE system"
-function get_lcc_sys(; ca_ss, ca_o, k_i, k_o, vm, name=:cicrsys)
+# 16-state CICR from Cortassa et al. (2006), including L-type calcium channels (LCC) and Ryanodine receptors (RyR)
+function get_cicr16_eqs(; vm, ca_ss, ca_i, ca_jsr, ca_nsr, k_i, ca_o, k_o)
     @parameters begin
+        R_TR = inv(9.09ms)          # Diffusion rate between JSR and NSR
+        R_XFER= inv(0.5747ms)       # Diffusion rate between subspace and cytosol
+        ## L-type calcium channels
         A_LCC = 2
         B_LCC = 1/2
         ω_LCC = 0.01kHz
@@ -11,6 +13,114 @@ function get_lcc_sys(; ca_ss, ca_o, k_i, k_o, vm, name=:cicrsys)
         P_CA_LCC = 1E-3cm * Hz # 1.24E-3cm * Hz
         P_K_LCC = 1.11E-11cm * Hz
         I_CA_HALF_LCC = -0.4583μAcm⁻²
+        ## Ryanodine receptors
+        R_RYR = 3.6kHz
+        KA_P_RYR = 1.125E10 / (mM^4 * ms)
+        KA_M_RYR = 0.576kHz
+        KB_P_RYR = 4.05E6 / (mM^3 * ms)
+        KB_M_RYR = 1.93kHz
+        KC_P_RYR = 0.1kHz
+        KC_M_RYR = 8E-4kHz
+    end
+
+    @variables begin
+        # State transition rates
+        α_lcc(t)
+        β_lcc(t)
+        # Closed LCC state
+        c0_lcc(t) # Conserved
+        c1_lcc(t) = 0
+        c2_lcc(t) = 0
+        c3_lcc(t) = 0
+        c4_lcc(t) = 0
+        # Opened LCC state
+        o_lcc(t) = 0
+        # Ca-inhibited LCC states
+        cca0_lcc(t) = 0
+        cca1_lcc(t) = 0
+        cca2_lcc(t) = 0
+        cca3_lcc(t) = 0
+        cca4_lcc(t) = 0
+        # Voltage-gated LCC
+        x_yca(t) = 1
+        y_inf(t)
+        τ_yca(t)
+        # RyR
+        po1_ryr(t) = 0
+        po2_ryr(t) = 0
+        pc1_ryr(t) ## Conserved
+        pc2_ryr(t) = 0
+        # Currents
+        ICaMax(t)
+        ICaK(t)
+        ICaL(t)
+        Jrel(t)
+        Jtr(t)
+        Jxfer(t)
+    end
+
+    ω, f, g = ω_LCC, f_LCC, g_LCC
+    α = α_lcc
+    β = β_lcc
+    α′ = A_LCC * α
+    β′ = B_LCC * β
+    γ = γ_LCC * ca_ss
+    rs_lcc = Dict()
+    add_rate!(rs_lcc, 4α, c0_lcc, β, c1_lcc)
+    add_rate!(rs_lcc, 3α, c1_lcc, 2β, c2_lcc)
+    add_rate!(rs_lcc, 2α, c2_lcc, 3β, c3_lcc)
+    add_rate!(rs_lcc, α, c3_lcc, 4β, c4_lcc)
+    add_rate!(rs_lcc, f, c4_lcc, g, o_lcc)
+    add_rate!(rs_lcc, 4α′, cca0_lcc, β′, cca1_lcc)
+    add_rate!(rs_lcc, 3α′, cca1_lcc, 2β′, cca2_lcc)
+    add_rate!(rs_lcc, 2α′, cca2_lcc, 3β′, cca3_lcc)
+    add_rate!(rs_lcc, α′, cca3_lcc, 4β′, cca4_lcc)
+    add_rate!(rs_lcc, γ, c0_lcc, ω, cca0_lcc)
+    add_rate!(rs_lcc, A_LCC * γ, c1_lcc, B_LCC * ω, cca1_lcc)
+    add_rate!(rs_lcc, A_LCC^2 * γ, c2_lcc, B_LCC^2 * ω, cca2_lcc)
+    add_rate!(rs_lcc, A_LCC^3 * γ, c3_lcc, B_LCC^3 * ω, cca3_lcc)
+    add_rate!(rs_lcc, A_LCC^4 * γ, c4_lcc, B_LCC^4 * ω, cca4_lcc)
+    lcceqs = [ D(x) ~ rs_lcc[x] for x in (c1_lcc, c2_lcc, c3_lcc, c4_lcc, o_lcc, cca0_lcc, cca1_lcc, cca2_lcc, cca3_lcc, cca4_lcc) ]
+
+    rs_ryr = Dict()
+    add_rate!(rs_ryr, KA_P_RYR * ca_ss^4, pc1_ryr, KA_M_RYR, po1_ryr)
+    add_rate!(rs_ryr, KB_P_RYR * ca_ss^3, po1_ryr, KB_M_RYR, po2_ryr)
+    add_rate!(rs_ryr, KC_P_RYR, po1_ryr, KC_M_RYR, pc2_ryr)
+    ryreqs = [D(x) ~ rs_ryr[x] for x in (po1_ryr, po2_ryr, pc2_ryr)]
+
+    eqs = [
+        1 ~ c0_lcc + c1_lcc + c2_lcc + c3_lcc + c4_lcc + o_lcc + cca0_lcc + cca1_lcc + cca2_lcc + cca3_lcc + cca4_lcc,
+        1 ~ pc1_ryr + po1_ryr + po2_ryr + pc2_ryr,
+        α_lcc ~ 0.4kHz * exp((vm + 2mV) * inv(10mV)),
+        β_lcc ~ 0.05kHz * exp(-(vm + 2mV) * inv(13mV)),
+        y_inf ~ expit(-(vm + 55mV) * inv(7.5mV)) + 0.5 * expit((vm - 21mV) * inv(6mV)),
+        τ_yca ~ 20ms + 600ms * expit(-(vm + 30mV) * inv(9.5mV)),
+        D(x_yca) ~ (y_inf - x_yca) / τ_yca,
+        ICaMax ~ ghk(P_CA_LCC, vm, 1μM, 0.341 * ca_o, 2) ,
+        ICaL ~ 6 * x_yca * o_lcc * ICaMax,
+        ICaK ~ hil(I_CA_HALF_LCC, ICaMax) * x_yca * o_lcc * ghk(P_K_LCC, vm, k_i, k_o),
+        Jrel ~ R_RYR * (po1_ryr + po2_ryr) * (ca_jsr - ca_ss),
+        Jtr ~ R_TR * (ca_nsr - ca_jsr),
+        Jxfer ~ R_XFER * (ca_ss - ca_i),
+    ]
+
+    eqs_cicr16 = [ryreqs; lcceqs; eqs]
+
+    return (; eqs_cicr16, ICaL, ICaK, Jrel, Jtr, Jxfer)
+end
+
+"LCC ODE system"
+function get_lcc_sys(; ca_ss, ca_o, k_i, k_o, vm, name=:cicrsys)
+    @parameters begin
+        A_LCC = 2
+        B_LCC = 1/2
+        ω_LCC = 10Hz
+        f_LCC = 300Hz
+        g_LCC = 2kHz
+        γ_LCC = 0.1875 / (ms * mM)
+        P_CA_LCC = 1E-3cm * Hz # 1.24E-3cm * Hz
+        P_K_LCC = 1.11E-11cm * Hz
+        I_CA_HALF_LCC = -0.4583μAcm⁻²
     end
 
     @variables begin
@@ -35,13 +145,18 @@ function get_lcc_sys(; ca_ss, ca_o, k_i, k_o, vm, name=:cicrsys)
         x_yca(t) = 1
         y_inf(t)
         τ_yca(t)
+        # RyR
+        po1_ryr(t) = 0
+        po2_ryr(t) = 0
+        pc1_ryr(t) ## Conserved
+        pc2_ryr(t) = 0
         # Currents
         ICaMax(t)
         ICaK(t)
         ICaL(t)
+        Jrel(t)
     end
 
-    v = vm / mV
     ω, f, g = ω_LCC, f_LCC, g_LCC
     α = α_lcc
     β = β_lcc
@@ -64,8 +179,8 @@ function get_lcc_sys(; ca_ss, ca_o, k_i, k_o, vm, name=:cicrsys)
     vc4ca4 = 16γ * c4_lcc - ω / 16 * cca4_lcc
 
     eqs = [
-        α_lcc ~ 0.4kHz * exp((v + 2mV) * inv(10mV)),
-        β_lcc ~ 0.05kHz * exp(-(v + 2mV) * inv(13mV)),
+        α_lcc ~ 0.4kHz * exp((vm + 2mV) * inv(10mV)),
+        β_lcc ~ 0.05kHz * exp(-(vm + 2mV) * inv(13mV)),
         1 ~ c0_lcc + c1_lcc + c2_lcc + c3_lcc + c4_lcc + o_lcc + cca0_lcc + cca1_lcc + cca2_lcc + cca3_lcc + cca4_lcc,
         # D(c0_lcc) ~ -vc0c1 - vc0ca0,
         D(c1_lcc) ~ vc0c1 - vc1c2 - vc1ca1,
@@ -78,8 +193,8 @@ function get_lcc_sys(; ca_ss, ca_o, k_i, k_o, vm, name=:cicrsys)
         D(cca2_lcc) ~ vc2ca2 + vca1ca2 - vca2ca3,
         D(cca3_lcc) ~ vc3ca3 + vca2ca3 - vca3ca4,
         D(cca4_lcc) ~ vca3ca4 + vc4ca4,
-        y_inf ~ expit(-(v + 55mV) * inv(7.5mV)) + 0.5 * expit((v - 21mV) * inv(6mV)),
-        τ_yca ~ 20ms + 600ms * expit(-(v + 30mV) * inv(9.5mV)),
+        y_inf ~ expit(-(vm + 55mV) * inv(7.5mV)) + 0.5 * expit((vm - 21mV) * inv(6mV)),
+        τ_yca ~ 20ms + 600ms * expit(-(vm + 30mV) * inv(9.5mV)),
         D(x_yca) ~ (y_inf - x_yca) / τ_yca,
         ICaMax ~ ghk(P_CA_LCC, vm, 1μM, 0.341 * ca_o, 2) ,
         ICaL ~ 6 * x_yca * o_lcc * ICaMax,

src/cicr16.jl

@@ -1,5 +1,4 @@
-"Creatine kinase (CK)"
-function get_ck_sys(; atp_i, adp_i, name=:cksys)
+function get_ck_eqs(; atp_i, adp_i)
     @parameters begin
         KCK_IN = 0.14Hz     # Rate constant of creatine kinase (cytosolic)
         KCK_MT = 1.33E-3Hz  # Rate constant of creatine kinase (juxta-mitochondrial)
@@ -9,7 +8,6 @@ function get_ck_sys(; atp_i, adp_i, name=:cksys)
         ΣCR = 25mM          # Pool of creatine and creatine phosphate
         ΣA_ic = 8mM         # Pool of cytosolic ATP + ADP
     end
-
     @variables begin
         cr_ic(t)  ## Conserved
         crp_ic(t) = 5.1291mM
@@ -21,7 +19,8 @@ function get_ck_sys(; atp_i, adp_i, name=:cksys)
         vCK_cyto(t)
         vTR_crp(t)
     end
-    eqs = [
+
+    eqs_ck = [
         ΣCR ~ cr_ic + crp_ic,
         ΣCR ~ cr_i + crp_i,
         ΣA_ic ~ adp_ic + atp_ic,
@@ -32,5 +31,12 @@ function get_ck_sys(; atp_i, adp_i, name=:cksys)
         D(crp_i) ~ vCK_mito - vTR_crp,
         D(crp_ic) ~ vTR_crp + vCK_cyto
     ]
-    return System(eqs, t; name)
+
+    return (; eqs_ck, vCK_mito, vCK_cyto, vTR_crp)
+end
+
+"Creatine kinase (CK)"
+function get_ck_sys(; atp_i, adp_i, name=:cksys)
+    @unpack eqs_ck = get_ck_eqs(; atp_i, adp_i)
+    return System(eqs_ck, t; name)
 end

src/ck.jl

@@ -1,5 +1,4 @@
-"Sarcomere force generation and ATP consumption"
-function get_force_sys(; atp_i, adp_i, ca_i, name=:forcesys)
+function get_force_eqs(; atp_i, adp_i, ca_i)
     @parameters begin
         ΣLTRPN = 70μM               # Total pool of low affinity troponin ca binding sites
         ΣHTRPN = 140μM              # Total pool of high affinity troponin ca binding sites
@@ -36,7 +35,7 @@ function get_force_sys(; atp_i, adp_i, ca_i, name=:forcesys)
         G23_SL_AM = Φ_AM * G23_AM
         G01_OFF_AM = Φ_AM * G_OFF
         N_TROP = 3.5 * SL - 2.0   ## Hill coefficient ~ 5.525
-        K_CA_TRPN = K_M_LTRPN / K_P_LTRPN  # Activation factor for calcium of low affinity troponin sites
+        K_CA_TRPN = K_M_LTRPN / K_P_LTRPN  ## Activation factor for calcium of low affinity troponin sites
         K½_TRPN = hil((1.7 - 0.8 * (SL - 1.7) / 0.6), K_CA_TRPN)
     end
 
@@ -58,13 +57,16 @@ function get_force_sys(; atp_i, adp_i, ca_i, name=:forcesys)
     end
 
     k_np_trop = K_PN_TROP * NaNMath.pow(ltr_ca / (ΣLTRPN * K½_TRPN), N_TROP)
-    v01 = F01_AM * x_p0 - G01_SL_AM * x_p1
-    v12 = F12_AM * x_p1 - G12_SL_AM * x_p2
-    v23 = F23_AM * x_p2 - G23_SL_AM * x_p3
-    v04 = K_PN_TROP * x_p0 - k_np_trop * x_n0
-    v15 = K_PN_TROP * x_p1 - k_np_trop * x_n1
-    v54 = G01_OFF_AM * x_n1
 
+    rs = Dict()
+    add_rate!(rs, F01_AM, x_p0, G01_SL_AM, x_p1)
+    add_rate!(rs, F12_AM, x_p1, G12_SL_AM, x_p2)
+    add_rate!(rs, F23_AM, x_p2, G23_SL_AM, x_p3)
+    add_rate!(rs, K_PN_TROP, x_p0, k_np_trop, x_n0)
+    add_rate!(rs, K_PN_TROP, x_p1, k_np_trop, x_n1)
+    add_rate!(rs, G01_OFF_AM, x_n1, 0, x_n0)
+
+    des = [D(x) ~ rs[x] for x in (x_p0, x_p1, x_p2, x_p3, x_n1)]
     eqs = [
         force ~ ζ_AM * (x_p1 + x_n1 + 2 * x_p2 + 3 * x_p3) * iDEN_FORCE,
         force_normal ~ (x_p1 + x_p2 + x_p3 + x_n1) * iDEN_FORCE_N,
@@ -75,11 +77,13 @@ function get_force_sys(; atp_i, adp_i, ca_i, name=:forcesys)
         D(htr_ca) ~ K_P_HTRPN * ca_i * htr_free - K_M_HTRPN * htr_ca,
         Jtrpn ~ D(ltr_ca) + D(htr_ca),
         1 ~ x_p0 + x_p1 + x_p2 + x_p3 + x_n0 + x_n1,
-        D(x_p0) ~ -v01 - v04,
-        D(x_p1) ~ v01 - v12 - v15,
-        D(x_p2) ~ v12 - v23,
-        D(x_p3) ~ v23,
-        D(x_n1) ~ v15 - v54,
     ]
-    return System(eqs, t; name)
+    eqs_force = [des; eqs]
+    return (; eqs_force, vAm, Jtrpn)
+end
+
+"Sarcomere force generation and ATP consumption"
+function get_force_sys(; atp_i, adp_i, ca_i, name=:forcesys)
+    @unpack eqs_force = get_force_eqs(; atp_i, adp_i, ca_i)
+    return System(eqs_force, t; name)
 end