Population updates (#62)

Author: wsphillips

demo/pinsky_network.jl

@@ -1,5 +1,5 @@
 
-using LinearAlgebra, OrdinaryDiffEq, ModelingToolkit
+using LinearAlgebra, OrdinaryDiffEq, ModelingToolkit, Graphs
 LinearAlgebra.BLAS.set_num_threads(6)
 
 include(joinpath(@__DIR__, "traub_kinetics.jl"))
@@ -48,67 +48,63 @@ import Conductor: NMDA, AMPA
                                   [Gate(inv(1 + 0.28 * exp(-0.062(Vₘ - 60.0))); name = :e),
                                    Gate(S, [D(S) ~ -S/τNMDA]),
                                    Gate(inv(1 - p), name = :pnmda)],
-                                  max_s = 0.027mS)
+                                  max_s = 0.02mS)
 
 @named AMPAChan = SynapticChannel(ConstantValueEvent(S; threshold = 20mV), AMPA,
                                   [Gate(S, [D(S) ~ -S/τAMPA]),
                                    Gate(inv(1 - p), name = :pampa)],
-                                  max_s = 0.008mS)
+                                  max_s = 0.04mS)
 
 ENMDA = EquilibriumPotential(NMDA, 60mV, name = :NMDA)
 EAMPA = EquilibriumPotential(AMPA, 60mV, name = :AMPA)
 revmap = Dict([NMDAChan => ENMDA, AMPAChan => EAMPA])
 
-# A single neuron was "briefly" stimulated to trigger the network
+# Create a neuron group
+n_neurons = 100
+@named neurons = Population(mcneuron, n_neurons);
+
+# A brief stimulation to trigger the network
 @named I_pulse = PulseTrain(amplitude = 5.0µA / 0.5,
                             duration = 50ms,
                             delay = 150.0ms,
                             offset = -0.5µA / 0.5)
 
-soma_stimulated = Compartment(Vₘ, deepcopy(soma_dynamics);
-                              geometry = Unitless(0.5), # FIXME: should take p param
-                              capacitance = capacitance,
-                              stimuli = [I_pulse],
-                              name = :soma)
-mcstim_topology = MultiCompartmentTopology([soma_stimulated, dendrite]);
-add_junction!(mcstim_topology, soma_stimulated, dendrite, (gc_soma, gc_dendrite))
-@named mcneuron_stim = MultiCompartment(mcstim_topology)
-
-# Need to introduce 10% gca variance as per Pinsky/Rinzel
-n_neurons = 100
-neuronpopulation = [Conductor.replicate(mcneuron) for _ in 1:n_neurons];
-neuronpopulation[4] = mcneuron_stim
-topology = NetworkTopology(neuronpopulation, [NMDAChan, AMPAChan]);
+# apply the stimulus to the root compartment of the 4th neuron
+add_stimuli!(neurons, [I_pulse], (4,1))
+topology = NetworkTopology(neurons, [NMDAChan, AMPAChan]);
 
-using Graphs
+# Generate connectivity graphs
 nmda_g = random_regular_digraph(n_neurons, fld(n_neurons, 5), dir = :in)
 ampa_g = random_regular_digraph(n_neurons, fld(n_neurons, 5), dir = :in)
 
-# We could allow users to supply a lambda/function to map in order to get this behavior
+# Use graphs to create synapses from soma (output) to dendrites (inputs) for neurons
 for (i, e) in enumerate(edges(nmda_g))
-    add_synapse!(topology, neuronpopulation[src(e)].soma, neuronpopulation[dst(e)].dendrite,
+    add_synapse!(topology, src(e), (dst(e),2),
                  NMDAChan, 1.0)
 end
 
 for (i, e) in enumerate(edges(ampa_g))
-    add_synapse!(topology, neuronpopulation[src(e)].soma, neuronpopulation[dst(e)].dendrite,
+    add_synapse!(topology, src(e), (dst(e),2),
                  AMPAChan, 1.0)
 end
 
 @named net = NeuronalNetworkSystem(topology, revmap);
 t_total = 2000.0
-simp = Simulation(net, t_total*ms, return_system = true)
 prob = Simulation(net, t_total*ms)
-sol = solve(prob, RK4(); adaptive = false, dt = 0.05);
+@time sol = solve(prob, RK4(); adaptive = false, dt = 0.05);
 
+############################################################################################
 # Pinsky and Rinzel displayed their results as a plot of N neurons over 20mV
+using Statistics, Plots
+
+# helper function
 indexof(sym, syms) = findfirst(isequal(sym), syms)
+
+simp = ODESystem(net)
 dvs = states(simp)
 interpolated = sol(1:0.2:t_total,
-                   idxs = [indexof(x.soma.Vₘ, dvs) for x in neuronpopulation[5:end]])
+                   idxs = [indexof(x.soma.Vₘ, dvs) for x in neurons[5:end]])
 abovethold = reduce(hcat, interpolated.u) .> 10.0
-using Statistics
 final = mean(abovethold, dims = 1)'
-using Plots
 plot(final)
 

demo/simplesynapse.jl

@@ -1,5 +1,5 @@
 # Example of synapses
-using Conductor, OrdinaryDiffEq, Unitful, ModelingToolkit, Plots
+using Conductor, Unitful, ModelingToolkit, OrdinaryDiffEq, Plots
 import Unitful: mV, mS, cm, µm, pA, nA, mA, µA, ms, nS, pS
 import Conductor: Na, K
 
@@ -28,18 +28,15 @@ kdr_kinetics = [Gate(AlphaBeta,
 @named Kdr = IonChannel(Potassium, kdr_kinetics, max_g = 36mS / cm^2);
 @named leak = IonChannel(Leak, max_g = 0.3mS / cm^2);
 
-channels = [NaV, Kdr, leak];
 reversals = Equilibria([Na => 50.0mV, K => -77.0mV, Leak => -54.4mV]);
-dynamics = HodgkinHuxley(channels, reversals);
+dynamics = HodgkinHuxley([NaV, Kdr, leak], reversals);
+geometry = Cylinder(radius = 25µm, height = 400µm)
 
 # Stimulate first neuron with a holding current to drive APs
 @named holding_current = Bias(5000pA);
-@named neuron1 = Compartment(Vₘ, dynamics;
-                             geometry = Cylinder(radius = 25µm, height = 400µm),
-                             stimuli = [holding_current]);
+@named neurons = Population(Compartment(Vₘ, dynamics; geometry), 2);
 
-@named neuron2 = Compartment(Vₘ, dynamics;
-                             geometry = Cylinder(radius = 25µm, height = 400µm));
+add_stimuli!(neurons, [holding_current], 1)
 
 ############################################################################################
 # Event-based Synaptic model
@@ -52,9 +49,9 @@ syn_kinetics = Gate(SimpleGate, m, [D(m) ~ -m/τsyn])
 event_model = ConstantValueEvent(m) # increment `m` gate by 1 for each incoming AP
 @named ExpAMPA = SynapticChannel(event_model, Cation, [syn_kinetics]; max_s = 30nS)
 
-top = NetworkTopology([neuron1, neuron2], [ExpAMPA]);
+top = NetworkTopology(neurons, [ExpAMPA]);
 
-top[neuron1, neuron2] = ExpAMPA(10nS)
+top[neurons[1], neurons[2]] = ExpAMPA(10nS)
 rev_map = Dict([ExpAMPA => EGlut])
 
 @named net = NeuronalNetworkSystem(top, rev_map)
@@ -63,8 +60,8 @@ total_time = 250.0
 sim = Simulation(net, total_time * ms)
 sol = solve(sim, Rosenbrock23(), abstol=1e-3, reltol=1e-3);
 
-plot(plot(sol, idxs = [neuron1.Vₘ]),
-     plot(sol, idxs = [neuron2.Vₘ]),
+plot(plot(sol, idxs = [neurons[1].Vₘ]),
+     plot(sol, idxs = [neurons[2].Vₘ]),
      layout=(2,1))
 
 ############################################################################################
@@ -78,18 +75,18 @@ syn_kinetics2 = Gate(SteadyStateTau, syn∞, tausyn, name = :z)
 
 @named IntAMPA = SynapticChannel(IntegratedSynapse(), Cation, [syn_kinetics2]; max_s = 30nS)
 
-top2 = NetworkTopology([neuron1, neuron2], [IntAMPA]);
+top2 = NetworkTopology(neurons, [IntAMPA]);
 
-top2[neuron1, neuron2] = IntAMPA(30nS)
+top2[neurons[1], neurons[2]] = IntAMPA(40nS)
 rev_map2 = Dict([IntAMPA => EGlut])
 
 @named net2 = NeuronalNetworkSystem(top2, rev_map2)
 
 total_time = 250.0
-sim2 = Simulation(net2, total_time * ms)
+sim2 = Simulation(net2, (0.0ms, total_time * ms))
 sol2 = solve(sim2, Rosenbrock23(), abstol = 1e-3, reltol = 1e-3);
 
-plot(plot(sol2, idxs = [neuron1.Vₘ]),
-     plot(sol2, idxs = [neuron2.Vₘ]),
+plot(plot(sol2, idxs = [neurons[1].Vₘ]),
+     plot(sol2, idxs = [neurons[2].Vₘ]),
      layout=(2,1))
 

src/Conductor.jl

@@ -107,7 +107,7 @@ export output, get_output, timeconstant, steadystate, forward_rate, reverse_rate
 export Sphere, Cylinder, Point, Unitless, area, radius, height
 
 export HodgkinHuxley
-export Bias, PulseTrain
+export Bias, PulseTrain, add_stimuli!
 
 # Metadata IDs
 struct ConductorMaxConductance end
@@ -137,6 +137,7 @@ struct IntegratedSynapse <: SynapticModel end
 
 export ConstantValueEvent, IntegratedSynapse
 
+include("populations.jl")
 include("util.jl")
 include("primitives.jl")
 include("stimulus.jl")
@@ -148,6 +149,5 @@ include("compartment.jl")
 include("multicompartment.jl")
 include("network.jl")
 include("simulation.jl")
-include("populations.jl")
 include("callbacks.jl")
 end # module

src/multicompartment.jl

@@ -186,3 +186,17 @@ function Base.convert(::Type{ODESystem}, mcsys::MultiCompartmentSystem)
                        systems = systems, checks = CheckComponents)
     return odesys
 end
+
+function add_stimuli!(pop::Population{MultiCompartmentSystem},
+                      stimuli::Vector{<:StimulusModel}, idxs)
+    nrn = pop.neurons[idxs[1]]
+    subcomps = get_compartments(nrn)
+    subcomps[idxs[2]] = remake(subcomps[idxs[2]]; stimuli = stimuli)
+    top = get_topology(nrn)
+    @set! top.compartments = subcomps
+    nrn = remake(nrn, topology = top)
+    pop.neurons[idxs[1]] = nrn
+    return nothing
+end
+
+

src/network.jl

@@ -2,21 +2,27 @@
 struct NetworkTopology
     multigraph::Dict{Any, SparseMatrixCSC{Float64, Int64}}
     neurons::Vector
+    lens::Vector{Int64}
+    function NetworkTopology(multigraph::Dict, neurons::Vector, lens::Vector{Int64})
+        new(multigraph, neurons, lens)
+    end
 end
 
 synaptic_systems(topology::NetworkTopology) = keys(graph(topology))
 
-function NetworkTopology(neurons::Vector, synaptic_systems::Vector)
+function NetworkTopology(neurons, synaptic_systems::Vector)
     m = length(neurons)
-    n = sum(length(neuron) for neuron in neurons)
+    lens = [length(neuron) for neuron in neurons]
+    n = sum(lens)
     multigraph = Dict(x => sparse(Int64[], Int64[], Float64[], n, m) for x in synaptic_systems)
-    return NetworkTopology(multigraph, neurons)
+    return NetworkTopology(multigraph, collect(neurons), lens)
 end
 
 function NetworkTopology(g::SimpleDiGraph, neurons, synaptic_system,
                          default_weight = get_gbar(synaptic_system))
     multigraph = Dict(synaptic_system => adjacency_matrix(g) * default_weight)
-    return NetworkTopology(multigraph, neurons)
+    lens = [length(neuron) for neuron in neurons]
+    return NetworkTopology(multigraph, neurons, lens)
 end
 
 neurons(topology::NetworkTopology) = getfield(topology, :neurons)
@@ -32,29 +38,25 @@ graph(topology::NetworkTopology) = getfield(topology, :multigraph)
 find_source(pre, topology) = findfirst(isequal(first(pre)), root_compartments(topology))
 find_target(post, topology) = findfirst(isequal(post), compartments(topology))
 
-function add_synapse!(topology, pre, post, synaptic_system, weight)
+function add_synapse!(topology, pre::AbstractCompartmentSystem,
+                      post::AbstractCompartmentSystem, synaptic_system, weight)
     src = find_source(pre, topology)
     dst = find_target(post, topology)
     g = graph(topology)[synaptic_system]
     g[dst, src] = weight # underlying matrix is transposed
 end
 
-# FIXME: need updates
-#=
-function remove_synapse!(topology, pre, post, synaptic_system)
-    src = find_source(pre, topology)
-    dst = find_target(post, topology)
+function add_synapse!(topology, pre::Integer, post::Integer, synaptic_system, weight)
     g = graph(topology)[synaptic_system]
-    g[dst, src] = zero(Num)
-    dropzeros!(g)
+    g[post, pre] = weight # underlying matrix is transposed
 end
 
-function add_layer!(topology, synaptic_system::ConductanceSystem,
-                    g = SimpleDiGraph(length(vertices(topology))),
-                    default_weight = get_gbar(synaptic_system))
-    push!(graph(topology), synaptic_system => adjacency_matrix(g) * default_weight)
+function add_synapse!(topology, pre::Integer, post::Tuple{Integer,Integer},
+                      synaptic_system, weight)
+    g = graph(topology)[synaptic_system]
+    post_idx = sum(topology.lens[1:post[1]-1]) + post[2]
+    g[post_idx, pre] = weight # underlying matrix is transposed
 end
-=#
 
 Base.eltype(nt::NetworkTopology) = AbstractCompartmentSystem
 Base.length(nt::NetworkTopology) = length(neurons(nt))
@@ -146,7 +148,7 @@ synaptic_systems(sys::NeuronalNetworkSystem) = synaptic_systems(get_topology(sys
 compartments(sys::NeuronalNetworkSystem) = compartments(get_topology(sys))
 
 function connect_synapses!(gen, syn_model, comps, topology, reversal_map)
-        # this method assumes weights scale the size of the event/alpha)
+    # this method assumes weights scale the size of the event/alpha)
     new_compartments = deepcopy(comps)
     reversal = reversal_map[syn_model]
     for (i,comp) in enumerate(new_compartments)
@@ -226,7 +228,6 @@ function NeuronalNetworkSystem(topology::NetworkTopology, reversal_map,
     for sys in synaptic_systems(topology)
         comps = connect_synapses!(gen, sys, comps, topology, reversal_map)
     end
-
     for neuron in neurons(topology)
         if typeof(neuron) <: MultiCompartmentSystem
             mctop = get_topology(neuron)