separate out ODESystem construction

Author: wsphillips

demo/prinz_kinetics.jl

@@ -1,11 +1,16 @@
 using Conductor, Unitful, ModelingToolkit
-import Unitful: mV, mS, cm, µm, ms, mM, µM, µA
-import Conductor: Na, K, Ca, Cation, Leak
+import Unitful: mV, mS, cm, µm, ms, mM, µM, µA, nM, K
+import Conductor: mR, ℱ
 
 Vₘ = ParentScope(MembranePotential(-50mV))
+
 Caᵢ = Concentration(Calcium, 0.05µM, dynamic = true)
+Caₒ = Concentration(Calcium, 3.0mM, location = Conductor.Outside)
+minCa = Concentration(Calcium, 1.0nM, name=:minCa)
 ICa = IonCurrent(Calcium, aggregate = true)
 
+T = Temperature(283.15K)
+
 nav_kinetics = [Gate(SteadyStateTau,
                      1.0 / (1.0 + exp((Vₘ + 25.5) / -5.29)),
                      2.64 - (2.52 / (1 + exp((Vₘ + 120.0) / -25.0))),
@@ -57,21 +62,14 @@ h_kinetics = [
          1.0 / (1.0 + exp((Vₘ + 75.0) / 5.5)),
          2 / (exp((Vₘ + 169.7) / (-11.6)) + exp((Vₘ - 26.7) / (14.3))), name = :m)]
 
-Ca_Nernst(Ca) = 1000 * ((-8.314 * 283.15) / (2 * 96485.365)) * log(max(Ca, 0.001) / 3000.0)
-@register_symbolic Ca_Nernst(Ca)
-ModelingToolkit.get_unit(op::typeof(Ca_Nernst), args) = mV
-
 gradients = Equilibria(Pair[Sodium => 50.0mV,
                             Potassium => -80.0mV,
                             Cation => (-20mV, :H),
                             Leak => -50mV,
-                            # x1000 for mV; thresholding Cai in case negative conc
-                            # TODO: write nernst function + use Unitful constants (RT/zℱ)
-                            Calcium => Ca_Nernst(Caᵢ)]);
+                            Calcium => (-mR*T/2ℱ) * log(max(Caᵢ, minCa) / Caₒ)]);
 
 # Reported area value in Prinz 2003; the area of a cylinder without the ends (2πrh).
 # Dimensions r = 25µm, h = 400µm given by Liu et al 1998.
-
 #r = 25µm; h = 400µm; area = round(ustrip(Float64, cm^2, 2π*r*h), sigdigits=3)
 
 geo = Cylinder(radius = 25µm, height = 400µm)

docs/src/core/primitives.md

@@ -12,6 +12,7 @@ ExtrinsicPotential
 IonConcentration
 IonCurrent
 EquilibriumPotential
+Conductor.Temperature
 ```
 
 ## Constants
@@ -21,6 +22,8 @@ Conductor.t
 Conductor.D
 Conductor.ℱ
 Conductor.IonSpecies
+Conductor.R
+Conductor.mR
 ```
 ### Quantities
 

src/Conductor.jl

@@ -95,7 +95,7 @@ export Simulation
 export @named
 
 export Calcium, Sodium, Potassium, Chloride, Cation, Anion, Leak, Ion, NonIonic
-
+export Temperature
 export t
 
 export CompartmentSystem, Compartment, ConductanceSystem, Conductance,

src/primitives.jl

@@ -19,7 +19,25 @@ julia> Conductor.ℱ
 96485.33212331001 C mol^-1
 ```
 """
-const ℱ = Unitful.q * Unitful.Na # Faraday's constant
+const ℱ = let name = :ℱ
+    only(@constants $name=ustrip(Unitful.q * Unitful.Na) [unit = Unitful.C*Unitful.mol^-1])
+end
+
+"""
+Universal Gas Constant (base units)
+"""
+const R = let name = :R
+    only(@constants $name=ustrip(Unitful.R) [unit = Unitful.unit(Unitful.R)])
+end
+
+"""
+Universal Gas Constant (milli-)
+
+Scaled by 1000x. Results in output 
+"""
+const mR = let name = :mR
+    only(@constants $name=(ustrip(Unitful.R)*1000) [unit = u"mJ*K^-1*mol^-1"])
+end
 
 """
 The independent variable for time, ``t``.
@@ -278,6 +296,28 @@ function EquilibriumPotential(ion::IonSpecies, eqv::Union{Nothing, Real, Voltage
     return setmetadata(ret, EquilibriumPotential, EquilibriumPotential(ion))
 end
 
+"""
+    Temperature(temp; <keyword arguments>)
+
+Temperature (in Kelvin)
+
+# Arguments
+- `dynamic::Bool = false`: a dynamic `Temperature` is assumed to vary with time.
+- `name::Symbol = :T`: the symbol to use for the symbolic variable.
+"""
+function Temperature(temp::Union{Nothing, Real, Unitful.Temperature};
+                     dynamic::Bool = false, name::Symbol = :T)
+    temp_val = temp isa Unitful.Temperature ? ustrip(Unitful.K, temp) : temp
+    if isnothing(temp_val)
+        ret = dynamic ? only(@variables $name(t) [unit = Unitful.K]) :
+              only(@parameters $name [unit = Unitful.K])
+    else
+        ret = dynamic ? only(@variables($name(t)=temp_val, [unit = Unitful.K])) :
+              only(@parameters($name=temp_val, [unit = Unitful.K]))
+    end
+    return ret
+end
+
 # Internal API: EquilibriumPotential trait queries
 isreversal(x) = hasmetadata(value(x), EquilibriumPotential)
 getreversal(x) = isreversal(x) ? getmetadata(value(x), EquilibriumPotential) : nothing

test/prinzneuron.jl

@@ -22,13 +22,13 @@ include(normpath(@__DIR__, "..", "demo", "prinz_kinetics.jl"))
                                       extensions = [calcium_conversion]);
     @test length.([equations(neuron),
                       states(neuron),
-                      parameters(neuron)]) == [31, 31, 17]
+                      parameters(neuron)]) == [31, 31, 20]
 
     simul_sys = ODESystem(neuron)
 
     @test length.([equations(simul_sys),
                       states(simul_sys),
-                      parameters(simul_sys)]) == [13, 13, 17]
+                      parameters(simul_sys)]) == [13, 13, 20]
 
     # Prinz STG neuron hand-written reference implementation