Skip to content

Commit

Permalink
simple pulse stimulus
Browse files Browse the repository at this point in the history
wsphillips committed Jun 26, 2021
1 parent b1c8624 commit 0f7ebe7
Showing 5 changed files with 46 additions and 41 deletions.
57 changes: 27 additions & 30 deletions src/Conductor.jl
Original file line number Diff line number Diff line change
@@ -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 = []

6 changes: 3 additions & 3 deletions test/STGnetwork.jl
Original file line number Diff line number Diff line change
@@ -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)),
2 changes: 1 addition & 1 deletion test/STGneuron.jl
Original file line number Diff line number Diff line change
@@ -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)
15 changes: 12 additions & 3 deletions test/hodgkinhuxley.jl
Original file line number Diff line number Diff line change
@@ -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())
7 changes: 3 additions & 4 deletions test/simplesynapse.jl
Original file line number Diff line number Diff line change
@@ -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)

0 comments on commit 0f7ebe7

Please sign in to comment.