Stimulus objects

demo/hodgkinhuxley.jl

@@ -26,18 +26,14 @@ kdr_kinetics = [
 channels = [NaV, Kdr, leak];
 reversals = Equilibria([Na => 50.0mV, K => -77.0mV, Leak => -54.4mV])
 
-@named Iₑ = IonCurrent(NonIonic)
-@named I_rest = IonCurrent(NonIonic, 0.0µA, dynamic = false)
-@named I_step = IonCurrent(NonIonic, 400.0pA, dynamic = false)
-@parameters tstart = 100.0 [unit=ms] tstop = 200.0 [unit=ms]
-electrode_pulse = Iₑ ~ ifelse((t > tstart) & (t < tstop), I_step, I_rest)
+@named Iₑ = PulseTrain(amplitude = 400.0pA, duration = 100ms, delay = 100ms)
 
 dynamics = HodgkinHuxley(Vₘ, channels, reversals;
                          geometry = Sphere(radius = 20µm),
-                         stimuli = [electrode_pulse]);
+                         stimuli = [Iₑ]);
 
 @named neuron = Compartment(dynamics)
 sim = Simulation(neuron, time = 300ms)
-solution = solve(sim, Rosenbrock23(), abstol=0.01, reltol=0.01, saveat=0.2);
+solution = solve(sim, Rosenbrock23(), abstol = 0.01, reltol = 0.01, dtmax = 100.0);
 plot(solution; size=(1200,800))
 

demo/pinskyrinzel.jl

@@ -6,22 +6,17 @@ include(joinpath(@__DIR__, "pinsky_setup.jl"))
 sim_time = 1500.0
 
 # Compartment holding currents
-@named I_s_holding = IonCurrent(NonIonic, -0.5µA, dynamic = false)
-@named Iₛ = IonCurrent(NonIonic, -0.5µA, dynamic = false)
-soma_holding = Iₛ ~ I_s_holding/p
-
-@named I_d_holding = IonCurrent(NonIonic, 0.0µA, dynamic = false)
-@named I_d = IonCurrent(NonIonic, 0.0µA, dynamic = false)
-dendrite_holding = I_d ~ I_d_holding/(1-p)
+@named Is_holding = Bias(-0.5µA / 0.5) # FIXME: we should be able to use 'p' parameter
+@named Id_holding = Bias(0.0µA / (1-0.5)) 
 
 soma_dynamics = HodgkinHuxley(Vₘ,
                          [NaV(30mS/cm^2),
                           Kdr(15mS/cm^2),
                           leak(0.1mS/cm^2)],
                           reversals[1:3];
-                          geometry = Unitless(0.5),
+                          geometry = Unitless(0.5), # FIXME: support 'p' parameter value
                           capacitance = capacitance,
-                          stimuli = [soma_holding]);
+                          stimuli = [Is_holding]);
 
 @named soma = Compartment(soma_dynamics)
 
@@ -33,9 +28,10 @@ dendrite_dynamics = HodgkinHuxley(Vₘ,
                              reversals[2:4],
                              geometry = Unitless(0.5),
                              capacitance = capacitance,
-                             stimuli = [dendrite_holding]);
+                             stimuli = [Id_holding]);
 
 @named dendrite = Compartment(dendrite_dynamics, extensions = [calcium_conversion])
+
 @named gc_soma = AxialConductance([Gate(SimpleGate, inv(p), name = :ps)], max_g = gc_val)
 @named gc_dendrite = AxialConductance([Gate(SimpleGate, inv(1-p), name = :pd)], max_g = gc_val)
 

demo/simplesynapse.jl

@@ -27,13 +27,11 @@ kdr_kinetics = [
 channels = [NaV, Kdr, leak];
 reversals = Equilibria([Na => 50.0mV, K => -77.0mV, Leak => -54.4mV])
 
-@named Iₑ = IonCurrent(NonIonic)
-@named I_hold = IonCurrent(NonIonic, 5000pA, dynamic = false)
-holding_current = Iₑ ~ I_hold
+@named Iₑ = Bias(5000pA)
 
 dynamics_1 = HodgkinHuxley(Vₘ, channels, reversals;
                            geometry = Cylinder(radius = 25µm, height = 400µm),
-                           stimuli = [holding_current]);
+                           stimuli = [Iₑ]);
 dynamics_2 = HodgkinHuxley(Vₘ, channels, reversals;
                            geometry = Cylinder(radius = 25µm, height = 400µm));
 

docs/src/basics.md

@@ -55,13 +55,11 @@ kdr_kinetics = [
 channels = [NaV, Kdr, leak];
 reversals = Equilibria([Na => 50.0mV, K => -77.0mV, Leak => -54.4mV])
 
-@named Iₑ = IonCurrent(NonIonic)
-@named I_hold = IonCurrent(NonIonic, 5000pA, dynamic = false)
-holding_current = Iₑ ~ I_hold
+@named Iₑ = Bias(5000pA)
 
 dynamics_1 = HodgkinHuxley(Vₘ, channels, reversals;
                          geometry = Cylinder(radius = 25µm, height = 400µm),
-                         stimuli = [holding_current])
+                         stimuli = [Iₑ])
 dynamics_2 = HodgkinHuxley(Vₘ, channels, reversals;
                            geometry = Cylinder(radius = 25µm, height = 400µm))
 
@@ -154,9 +152,8 @@ current and use the resulting symbolic to write an equation describing what the
 electrode current should be.
 
 ```@example gate_example; continued = true
-@named Iₑ = IonCurrent(NonIonic)
-@named I_hold = IonCurrent(NonIonic, 5000pA, dynamic = false)
-holding_current = Iₑ ~ I_hold
+@named Iₑ = Bias(5000pA)
+nothing # hide
 ```
 
 Finally, we construct our two neurons, providing the holding current stimulus to `neuron1`,
@@ -165,13 +162,14 @@ which will be our presynaptic neuron.
 ```@example gate_example
 dynamics_1 = HodgkinHuxley(Vₘ, channels, reversals;
                          geometry = Cylinder(radius = 25µm, height = 400µm),
-                         stimuli = [holding_current])
+                         stimuli = [Iₑ])
 
 dynamics_2 = HodgkinHuxley(Vₘ, channels, reversals;
                            geometry = Cylinder(radius = 25µm, height = 400µm))
 
 @named neuron1 = Compartment(dynamics_1)
 @named neuron2 = Compartment(dynamics_2)
+nothing # hide
 ``` 
 For our neurons to talk to each other, we'll need a model for a synaptic conductance. This
 time we'll use a model that's presented in a different form in the literature. A model of a
@@ -213,6 +211,7 @@ add_synapse!(topology, neuron1, neuron2, Glut)
 reversal_map = Dict([Glut => EGlut])
 
 @named net = NeuronalNetworkSystem(topology, reversal_map)
+nothing # hide
 ```
 Now we're ready to run our simulation.
 

docs/src/core/compartments.md

@@ -72,21 +72,17 @@ equation for an electrode current and provide it to `CompartmentSystem`:
 
 ```@example compartment_example
 import Unitful: µA, pA
-@named Iₑ = IonCurrent(NonIonic)
 
 # A 400 picoamp squarewave pulse when 100ms > t > 200ms
-@named I_rest = IonCurrent(NonIonic, 0.0µA, dynamic = false)
-@named I_step = IonCurrent(NonIonic, 400.0pA, dynamic = false)
-@parameters tstart = 100.0 [unit=ms] tstop = 200.0 [unit=ms]
-electrode_pulse = Iₑ ~ ifelse((t > tstart) & (t < tstop), I_step, I_rest)
+@named Iₑ = PulseTrain(amplitude = 400.0pA, duration = 100ms, delay = 100ms)
 
 stim_dynamics = HodgkinHuxley(Vₘ, channels, reversals;
                               geometry = soma_shape,
-                              stimuli = [electrode_pulse])
+                              stimuli = [Iₑ])
 
 @named neuron_stim = CompartmentSystem(stim_dynamics)
-@assert length.((equations(neuron_stim), states(neuron_stim), parameters(neuron_stim))) == (13,13,12); # hide
-neuron_stim # hide
+@assert length.((equations(neuron_stim), states(neuron_stim), parameters(neuron_stim))) == (13,13,8); # hide
+nothing # hide
 ```
 Putting it all together, our neuron simulation now produces a train of action potentials:
 

src/Conductor.jl

@@ -101,6 +101,7 @@ export output, get_output, timeconstant, steadystate, forward_rate, reverse_rate
 export Sphere, Cylinder, Point, Unitless, area, radius, height
 
 export HodgkinHuxley
+export Bias, PulseTrain
 
 # Metadata IDs
 struct ConductorMaxConductance end
@@ -111,6 +112,7 @@ abstract type AbstractNeuronalNetworkSystem <: AbstractTimeDependentSystem end
 
 include("util.jl")
 include("primitives.jl")
+include("stimulus.jl")
 include("gate.jl")
 include("conductance.jl")
 include("geometry.jl")

src/compartment.jl

@@ -18,7 +18,7 @@ struct HodgkinHuxley <: CompartmentForm
     "Axial (intercompartmental) conductances."
     axial_conductances::Vector{Tuple{AbstractConductanceSystem,Num}}
     "Experimental stimuli (for example, current injection)."
-    stimuli::Vector{Equation}
+    stimuli::Vector
 end
 
 # basic constructor
@@ -28,7 +28,7 @@ function HodgkinHuxley(
     channel_reversals;
     capacitance = 1µF/cm^2,
     geometry::Geometry = Point(),
-    stimuli::Vector{Equation} = Equation[])
+    stimuli::Vector = [])
 
     @parameters cₘ = ustrip(Float64, µF/cm^2, capacitance) [unit=µF/cm^2]
     synaptic_channels = Vector{AbstractConductanceSystem}()
@@ -123,21 +123,42 @@ function generate_currents!(gen, dynamics::HodgkinHuxley, Vₘ, aₘ)
     return currents
 end
 
-function process_stimuli!(currents, gen, dynamics::HodgkinHuxley)
+function process_stimulus!(currents, gen, stimulus::PulseTrain)
     (; eqs, dvs, ps, defs) = gen
-    for stimulus in dynamics.stimuli
-        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
-            validate(stimulus) && push!(eqs, stimulus)
-            push!(dvs, I)
-            push!(currents, -I)
-        end
+
+    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
+        validate(stimulus) && push!(eqs, stimulus)
+        push!(dvs, I)
+        push!(currents, -I)
+    end
+    return
+end
+
+# current bias stimulus
+function process_stimulus!(currents, gen, sym, stimulus::Bias{T}) where {T<:Current}
+    ps = gen.ps
+    push!(ps, sym)
+    push!(currents, -sym)
+    return
+end
+
+function process_stimulus!(currents, gen, sym, stimulus::PulseTrain{T}) where {T<:Current}
+    (; dvs, eqs) = gen
+    push!(dvs, sym)
+    push!(eqs, sym ~ current_pulses(t, stimulus))
+    push!(currents, -sym)
+end
+
+function process_stimuli!(currents, gen, dynamics::HodgkinHuxley)
+    for sym in dynamics.stimuli
+        process_stimulus!(currents, gen, sym, get_stimulus(sym))
     end
 end
 

src/stimulus.jl

@@ -1,18 +1,67 @@
 
-abstract type AbstractStimulus end
+abstract type Stimulus end
+get_stimulus(x::Num) = getmetadata(x, Stimulus)
 
-struct Stimulator{T<:AbstractStimulus}
-    channels::Vector{T} 
+# Fixed value
+struct Bias{T} <: Stimulus
+    val::T
 end
 
-struct PulseTrain{T<:Real} <: AbstractStimulus
-    interval::T
-    duration::T
-    delay::T
+function Bias(amplitude::T; name=Base.gensym("bias")) where {T<:Current}
+    val = ustrip(µA, amplitude)
+    sym = only(@parameters $name=val [unit=µA])
+    return setmetadata(sym, Stimulus, Bias{T}(amplitude))
+end
+
+struct PulseTrain{U} <: Stimulus
+    interval
+    duration
+    delay
     npulses::Int
-    offset::T
-    amplitude::T
+    offset
+    amplitude
+end
+
+function current_pulses(t, x)
+   (; interval, duration, delay, npulses, offset, amplitude) = x
+   if t > delay
+        telapsed = t - delay
+        pulse_number = (telapsed ÷ interval) + 1
+        pulse_number > npulses && return offset
+        t_in_period = telapsed - (pulse_number - 1)*interval
+        t_in_period < duration && return amplitude
+    end
+    return offset
+end
+
+@register_symbolic current_pulses(t, x::PulseTrain)
+
+function PulseTrain(; amplitude::T1, 
+                      duration::Time, 
+                      offset::T2 = 0.0µA,
+                      delay::Time, 
+                      npulses = 1,
+                      interval::Time = duration,
+                      name = Base.gensym("pulsetrain")) where {T1<:Current, T2<:Current}
+
+    val = ustrip(µA, offset)
+    sym = only(@variables $name(t)=val [unit=µA])
+    return setmetadata(sym, Stimulus, PulseTrain{T1}(ustrip(ms, interval),
+                                                             ustrip(ms, duration),
+                                                             ustrip(ms, delay),
+                                                             npulses,
+                                                             val, # offset
+                                                             ustrip(µA, amplitude)))
 end
 
-struct Waveform{T<:Real} <: AbstractStimulus end
+
+# Arbitrary input
+#struct CommandCurrent end
+#struct CommandPotential end
+#struct Stimulator{T<:Stimulus}
+#    channels::Vector{T} 
+#end
+#
+#struct Waveform{T} <: Stimulus end
+
 

test/hodgkinhuxley.jl

@@ -41,28 +41,25 @@ alphamss = ((4.0exp(-3.6111111111111107 - (0.05555555555555555Vₘ)) +
 channels = [NaV, Kdr, leak]
 @test [length(equations(x)) for x in channels] == [3,2,1]
 reversals = Equilibria([Na => 50.0mV, K => -77.0mV, Leak => -54.4mV])
-@named Iₑ = IonCurrent(NonIonic)
-@named I_rest = IonCurrent(NonIonic, 0.0µA, dynamic = false)
-@named I_step = IonCurrent(NonIonic, 400.0pA, dynamic = false)
-@parameters tstart = 100.0 [unit=ms] tstop = 200.0 [unit=ms]
-electrode_pulse = Iₑ ~ ifelse((t > tstart) & (t < tstop), I_step, I_rest)
+
+@named Iₑ = PulseTrain(amplitude = 400.0pA, duration = 100ms, delay = 100ms)
 
 dynamics = HodgkinHuxley(Vₘ, channels, reversals;
                          geometry = Sphere(radius = 20µm),
-                         stimuli = [electrode_pulse])
+                         stimuli = [Iₑ])
 
 @named neuron = Compartment(dynamics)
 
 @test length.([equations(neuron),
                states(neuron),
-               parameters(neuron)]) == [13,13,12]
+               parameters(neuron)]) == [13,13,8]
 
 time = 300.
 sim_sys = Simulation(neuron, time = time*ms, return_system = true)
 
 @test length.([equations(sim_sys),
                states(sim_sys),
-               parameters(sim_sys)]) == [4,4,12]
+               parameters(sim_sys)]) == [4,4,8]
 
 expect_syms = [:Vₘ, :NaV₊m, :NaV₊h, :Kdr₊n]
 @test all(x -> hasproperty(sim_sys, x), expect_syms)

test/simplesynapse.jl

@@ -31,13 +31,11 @@ kdr_kinetics = [
 channels = [NaV, Kdr, leak];
 reversals = Equilibria([Na => 50.0mV, K => -77.0mV, Leak => -54.4mV])
 
-@named Iₑ = IonCurrent(NonIonic)
-@named I_hold = IonCurrent(NonIonic, 5000pA, dynamic = false)
-holding_current = Iₑ ~ I_hold
+@named Iₑ = Bias(5000pA)
 
 geo = Cylinder(radius = 25µm, height = 400µm)
 
-dynamics_1 = HodgkinHuxley(Vₘ, channels, reversals; geometry = geo, stimuli = [holding_current]);
+dynamics_1 = HodgkinHuxley(Vₘ, channels, reversals; geometry = geo, stimuli = [Iₑ]);
 dynamics_2 = HodgkinHuxley(Vₘ, channels, reversals; geometry = geo);
 @named neuron1 = Compartment(dynamics_1)
 @named neuron2 = Compartment(dynamics_2)
@@ -62,7 +60,7 @@ reversal_map = Dict([Glut => EGlut])
 
 @test length.([equations(network),
                states(network),
-               parameters(network)]) == [29,29,19]
+               parameters(network)]) == [28,28,19]
 
 ttot = 250.
 simul_sys = Simulation(network, time = ttot*ms, return_system = true)