diff --git a/demo/pinskyrinzel.jl b/demo/pinskyrinzel.jl new file mode 100644 index 0000000..b7b9a82 --- /dev/null +++ b/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ₘ]) diff --git a/demo/traub_kinetics.jl b/demo/traub_kinetics.jl new file mode 100644 index 0000000..6ca1274 --- /dev/null +++ b/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) + diff --git a/src/Conductor.jl b/src/Conductor.jl index fdf0a12..5da4f61 100644 --- a/src/Conductor.jl +++ b/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) diff --git a/src/channels.jl b/src/channels.jl index 2e24966..ecca236 100644 --- a/src/channels.jl +++ b/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 diff --git a/src/compartments.jl b/src/compartments.jl index 5199f38..f03377c 100644 --- a/src/compartments.jl +++ b/src/compartments.jl @@ -11,6 +11,10 @@ end struct Point <: Geometry end +struct Unitless <: Geometry + value +end + struct Cylinder <: Geometry radius height @@ -23,34 +27,13 @@ end height(x::Geometry) = isdefined(x, :height) ? getfield(x, :height) : nothing radius(x::Geometry) = getfield(x, :radius) -radius(::Point) = 0.0 +radius(::Union{Point,Unitless}) = 0.0 # TODO: Remove unit stripping when we have proper unit checking implemented area(x::Sphere) = ustrip(Float64, cm^2, 4*π*radius(x)^2) area(x::Cylinder) = ustrip(Float64, cm^2, 2*π*radius(x)*(height(x) + (x.open_ends ? 0µm : radius(x)))) area(::Point) = 1.0 - -#= -@enum StimulusTrigger ContinuousTrigger EdgeTriggered - -abstract type Stimulus end - -struct CurrentStimulus - waveform - offset::Current - trigger::StimulusTrigger - delay::TimeF64 - function CurrentStimulus(waveform; offset::Current = 0nA, - trigger::StimulusTrigger = ContinuousTrigger, - delay::Time = 0ms) - return new(waveform, offset, trigger, delay) - end -end - -struct VoltageStimulus end -=# - -# IfElse.ifelse(t > 100.0, IfElse.ifelse(t <= 250.0, amplitude, 0.0), 0.0) +area(x::Unitless) = x.value abstract type AbstractCompartmentSystem <: AbstractTimeDependentSystem end @@ -65,6 +48,7 @@ struct CompartmentSystem <: AbstractCompartmentSystem channel_reversals::Set{Num} synapses::Set{AbstractConductanceSystem} # synaptic conductance systems synaptic_reversals::Set{Num} + axial_conductance::Set{Tuple{AbstractConductanceSystem,Num}} stimuli::Vector{Equation} extensions::Vector{ODESystem} defaults::Dict @@ -73,13 +57,13 @@ struct CompartmentSystem <: AbstractCompartmentSystem systems::Vector{AbstractTimeDependentSystem} observed::Vector{Equation} function CompartmentSystem(iv, voltage, capacitance, geometry, chans, channel_reversals, - synapses, synaptic_reversals, stimuli, extensions, defaults, + synapses, synaptic_reversals, axial_conductance, stimuli, extensions, defaults, name, eqs, systems, observed; checks = false) if checks # placeholder end new(iv, voltage, capacitance, geometry, chans, channel_reversals, synapses, - synaptic_reversals, stimuli, extensions, defaults, name, eqs, systems, observed) + synaptic_reversals, axial_conductance, stimuli, extensions, defaults, name, eqs, systems, observed) end end @@ -102,7 +86,7 @@ function CompartmentSystem( systems = AbstractSystem[] observed = Equation[] return CompartmentSystem(t, Vₘ, cₘ, geometry, Set(channels), Set(reversals), Set(), - Set(), stimuli, extensions, Dict(), name, eqs, systems, observed) + Set(), Set(), stimuli, extensions, Dict(), name, eqs, systems, observed) end # AbstractSystem interface extensions @@ -163,6 +147,7 @@ function build_toplevel!(dvs, ps, eqs, defs, comp_sys::CompartmentSystem) end end + filter!(x -> !isequal(x, t), rev_eq_vars) foreach(x -> isparameter(x) && push!(ps, x), rev_eq_vars) # Gather channel current equations @@ -189,6 +174,18 @@ function build_toplevel!(dvs, ps, eqs, defs, comp_sys::CompartmentSystem) push!(syncurrents, I) merge!(defs, defaults(synapse)) end + + # Gather axial current equations + for (ax, childvm) in get_axial_conductance(comp_sys) + I = IonCurrent(Leak, name = Symbol("I", nameof(ax))) + g = renamespace(ax, get_output(ax)) + # TODO: we should have a metadata check for area-dependent units + push!(eqs, I ~ g*aₘ*(childvm - Vₘ)) + push!(dvs, I) + push!(dvs, childvm) + push!(currents, -I) + merge!(defs, defaults(ax)) + end # Gather extension equations for extension in get_extensions(comp_sys) @@ -200,9 +197,12 @@ function build_toplevel!(dvs, ps, eqs, defs, comp_sys::CompartmentSystem) for stimulus in get_stimuli(comp_sys) I = only(get_variables(stimulus.lhs)) + _vars = filter!(x -> !isequal(x,t), get_variables(stimulus.rhs)) + foreach(x -> push!(isparameter(x) ? ps : dvs, x), _vars) hasdefault(I) || push!(defs, I => stimulus.rhs) if isparameter(I) push!(ps, I) + push!(currents, -I) else push!(eqs, stimulus) push!(dvs, I) @@ -243,12 +243,12 @@ function get_systems(x::AbstractCompartmentSystem; rebuild = false) # collect channels + synapses + input systems #if rebuild || isempty(getfield(x, :systems)) empty!(getfield(x, :systems)) - union!(getfield(x, :systems), getfield(x, :chans), getfield(x, :synapses)) + union!(getfield(x, :systems), getfield(x, :chans), getfield(x, :synapses), first.(getfield(x, :axial_conductance))) #end return getfield(x, :systems) end -function Base.:(==)(sys1::CompartmentSystem, sys2::CompartmentSystem) +function Base.:(==)(sys1::AbstractCompartmentSystem, sys2::AbstractCompartmentSystem) sys1 === sys2 && return true iv1 = get_iv(sys1) iv2 = get_iv(sys2) @@ -297,3 +297,23 @@ function Base.convert(::Type{ODESystem}, compartment::CompartmentSystem) return ODESystem(eqs, t, dvs, ps; systems = syss, defaults = defs, name = nameof(compartment)) end +#= +@enum StimulusTrigger ContinuousTrigger EdgeTriggered + +abstract type Stimulus end + +struct CurrentStimulus + waveform + offset::Current + trigger::StimulusTrigger + delay::TimeF64 + function CurrentStimulus(waveform; offset::Current = 0nA, + trigger::StimulusTrigger = ContinuousTrigger, + delay::Time = 0ms) + return new(waveform, offset, trigger, delay) + end +end + +struct VoltageStimulus end +=# + diff --git a/src/gates.jl b/src/gates.jl index b698376..51795f5 100644 --- a/src/gates.jl +++ b/src/gates.jl @@ -7,7 +7,7 @@ steadystate(x::AbstractGatingVariable) = getfield(x, :steadystate) forward_rate(x::AbstractGatingVariable) = getfield(x, :alpha) reverse_rate(x::AbstractGatingVariable) = getfield(x, :beta) hasexponent(x::AbstractGatingVariable) = hasfield(typeof(x), :p) ? getfield(x, :p) !== one(typeof(x.p)) : false -exponent(x::AbstractGatingVariable) = hasexponent(x) ? getfield(x, :p) : nothing +exponent(x::AbstractGatingVariable) = getfield(x, :p) abstract type GateVarForm end struct AlphaBeta <: GateVarForm end @@ -24,7 +24,7 @@ struct Gate{T<:GateVarForm} <: AbstractGatingVariable p::Real end -function Gate(::Type{AlphaBeta}, x::Num, y::Num, p::Real = 1; name = Base.gensym("GateVar")) +function Gate(::Type{AlphaBeta}, x, y, p = 1; name = Base.gensym("GateVar")) alpha, beta = x, y ss = alpha/(alpha + beta) tau = inv(alpha + beta) @@ -32,7 +32,7 @@ function Gate(::Type{AlphaBeta}, x::Num, y::Num, p::Real = 1; name = Base.gensym Gate{AlphaBeta}(out, alpha, beta, ss, tau, p) end -function Gate(::Type{SteadyStateTau}, x::Num, y::Num, p::Real = 1; name = Base.gensym("GateVar")) +function Gate(::Type{SteadyStateTau}, x, y, p = 1; name = Base.gensym("GateVar")) ss, tau = x, y alpha = ss/tau beta = inv(tau) - alpha @@ -40,12 +40,16 @@ function Gate(::Type{SteadyStateTau}, x::Num, y::Num, p::Real = 1; name = Base.g return Gate{SteadyStateTau}(out, alpha, beta, ss, tau, p) end -function Gate(::Type{SteadyState}, x::Num, p::Real = 1; name = Base.gensym("GateVar")) +function Base.convert(::Type{Gate{SteadyState}}, x::Union{Gate{AlphaBeta},Gate{SteadyStateTau}}) + return Gate(SteadyState, steadystate(x), exponent(x), name = Symbolics.tosymbol(output(x), escape=false)) +end + +function Gate(::Type{SteadyState}, x, p = 1; name = Base.gensym("GateVar")) out = only(@variables $name(t) = x) return Gate{SteadyState}(out, 0, 0, x, 0, p) end -function Gate(::Type{ConstantValue}, x, p::Real = 1; name = Base.gensym("GateVar")) +function Gate(::Type{ConstantValue}, x, p = 1; name = Base.gensym("GateVar")) out = only(@parameters $name = x) return Gate{ConstantValue}(out, 0, 0, x, 0, p) end diff --git a/src/multicompartment.jl b/src/multicompartment.jl new file mode 100644 index 0000000..5d258db --- /dev/null +++ b/src/multicompartment.jl @@ -0,0 +1,130 @@ +abstract type AbstractJunction end + +struct Junction <: AbstractJunction + parent::CompartmentSystem + child::CompartmentSystem + conductance::ConductanceSystem + symmetric::Bool +end + +Junction(x::Pair, cond; symmetric = true) = Junction(x.first, x.second, cond, symmetric) +get_conductance(x::Junction) = getfield(x, :conductance) +issymmetric(x::Junction) = getfield(x, :symmetric) + +struct MultiCompartmentSystem <: AbstractCompartmentSystem + iv::Num + junctions::Vector{AbstractJunction} + compartments::Vector{CompartmentSystem} + extensions::Vector{ODESystem} + eqs::Vector{Equation} + systems::Vector{AbstractTimeDependentSystem} + observed::Vector{Equation} + defaults::Dict + name::Symbol + function MultiCompartmentSystem(iv, junctions, compartments, extensions, eqs, systems, + observed, defaults, name; checks = false) + if checks + #placeholder + end + return new(iv, junctions, compartments, extensions, eqs, systems, observed, defaults, name) + end +end + +const MultiCompartment = MultiCompartmentSystem + +function MultiCompartment(junctions::Vector{<:AbstractJunction}; extensions = ODESystem[], + name = Base.gensym("MultiCompartment")) + + eqs = Equation[] + systems = AbstractTimeDependentSystem[] + observed = Equation[] + defaults = Dict() + all_comp = Set{CompartmentSystem}() + + for jxn in junctions + push!(all_comp, jxn.child) + push!(all_comp, jxn.parent) + end + + return MultiCompartment(t, junctions, collect(all_comp), extensions, eqs, systems, observed, + defaults, name) +end + +get_junctions(x::MultiCompartmentSystem) = getfield(x, :junctions) +get_axial_conductance(x::AbstractCompartmentSystem) = getfield(x, :axial_conductance) +get_compartments(x::MultiCompartmentSystem) = getfield(x, :compartments) + +#MTK.get_states(x::MultiCompartmentSystem) = collect(build_toplevel(x)[1]) +#MTK.has_ps(x::MultiCompartmentSystem) = !isempty(build_toplevel(x)[2]) +#MTK.get_ps(x::MultiCompartmentSystem) = collect(build_toplevel(x)[2]) +# +#function MTK.get_eqs(x::AbstractMultiCompartmentSystem) +# empty!(getfield(x, :eqs)) +# union!(getfield(x, :eqs), build_toplevel(x)[3]) +#end +# +#MTK.get_defaults(x::AbstractMultiCompartmentSystem) = build_toplevel(x)[4] +# +#function MTK.get_systems(x::AbstractMultiCompartmentSystem) +# empty!(getfield(x, :systems)) +# union!(getfield(x, :systems), build_toplevel(x)[5], get_extensions(x)) +#end +# +function get_systems(x::MultiCompartmentSystem; rebuild = false) + empty!(getfield(x, :systems)) + union!(getfield(x, :systems), getfield(x, :compartments), getfield(x, :extensions)) + return getfield(x, :systems) +end + + +function build_toplevel!(dvs, ps, eqs, defs, mcsys::MultiCompartmentSystem) + + junctions = get_junctions(mcsys) + compartments = get_compartments(mcsys) + jxn_types = Set{ConductanceSystem}() + + forwards = Set{Equation}() + + for junction in junctions + push!(jxn_types, get_conductance(junction)) + end + + # Reset subcompartment axial connections + foreach(x -> empty!(get_axial_conductance(x)), compartments) + + for jxn in junctions + axial = get_conductance(jxn) # maybe replicate? needs a toggle + parent = jxn.parent + child = jxn.child + childvm = MembranePotential(; name = Symbol(:V, nameof(child))) + push!(get_axial_conductance(parent), (axial, childvm)) + #FIXME: make this generic for any unresolved state in the conductance sys + if hasproperty(getproperty(parent, nameof(axial)), :Vₘ) + push!(forwards, child.Vₘ ~ getproperty(parent, nameof(axial)).Vₘ) + end + push!(forwards, child.Vₘ ~ getproperty(parent, tosymbol(childvm, escape=false))) + if issymmetric(jxn) + parentvm = MembranePotential(; name = Symbol(:V, nameof(parent))) + push!(get_axial_conductance(child), (axial, parentvm)) + #FIXME: make this generic ofr any unresolved state in the cond. sys + if hasproperty(getproperty(child, nameof(axial)), :Vₘ) + push!(forwards, parent.Vₘ ~ getproperty(child, nameof(axial)).Vₘ) + end + push!(forwards, parent.Vₘ ~ getproperty(child, tosymbol(parentvm, escape=false))) + end + end + + union!(eqs, forwards) + return dvs, ps, eqs, defs, collect(compartments) +end + +#get_geometry(x::MultiCompartmentSystem) +#area(x::MultiCompartmentSystem) + +function Base.convert(::Type{ODESystem}, mcsys::MultiCompartmentSystem) + states, params, eqs, defs, compartments = build_toplevel(mcsys) + + all_systems = map(x -> convert(ODESystem, x), compartments) + odesys = ODESystem(eqs, t, states, params; defaults = defs, name = nameof(mcsys)) + return compose(odesys, all_systems) +end diff --git a/src/networks.jl b/src/networks.jl index 499d46a..f2a5db7 100644 --- a/src/networks.jl +++ b/src/networks.jl @@ -62,7 +62,6 @@ end get_extensions(x::AbstractNeuronalNetworkSystem) = getfield(x, :extensions) get_synapses(x::AbstractNeuronalNetworkSystem) = getfield(x, :synapses) -# Forward getters to internal system MTK.get_states(x::AbstractNeuronalNetworkSystem) = collect(build_toplevel(x)[1]) MTK.has_ps(x::AbstractNeuronalNetworkSystem) = !isempty(build_toplevel(x)[2]) MTK.get_ps(x::AbstractNeuronalNetworkSystem) = collect(build_toplevel(x)[2]) @@ -76,7 +75,7 @@ MTK.get_defaults(x::AbstractNeuronalNetworkSystem) = build_toplevel(x)[4] function MTK.get_systems(x::AbstractNeuronalNetworkSystem) empty!(getfield(x, :systems)) - union!(getfield(x, :systems), build_toplevel(x)[5]) + union!(getfield(x, :systems), build_toplevel(x)[5], get_extensions(x)) end function Base.convert(::Type{ODESystem}, nnsys::NeuronalNetworkSystem) @@ -89,25 +88,19 @@ end function build_toplevel!(dvs, ps, eqs, defs, network_sys::NeuronalNetworkSystem) synapses = get_synapses(network_sys) - preneurons = Set() - postneurons = Set() - synapse_types = Set() - reversals = Set() + neurons = Set() voltage_fwds = Set{Equation}() # Bin the fields for synapse in synapses - push!(preneurons, presynaptic(synapse)) - push!(preneurons, postsynaptic(synapse)) - push!(synapse_types, class(synapse)) - push!(reversals, reversal(synapse)) + push!(neurons, presynaptic(synapse)) + push!(neurons, postsynaptic(synapse)) end # Reset all synaptic information - allneurons = union(preneurons, postneurons) - foreach(x -> empty!(get_synapses(x)), allneurons) - foreach(x -> empty!(get_synaptic_reversals(x)), allneurons) + foreach(x -> empty!(get_synapses(x)), neurons) + foreach(x -> empty!(get_synaptic_reversals(x)), neurons) # Push synaptic information to each neuron for synapse in synapses @@ -122,7 +115,7 @@ function build_toplevel!(dvs, ps, eqs, defs, network_sys::NeuronalNetworkSystem) end union!(eqs, voltage_fwds) - return dvs, ps, eqs, defs, collect(allneurons) + return dvs, ps, eqs, defs, collect(neurons) end #= diff --git a/src/util.jl b/src/util.jl new file mode 100644 index 0000000..7f71c32 --- /dev/null +++ b/src/util.jl @@ -0,0 +1,6 @@ + +using IfElse + +function heaviside(x) + IfElse.ifelse(x > zero(x), one(x), zero(x)) +end