diff --git a/src/Conductor.jl b/src/Conductor.jl index 77b006d..f72f381 100644 --- a/src/Conductor.jl +++ b/src/Conductor.jl @@ -35,7 +35,7 @@ getdefault(x::Num) = getdefault(ModelingToolkit.value(x)) # Basic symbols const t = toparam(Num(Sym{Real}(:t))) const D = Differential(t) -const MembranePotential() = symoft(:Vₘ) +MembranePotential() = symoft(:Vₘ) @enum Location Outside Inside # Custom Unitful.jl quantities @@ -241,7 +241,7 @@ function IonChannel(conducts::Type{I}, # retrieve all variables present in the RHS of kinetics equations # g = total conductance (e.g. g(m,h) ~ ̄gm³h) if passive - @variables g() # uncalled + @variables g() else for i in gate_vars syms = value.(get_variables(getequation(i))) @@ -259,7 +259,7 @@ function IonChannel(conducts::Type{I}, if passive eqs = [g ~ gbar] - push!(defaultmap, g => gbar_val) # this might be redundant + push!(defaultmap, g => gbar_val) # perhaps redundant else geq = [g ~ gbar * prod(hasexponent(x) ? getsymbol(x)^x.p : getsymbol(x) for x in gate_vars)] eqs = [[getequation(x) for x in gate_vars] @@ -290,13 +290,13 @@ function PassiveChannel(conducts::Type{I}, max_g::SpecificConductance = 0mS/cm^2 end struct SynapticChannel <: AbstractConductance - gbar::ElectricalConductance # scaling term - maximal conductance per unit area - conducts::DataType # ion permeability + gbar::ElectricalConductance + conducts::DataType reversal::Num - inputs::Vector{Num} # cell states dependencies (input args to kinetics); we can infer this + inputs::Vector{Num} params::Vector{Num} - kinetics::Vector{<:AbstractGatingVariable} # gating functions; none = passive channel - sys::Union{ODESystem, Nothing} # symbolic system + kinetics::Vector{<:AbstractGatingVariable} + sys::Union{ODESystem, Nothing} end function SynapticChannel(conducts::Type{I}, @@ -328,17 +328,12 @@ function SynapticChannel(conducts::Type{I}, gates inputs] - #if passive - # eqs = [g ~ gbar] - # push!(defaultmap, g => gbar_val) # this might be redundant - #else geq = [g ~ gbar * prod(hasexponent(x) ? getsymbol(x)^x.p : getsymbol(x) for x in gate_vars)] - eqs = [[getequation(x) for x in gate_vars] - geq] - for x in gate_vars - push!(defaultmap, getsymbol(x) => #= hassteadystate(x) ? x.ss :=# 0.0) # fallback to zero - end - #end + eqs = [[getequation(x) for x in gate_vars] + geq] + for x in gate_vars + push!(defaultmap, getsymbol(x) => 0.0) + end system = ODESystem(eqs, t, states, params; defaults = defaultmap, name = name) @@ -365,35 +360,39 @@ struct Compartment{G} sys::ODESystem end -# TODO: input validation - check channels/reversals for duplicates/conflicts function Compartment{Sphere}(channels::Vector{<:AbstractConductance}, gradients; name::Symbol, area::Float64 = 0.628e-3, #radius = 20µm, capacitance::SpecificCapacitance = 1µF/cm^2, V0::Voltage = -65mV, - applied::Current = 0nA, + holding::Current = 0nA, + stimulus::Union{Function,Nothing} = nothing, aux::Union{Nothing, Vector{AuxConversion}} = nothing) Vₘ = MembranePotential() @variables Iapp(t) Isyn(t) params = @parameters cₘ aₘ grad_meta = getmetadata.(gradients, ConductorEquilibriumCtx) - #r_val = ustrip(Float64, cm, radius) + #r_val = ustrip(Float64, cm, radius) # FIXME: make it so we calculate area from dims as needed systems = [] eqs = Equation[] # equations must be a vector required_states = [] # states not produced or intrinsic (e.g. not currents or Vm) states = Any[Vₘ, Iapp, Isyn] # grow this as we discover/generate new states currents = [] - defaultmap = Pair[Iapp => ustrip(Float64, µA, applied), + defaultmap = Pair[Iapp => ustrip(Float64, µA, holding), aₘ => area, Vₘ => ustrip(Float64, mV, V0), Isyn => 0, cₘ => ustrip(Float64, mF/cm^2, capacitance)] - # By default, we say applied current and synaptic current are constant - # They are both optionally revised to a new value during a "Stimulus" pass or during the - # Network constructor - append!(eqs, [D(Iapp) ~ 0]) + # By default, applied current is constant (just a bias/offset/holding current) + # TODO: lift "stimulus" to a pass that happens at a higher level (i.e. after neuron + # construction; with its own data type) + if stimulus == nothing + append!(eqs, [D(Iapp) ~ 0]) + else + push!(eqs, Iapp ~ stimulus(t,Iapp)) + end # auxillary state transformations (e.g. net calcium current -> Ca concentration) if aux !== nothing @@ -438,7 +437,6 @@ function Compartment{Sphere}(channels::Vector{<:AbstractConductance}, end # write the current equation state - # I = symoft(Symbol(:I,nameof(sys))) # alternatively "uncalled" term I = MembraneCurrent{ion}(name = nameof(sys), aggregate = false) push!(states, I) push!(currents, I) @@ -468,7 +466,7 @@ function Compartment{Sphere}(channels::Vector{<:AbstractConductance}, end end end - # FIXME: handle aggregate current, handle all states + required_states = unique(value.(required_states)) states = Any[unique(value.(states))...] @@ -491,7 +489,6 @@ function Compartment{Sphere}(channels::Vector{<:AbstractConductance}, push!(eqs, vm_eq) system = ODESystem(eqs, t, states, params; systems = systems, defaults = defaultmap, name = name) - # Switch to outputting a compartment/neuron when we build-out networks/synapses return Compartment{Sphere}(capacitance, channels, states, params, system) end @@ -507,7 +504,7 @@ function Network(neurons, topology; name = :Network) all_neurons = Set(getproperty.(neurons, :sys)) eqs = Equation[] params = Set() - states = Set() # placeholder for now + states = Set() # place holder defaultmap = Pair[] systems = [] diff --git a/test/STGnetwork.jl b/test/STGnetwork.jl index 0d5a963..c01b628 100644 --- a/test/STGnetwork.jl +++ b/test/STGnetwork.jl @@ -13,7 +13,7 @@ using OrdinaryDiffEq, Plots Kdr( 50mS/cm^2), H( .01mS/cm^2), leak( 0mS/cm^2)], - gradients, area = area, V0 = -50mV, applied = 0nA, aux = [calcium_conversion]); + gradients, area = area, V0 = -50mV, aux = [calcium_conversion]); # LP 4 @named LP = Soma([NaV( 100mS/cm^2), @@ -24,7 +24,7 @@ using OrdinaryDiffEq, Plots Kdr( 25mS/cm^2), H( .05mS/cm^2), leak(.03mS/cm^2)], - gradients, area = area, V0 = -50mV, applied = 0nA, aux = [calcium_conversion]); + gradients, area = area, V0 = -50mV, aux = [calcium_conversion]); # PY 1 @named PY = Soma([NaV( 100mS/cm^2), @@ -35,7 +35,7 @@ using OrdinaryDiffEq, Plots Kdr( 125mS/cm^2), H( .05mS/cm^2), leak(.01mS/cm^2)], - gradients, area = area, V0 = -50mV, applied = 0nA, aux = [calcium_conversion]); + gradients, area = area, V0 = -50mV, aux = [calcium_conversion]); topology = [ABPD => (LP, Glut(30nS)), ABPD => (LP, Chol(30nS)), diff --git a/test/STGneuron.jl b/test/STGneuron.jl index 3b1c37a..b0602c7 100644 --- a/test/STGneuron.jl +++ b/test/STGneuron.jl @@ -182,7 +182,7 @@ channels = [NaV(300mS/cm^2), =# @named neuron = Soma(channels, gradients, - area = area, V0 = -50mV, applied = 0nA, aux = [calcium_conversion]); + area = area, V0 = -50mV, aux = [calcium_conversion]); t = 5000 sim = Simulation(neuron, time = t*ms) diff --git a/test/hodgkinhuxley.jl b/test/hodgkinhuxley.jl index fa410b9..f6584a8 100644 --- a/test/hodgkinhuxley.jl +++ b/test/hodgkinhuxley.jl @@ -31,9 +31,18 @@ gradients = Equilibria([Na => 50.0mV, area = 4*pi*(20µm)^2 -@named neuron = Soma([NaV,Kdr,leak], gradients, applied = 400pA, area = ustrip(Float64, cm^2, area)) - -t = 250 +function pulse(t, current) + if 100 < t < 200 + return ustrip(Float64, µA, 400pA) + else + return 0.0 + end +end +@register pulse(a,b) + +@named neuron = Soma([NaV,Kdr,leak], gradients, stimulus = pulse, area = ustrip(Float64, cm^2, area)) + +t = 300 sim = Simulation(neuron, time = t*ms) solution = solve(sim, Rosenbrock23()) diff --git a/test/simplesynapse.jl b/test/simplesynapse.jl index 180b5c2..0af6d21 100644 --- a/test/simplesynapse.jl +++ b/test/simplesynapse.jl @@ -32,8 +32,8 @@ gradients = Equilibria([Na => 50.0mV, #area = 4*pi*(10µm)^2 area = 0.629e-3cm^2 -@named neuron1 = Soma([NaV,Kdr,leak], gradients, applied = 5000pA, area = ustrip(Float64, cm^2, area)) -@named neuron2 = Soma([NaV,Kdr,leak], gradients, applied = 0pA, area = ustrip(Float64, cm^2, area)) +@named neuron1 = Soma([NaV,Kdr,leak], gradients, holding = 5000pA, area = ustrip(Float64, cm^2, area)) +@named neuron2 = Soma([NaV,Kdr,leak], gradients, area = ustrip(Float64, cm^2, area)) # Synaptic model syn∞ = 1/(1 + exp((-35 - Vₘ)/5)) @@ -43,8 +43,7 @@ EGlut = Equilibrium{Conductor.Mixed}(0mV, :Glut) # NOTE: -70mV reversal => IPSP @named Glut = Conductor.SynapticChannel(Leak, [syn_kinetics], EGlut, 30nS) -topology = [neuron1 => (neuron2, Glut), - neuron1 => (neuron2, Glut)] +topology = [neuron1 => (neuron2, Glut)] network = Conductor.Network([neuron1, neuron2], topology)