WIP: MultiCompartmentSystem (#26)

Implemented `MultiCompartmentSystem` with demo of a single neuron based on Pinsky and Rinzel 1994 two-compartment model.

demo/pinskyrinzel.jl

@@ -0,0 +1,72 @@
+
+include(joinpath(@__DIR__, "traub_kinetics.jl")
+
+import Unitful: µF, pA, µA, nA
+
+@parameters ϕ = 0.13 β = 0.075 p = 0.5
+
+@named calcium_conversion = ODESystem([D(Caᵢ) ~ -ϕ*ICa - β*Caᵢ]);
+
+reversals = Equilibria(Pair[Sodium    =>  120.0mV,
+                            Potassium =>  -15.0mV,
+                            Leak      =>    0.0mV,
+                            Calcium   =>  140.0mV]);
+
+capacitance = 3.0µF/cm^2
+gc_val = 2.1mS/cm^2
+
+# Pinsky modifies NaV to have instantaneous activation, so we can ignore tau
+pinsky_nav_kinetics = [convert(Gate{SteadyState}, nav_kinetics[1]), nav_kinetics[2]]
+@named NaV = IonChannel(Sodium, pinsky_nav_kinetics) 
+
+# No inactivation term for calcium current in Pinsky model
+pinsky_ca_kinetics = [ca_kinetics[1]]
+@named CaS = IonChannel(Calcium, pinsky_ca_kinetics)
+
+is_val = ustrip(Float64, µA, 0.75µA)/p
+@named Iₛ = IonCurrent(NonIonic, is_val, dynamic = false)
+soma_holding = Iₛ ~ is_val
+
+id_val = ustrip(Float64, µA, 0.0µA)/(1-p)
+@named I_d = IonCurrent(NonIonic, id_val, dynamic = false)
+dendrite_holding = I_d ~ id_val
+
+@named soma = Compartment(Vₘ,
+                         [NaV(30mS/cm^2),
+                          Kdr(15mS/cm^2),
+                          leak(0.1mS/cm^2)],
+                          reversals[1:3],
+                          geometry = Unitless(0.5), # FIXME: should take p param
+                          capacitance = capacitance,
+                          stimuli = [soma_holding])
+
+@named dendrite = Compartment(Vₘ,
+                             [KAHP(0.8mS/cm^2),
+                              CaS(10mS/cm^2),
+                              KCa(15mS/cm^2),
+                              leak(0.1mS/cm^2)],
+                             reversals[2:4],
+                             geometry = Unitless(0.5), # FIXME: should take p param
+                             capacitance = capacitance,
+                             extensions = [calcium_conversion],
+                             stimuli = [dendrite_holding])
+
+@named gc_soma = AxialConductance([Gate(SteadyState, 1/p, name = :ps)],
+                                  max_g = gc_val)
+@named gc_dendrite = AxialConductance([Gate(SteadyState, 1/(1-p), name = :pd)],
+                                      max_g = gc_val)
+soma2dendrite = Junction(soma => dendrite, gc_soma)
+dendrite2soma = Junction(dendrite => soma, gc_dendrite)
+
+mc_neuron = MultiCompartment([soma2dendrite, dendrite2soma])
+
+simp = Simulation(mc_neuron, time=2000ms, return_system = true)
+prob = ODAEProblem(simp, [], (0., 2000), [])
+# Uncomment to explicitly use the same u0 as published
+# prob = ODAEProblem(simp, [-4.6, 0.999, 0.001, 0.2, -4.5, 0.01, 0.009, .007], (0., 2000), [])
+
+using OrdinaryDiffEq, Plots
+# Pinsky & Rinzel originally solved using RK4 and dt=0.05
+#sol = solve(prob, RK4(), dt=0.05, maxiters=1e9)
+sol = solve(prob, Rosenbrock23())
+plot(sol, vars=[soma.Vₘ])

demo/traub_kinetics.jl

@@ -0,0 +1,92 @@
+using Conductor, Unitful, ModelingToolkit, IfElse
+import Unitful: mV, mS, cm, µm, ms, mM, µM
+import Conductor: Na, K, Ca, Cation, Leak
+
+Vₘ = MembranePotential(-4.6mV)
+Caᵢ = Concentration(Calcium, 0.2µM, dynamic = true)
+ICa = IonCurrent(Calcium, aggregate = true)
+#@parameters ϕ β = 0.075
+#@named calcium_conversion = ODESystem([D(Caᵢ) ~ -ϕ*ICa - β*Caᵢ]);
+
+# Sodium channels
+nav_kinetics = [
+    # Activation
+    Gate(AlphaBeta,
+         (0.32(13.1 - Vₘ))/(exp((13.1 - Vₘ)/4.)-1.),
+         (0.28(Vₘ - 40.1))/(exp((Vₘ-40.1)/5.)-1.),
+         2,
+         name = :m)
+    # Inactivation
+    Gate(AlphaBeta,
+         0.128*exp((17. - Vₘ)/18.),
+         4. /(1. + exp((40. - Vₘ)/5.)),
+         name = :h)]
+
+# Calcium channels
+αᵣ = IfElse.ifelse(Vₘ <= zero(Float64), 0.005, exp(-Vₘ/20.0)/200.0)
+
+ca_kinetics = [
+    # Activation
+    Gate(AlphaBeta,
+         1.6/(1. + exp(-0.072*(Vₘ - 65.0))),
+         (0.02*(Vₘ - 51.1))/(exp((Vₘ - 51.1)/5.) - 1.),
+         2,
+         name = :s)
+    # Inactivation; unused in Pinsk & Rinzel reduction
+    Gate(AlphaBeta,
+         αᵣ,
+         IfElse.ifelse(Vₘ <= zero(Float64), zero(Float64), 0.005 - αᵣ),
+         name = :r)]
+# Delayed rectifier potassium
+kdr_kinetics = [
+    Gate(AlphaBeta,
+         (0.016*(35.1 - Vₘ))/(exp((35.1 - Vₘ)/5.)-1.),
+         0.25*exp((20.0 - Vₘ)/40.0),
+         name = :n)]
+
+# After-hyperpolarization potassium
+# Note: The Traub model calls Calcium concentration 'χ' 
+kahp_kinetics = [
+    Gate(AlphaBeta,
+         min(0.00002*Caᵢ, 0.01),
+         0.001,
+         name = :q)]
+
+αc = IfElse.ifelse(Vₘ <= 50,
+                   (exp((Vₘ-10.)/11.)-exp((Vₘ-6.5)/27.))/18.975,
+                   2. *exp((6.5-Vₘ)/27.))
+
+# Calcium-dependent potassium
+kca_kinetics = [
+    # Activation
+    Gate(AlphaBeta,
+         αc,
+         IfElse.ifelse(Vₘ <= 50,
+                       2. *exp((6.5-Vₘ)/27.) - αc,
+                       zero(Float64)),
+         name = :c),
+    # Calcium saturation term
+    Gate(SteadyState, min(Caᵢ/250., 1.),
+         name = :χ)]
+
+# A-type Potassium
+ka_kinetics = [
+    # Activation
+    Gate(AlphaBeta,
+         (0.02*(13.1 - Vₘ))/(exp((13.1-Vₘ)/10.)-1.),
+         (0.0175(Vₘ / 40.1))/(exp((Vₘ-40.1)/10.)-1.),
+         name = :a)
+    # Inactivation
+    Gate(AlphaBeta,
+         0.0016*exp((-13. - Vₘ)/18.),
+         0.05/(1+exp((10.1 - Vₘ)/5.)),
+         name = :b)]
+
+@named NaV = IonChannel(Sodium, nav_kinetics)
+@named CaS = IonChannel(Calcium, ca_kinetics)
+@named Kdr = IonChannel(Potassium, kdr_kinetics)
+@named KAHP = IonChannel(Potassium, kahp_kinetics)
+@named KA  = IonChannel(Potassium, ka_kinetics)
+@named KCa = IonChannel(Potassium, kca_kinetics)
+@named leak = IonChannel(Leak)
+

src/Conductor.jl

@@ -40,7 +40,6 @@ import ModelingToolkit:
     getdefault,
     setdefault,
     AbstractTimeDependentSystem,
-    AbstractSystem,
     independent_variables,
     get_variables!
 
@@ -57,19 +56,33 @@ import SymbolicUtils: FnType
 import Unitful: mV, mS, cm, µF, mF, µm, pA, nA, mA, µA, ms, mM, µM
 import Base: show, display
 
-export Gate, AlphaBeta, SteadyStateTau, IonChannel, PassiveChannel, SynapticChannel, Synapse
-export EquilibriumPotential, Equilibrium, Equilibria, MembranePotential, IonCurrent
-export AuxConversion, D, NeuronalNetwork
-export Simulation, Concentration, IonConcentration
+export Gate, AlphaBeta, SteadyStateTau, SteadyState, ConstantValue,
+       IonChannel, PassiveChannel, SynapticChannel, Synapse, Junction,
+       AxialConductance
+
+export EquilibriumPotential, Equilibrium, Equilibria, MembranePotential,
+       IonCurrent, IonConcentration, Concentration
+
+export AuxConversion, D
+export Simulation
 export @named
+
 export Calcium, Sodium, Potassium, Chloride, Cation, Anion, Leak, Ion, NonIonic
+
 export t
-export CompartmentSystem, ConductanceSystem, NeuronalNetworkSystem, Conductance, Compartment
-export output, get_output, timeconstant, steadystate, forward_rate, reverse_rate, hasexponent, exponent
-export Sphere, Cylinder, Point, area, radius, height
+
+export CompartmentSystem, Compartment,
+       ConductanceSystem, Conductance, 
+       NeuronalNetworkSystem, NeuronalNetwork,
+       MultiCompartmentSystem, MultiCompartment
+
+export output, get_output, timeconstant, steadystate, forward_rate,
+       reverse_rate, hasexponent, exponent
+
+export Sphere, Cylinder, Point, Unitless, area, radius, height
 
 const ℱ = Unitful.q*Unitful.Na # Faraday's constant
-const t = let name = :t; only(@parameters $name) end
+const t = let name = :t; only(@variables $name) end
 const D = Differential(t)
 const ExprValues = Union{Expr,Symbol,Number}  # For use in macros
 
@@ -106,10 +119,12 @@ include("ions.jl")
 include("gates.jl")
 include("channels.jl")
 include("compartments.jl")
+include("multicompartment.jl")
 include("networks.jl")
 include("io.jl")
+include("util.jl")
 
-function Simulation(neuron::CompartmentSystem; time::Time, return_system = false)
+function Simulation(neuron::AbstractCompartmentSystem; time::Time, return_system = false)
     odesys = convert(ODESystem, neuron)
     t_val = ustrip(Float64, ms, time)
     simplified = structural_simplify(odesys)

src/channels.jl

@@ -9,7 +9,6 @@ get_inputs(x::AbstractConductanceSystem) = getfield(x, :inputs)
 get_output(x::AbstractConductanceSystem) = getfield(x, :output)
 # Abstract types without parametrics
 struct ConductanceSystem{S<:AbstractTimeDependentSystem} <: AbstractConductanceSystem
-
     iv
     output::Num # 'g' by default 
     ion::IonSpecies # ion permeability
@@ -104,7 +103,8 @@ function get_eqs(var::Gate{<:Union{AlphaBeta,SteadyStateTau}})
     return [D(x) ~ inv(τₓ)*(x∞ - x)]
 end
 
-get_eqs(var::Gate{<:Union{SteadyState,ConstantValue}}) = [output(var) ~ steadystate(var)]
+get_eqs(var::Gate{SteadyState}) = [output(var) ~ steadystate(var)]
+get_eqs(var::Gate{ConstantValue}) = Equation[]
 
 function ConductanceSystem(g::Num, ion::IonSpecies, gate_vars::Vector{<:AbstractGatingVariable};
         gbar::Num, linearity::IVCurvature = Linear, extensions::Vector{ODESystem} = ODESystem[],
@@ -176,12 +176,22 @@ function IonChannel(ion::IonSpecies,
                       extensions = extensions, linearity = linearity)
 end
 
+function AxialConductance(gate_vars::Vector{<:AbstractGatingVariable} = AbstractGatingVariable[];
+                          max_g = 0mS/cm^2, extensions::Vector{ODESystem} = ODESystem[],
+                          name::Symbol = Base.gensym("Axial"), defaults = Dict())
+
+    IonChannel(Leak, gate_vars, max_g = max_g, extensions = extensions, name = name,
+               defaults = defaults)
+end
+
 function SynapticChannel(ion::IonSpecies,
-        gate_vars::Vector{<:AbstractGatingVariable} = AbstractGatingVariable[];
+                         gate_vars::Vector{<:AbstractGatingVariable} = AbstractGatingVariable[];
                          max_s::Union{Num, ElectricalConductance} = 0mS,
                          extensions::Vector{ODESystem} = ODESystem[],
                          name::Symbol = Base.gensym("SynapticChannel"),
                          linearity::IVCurvature = Linear, defaults = Dict())
+    # to make generic, check for <:Quantity then write a
+    # unit specific "strip" method 
     if max_s isa ElectricalConductance
         sbar_val = ustrip(Float64, mS, max_s)
         @parameters sbar