wip

src/CaMKIIModel.jl

@@ -1,12 +1,7 @@
 module CaMKIIModel
-
-export make_b1ar_rn
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+using NaNMath
 
 include("utils.jl")
-include("bar.jl")
-include("cam.jl")
-include("ecc/ica_markov.jl")
-include("ecc/ina_hh.jl")
-include("ecc/ina_markov.jl")
-
-end # module CaMKIIModel
+end

src/neonatal.jl

@@ -0,0 +1,18 @@
+# Neonatal mouse CMC model
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
+include("neonatal_ca_pde.jl")
+
+function build_neonatal_model()
+
+    @constants begin
+        Vmyo = 2.1454e-11           # [L]
+        Vdyad = 1.7790e-014         # [L]
+        VSL = 6.6013e-013           # [L]
+        kSLmyo = 8.587e-15          # [L/msec]
+        CaMKIItotDyad = 120         # [uM]
+        BtotDyad = 1.54 / 8.293e-4  # [uM]
+    end
+
+end

src/neonatal_ca_pde.jl

@@ -0,0 +1,33 @@
+# Cell geometry for calcium diffusion in neonatal mouse CMC model
+using ModelingToolkit
+
+# Calcium nuffer term
+function beta_cai(Cai, TnI_PKAp)
+    fracTnIpo = 0.062698  # Derived quantity (TnI_PKAp(baseline)/TnItot)
+    fPKA_TnI = (1.61 - 0.61 * (1 - TnI_PKAp) / (1 - fracTnIpo)) # Max effect +61#
+    TrpnTotal = 35
+    KmTrpn = 0.5 / fPKA_TnI
+    CmdnTotal = 50
+    KmCmdn = 2.38
+    return inv(1 + TrpnTotal * KmTrpn / (Cai + KmTrpn)^2 + CmdnTotal * KmCmdn / (Cai + KmCmdn)^2)
+end
+
+function get_capde_eqs(;dx = 0.1, rSR_true = 6, rSL_true = 10.5)
+    # Ca diffusion grid
+    rSR = rSR_true + 0.5 * dx
+    rSL = rSL_true - 0.5 * dx
+    j = round(rSR / dx):1:round(rSL / dx) # Spatial indices
+    m = length(j)
+    @variables t Cai(t)[1:m] Cai_mean(t) Cai_sub_SR(t) Cai_sub_SL(t)
+    @parameters Dca = 7 # mum^2/ms set to achive correct diff. speed 0.31 mum/ms
+
+    eqs = [
+        Cai_mean ~ sum(collect(Cai))/m,
+        Cai_sub_SR ~ Cai[1],
+        Cai_sub_SL ~ Cai[m],
+
+    ]
+
+
+    return eqs
+end

src/utils.jl

@@ -1,29 +1,71 @@
-# Constants
-const R = 8314      # [J/kmol*K]
-const Frdy = 96485  # [C/mol]
-const Temp = 310    # [K] 310 K (37 C) for BT / 295 K (22 C) for RT
-const FoRT = Frdy / R / Temp
-const VT = inv(FoRT) # Thermal temparature (mV)
-const Qpow = (Temp - 310) / 10 # Temp factor
+#===
+Constants
+===#
 
+# Units
+const second = 1           # second
+const minute = 60second    # minute
+const ms = 1e-3second      # millisecond
+const Hz = inv(second)     # Herz
+const kHz = 1e3Hz          # kilohertz
+const metre = 1            # meter
+const cm = 0.01metre       # centimeter
+const cm² = cm^2           # square centimeter
+const μm = 1E-6metre       # Micrometer
+const mL = cm^3            # milliliter = cubic centimeter
+const Liter = 1e3mL        # liter
+const μL = 1E-6Liter
+const pL = 1E-12Liter
+const mM = 1
+const Molar = 1000mM       # molar (1000 since the SI units is mM)
+const μM = 1E-3mM          # micromolar
+const nM = 1E-6mM          # nanomolar
+const Amp = 1              # ampere
+const mA = 1E-3Amp         # milliampere
+const μA = 1E-6Amp         # micrpampere
+const Volt = 1             # volt
+const mV = 1E-3Volt        # millivolt
+const mS = mA / Volt       # milliseimens
+const T₀ = 310.0           # Default temp (37C)
+const F = 96485.0          # Faraday constant (columb / mol)
+const μF = 1E-6F
+const R = 8.314            # Ideal gas constant
+const VT = R * T₀ / F      # Thermal voltage (@37C), around 26.7 mV
+const iVT = inv(VT)        # Reciprocal of thermal voltage (@37C)
 # Utility functions
 
-"""Hill function"""
-hil(x, k=one(x)) = x / (x + k)
-hil(x, k, n) = hil(x^n, k^n)
+"""
+Regular Hill/MM function
+"""
+hil(x, k = one(x)) = x / (x + k)
+hil(x, k, n) = hil(NaNMath.pow(x, n), NaNMath.pow(k, n))
 
 """
-Relative exponential function.
-Returns one when x is close to zero
+Repressive Hill/MM function
 """
-function exprel(x)
-    res = x / expm1(x)
-    return ifelse(isapprox(x, zero(x)), one(res), res)
-end
+hilr(x, k = one(x)) = hil(k, x)
+hilr(x, k, n) = hil(k, x, n)
 
-"""Logistic function"""
-expit(x) = inv(one(x) + exp(-x))
+"""
+Logistic sigmoid function.
+"""
+expit(x) = hilr(exp(-x))
+
+"""
+    exprel(x, em1 = expm1(x))
+"""
+exprel(x, em1 = expm1(x)) = x / em1
 
 """Nernst potential"""
-nernst(xo, xi) = VT * log(xo / xi)
-nernst(xo, xi, z) = nernst(xo, xi) / z
+nernst(x_out, x_in) = VT * NaNMath.log(x_out / x_in)
+nernst(x_out, x_in, z) = nernst(x_out, x_in) / z
+
+"GHK flux equation"
+ghk(px, x_i, x_o, zvfrt, ezvfrtm1 = expm1(zvfrt), z = 1) = px * z * F * ((ezvfrtm1 + 1) * x_i - x_o) * exprel(zvfrt, ezvfrtm1)
+
+"GHK flux equation from voltage across the membrane"
+function ghkVm(px, vm, x_i, x_o, z = 1)
+    zvfrt = z * vm * iVT
+    em1 = expm1(zvfrt)
+    return ghk(px, x_i, x_o, zvfrt, em1, z)
+end