From dd578567654d9cf4a067a85cf00a21364f41c7a9 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Wed, 26 Feb 2025 10:55:54 -0500 Subject: [PATCH 01/18] Some temporary commits --- examples/adams_example_for_brain_r01.jl | 2 +- examples/qif_ngnmm_learning_demo.jl | 112 +++++++++++ .../GraphDynamicsInterop.jl | 5 +- .../connection_interop.jl | 35 ++++ src/GraphDynamicsInterop/neuron_interop.jl | 29 +++ src/Neuroblox.jl | 2 + src/blox/adam_examples.jl | 174 ++++++++++++++++++ src/blox/adam_neurons.jl | 105 +++++++++++ src/blox/connections.jl | 35 ++++ 9 files changed, 497 insertions(+), 2 deletions(-) create mode 100644 examples/qif_ngnmm_learning_demo.jl create mode 100644 src/blox/adam_examples.jl create mode 100644 src/blox/adam_neurons.jl diff --git a/examples/adams_example_for_brain_r01.jl b/examples/adams_example_for_brain_r01.jl index 29158a7b..e2c9ee37 100644 --- a/examples/adams_example_for_brain_r01.jl +++ b/examples/adams_example_for_brain_r01.jl @@ -19,7 +19,7 @@ add_edge!(g, exci_PING => inhi_PING; weight=10.0) add_edge!(g, inhi_PING => exci_PING; weight=10.0) @named sys = system_from_graph(g) -sys = structural_simplify(sys) +#sys = structural_simplify(sys) sim_dur = 1000.0 prob = SDEProblem(sys, [], (0.0, sim_dur)) diff --git a/examples/qif_ngnmm_learning_demo.jl b/examples/qif_ngnmm_learning_demo.jl new file mode 100644 index 00000000..23d24b36 --- /dev/null +++ b/examples/qif_ngnmm_learning_demo.jl @@ -0,0 +1,112 @@ +using DifferentialEquations +using Distributions +using Statistics +using Random +using Plots + +function qif_ngnmm_params(;Δₑ=1.0, + τₘₑ=20.0, + Hₑ=1.3, + Jₑₑ=8.0, + Jₑᵢ=10.0, + Δᵢ=1.0, + τₘᵢ=10.0, + Hᵢ=-5.0, + Jᵢₑ=10.0, + Jᵢᵢ=0.0, + w₁¹=1.0, + w₂¹=5.0, + w₁²=1.0, + w₂²=5.0, + curr_stim=0.0) + + return [Δₑ, τₘₑ, Hₑ, Jₑₑ, Jₑᵢ, Δᵢ, τₘᵢ, Hᵢ, Jᵢₑ, Jᵢᵢ, w₁¹, w₂¹, w₁², w₂², curr_stim] +end + +p = qif_ngnmm_params() + +function simple_cs_learning!(du, u, p, t) + Δₑ, τₘₑ, Hₑ, Jₑₑ, Jₑᵢ, Δᵢ, τₘᵢ, Hᵢ, Jᵢₑ, Jᵢᵢ, w₁¹, w₂¹, w₁², w₂², curr_stim = p + + rₑ¹, Vₑ¹, rᵢ¹, Vᵢ¹, rₑ², Vₑ², rᵢ², Vᵢ², I₁, I₂ = u + + du[1] = Δₑ/(π*τₘₑ^2) + 2*rₑ¹*Vₑ¹/τₘₑ + du[2] = (Vₑ¹^2 + Hₑ + 4*sin(5*2*π*t/1000))/τₘₑ - τₘₑ*(π*rₑ¹)^2 + Jₑₑ*rₑ¹ + Jᵢₑ*rᵢ¹ + w₁¹*I₁ + w₂¹*I₂ - 1*rₑ² + du[3] = Δᵢ/(π*τₘᵢ^2) + 2*rᵢ¹*Vᵢ¹/τₘᵢ + du[4] = (Vᵢ¹^2 + Hᵢ + 2*sin(5*2*π*t/1000))/τₘᵢ - τₘᵢ*(π*rᵢ¹)^2 + Jₑᵢ*rₑ¹ + Jᵢᵢ*rᵢ¹ + du[5] = Δₑ/(π*τₘₑ^2) + 2*rₑ²*Vₑ²/τₘₑ + du[6] = (Vₑ²^2 + Hₑ + 4*sin(5*2*π*(t-100)/1000))/τₘₑ - τₘₑ*(π*rₑ²)^2 + Jₑₑ*rₑ² + Jᵢₑ*rᵢ² + w₁²*I₁ + w₂²*I₂ - 1*rₑ¹ + du[7] = Δᵢ/(π*τₘᵢ^2) + 2*rᵢ²*Vᵢ²/τₘᵢ + du[8] = (Vᵢ²^2 + Hᵢ + 2*sin(5*2*π*(t-100)/1000))/τₘᵢ - τₘᵢ*(π*rᵢ²)^2 + Jₑᵢ*rₑ² + Jᵢᵢ*rᵢ² + du[9] = -u[9]/30 + du[10] = -u[10]/30 +end + + +max_time = 100000.0 +stim1_times = collect(250:500:max_time) +stim2_times = collect(500:500:max_time) +all_stim_times = sort(vcat(stim1_times, stim2_times)) + +condtion_s1(u, t, integrator) = t ∈ stim1_times +condtion_s2(u, t, integrator) = t ∈ stim2_times + +function affect_s1!(integrator) + integrator.u[9] += 30.0 + integrator.p[15] = 1.0 +end + +function affect_s2!(integrator) + integrator.u[10] += 30.0 + integrator.p[15] = 2.0 +end + +cb_s1 = DiscreteCallback(condtion_s1, affect_s1!) +cb_s2 = DiscreteCallback(condtion_s2, affect_s2!) + +eval_times = all_stim_times[2:end] +eval_times .-= 50.0 + +condition_eval(u, t, integrator) = t ∈ eval_times + +all_choices = [] + +function learn!(integrator) + u = integrator.u + p = integrator.p + + str1 = u[2] + str2 = u[6] + + choice = str1 > str2 + + if p[15] == 1 && choice + p[11] += 0.1*rand() + p[11] = min(p[11], 5.0) + push!(all_choices, 1.0) + elseif p[15] == 1 && !choice + p[13] -= 0.1*rand() + p[13] = max(p[13], -5.0) + push!(all_choices, 0.0) + elseif p[15] == 2 && !choice + p[14] += 0.1*rand() + p[14] = min(p[14], 5.0) + push!(all_choices, 1.0) + elseif p[15] == 2 && choice + p[12] -= 0.1*rand() + p[12] = max(p[12], -5.0) + push!(all_choices, 0.0) + end +end + +cb_learn = DiscreteCallback(condition_eval, learn!) + +tstop = sort(vcat(all_stim_times, eval_times)) + +u₀ = zeros(10) + +tspan = (0.0, max_time) +cbs = CallbackSet(cb_s1, cb_s2, cb_learn) +prob = ODEProblem(simple_cs_learning!, u₀, tspan, p) +sol = solve(prob, Tsit5(), callback = cbs, tstops = tstop, saveat=1.0) +plot(all_choices) \ No newline at end of file diff --git a/src/GraphDynamicsInterop/GraphDynamicsInterop.jl b/src/GraphDynamicsInterop/GraphDynamicsInterop.jl index 6b945f5c..e590bd96 100644 --- a/src/GraphDynamicsInterop/GraphDynamicsInterop.jl +++ b/src/GraphDynamicsInterop/GraphDynamicsInterop.jl @@ -51,7 +51,10 @@ using ..Neuroblox: PINGNeuronInhib, AbstractPINGNeuron, Connector, - VanDerPol + VanDerPol, + AdamPYR, + AdamINP, + AbstractAdamNeuron using GraphDynamics: GraphDynamics, diff --git a/src/GraphDynamicsInterop/connection_interop.jl b/src/GraphDynamicsInterop/connection_interop.jl index 86d364f3..508955b5 100644 --- a/src/GraphDynamicsInterop/connection_interop.jl +++ b/src/GraphDynamicsInterop/connection_interop.jl @@ -1047,3 +1047,38 @@ function (c::PINGConnection)(blox_src::Subsystem{PINGNeuronInhib}, blox_dst::Sub (; jcn = w * s * (V_I - V)) end +# #------------------------- +# Adam Network +# #------------------------- +# PING Network +struct AdamConnection <: ConnectionRule + w::Float64 + V_E::Float64 + V_I::Float64 +end +Base.zero(::Type{AdamConnection}) = AdamConnection(0.0, 0.0, 0.0) + +function get_connection(blox_src::AdamPYR, blox_dst::AbstractAdamNeuron, kwargs) + (;w_val, name) = generate_weight_param(blox_src, blox_dst, kwargs) + V_E = get(kwargs, :V_E, 0.0) + (; conn = AdamConnection(w_val, V_E, 0.0), names=[name]) +end +function get_connection(blox_src::AdamINP, blox_dst::AbstractAdamNeuron, kwargs) + (;w_val, name) = generate_weight_param(blox_src, blox_dst, kwargs) + V_I = get(kwargs, :V_I, 0.0) + (; conn = AdamConnection(w_val, V_I, -80.0), names=[name]) +end + +function (c::AdamConnection)(blox_src::Subsystem{AdamPYR}, blox_dst::Subsystem{<:AbstractAdamNeuron}) + (; w, V_E) = c + (;sₐₘₚₐ) = blox_src + (;V) = blox_dst + (; jcn = w * sₐₘₚₐ * (V - V_E)) +end + +function (c::AdamConnection)(blox_src::Subsystem{AdamINP}, blox_dst::Subsystem{<:AbstractAdamNeuron}) + (; w, V_I) = c + (;sᵧ) = blox_src + (;V) = blox_dst + (; jcn = w * sᵧ * (V - V_I)) +end diff --git a/src/GraphDynamicsInterop/neuron_interop.jl b/src/GraphDynamicsInterop/neuron_interop.jl index 4fedcf6c..55d5ec98 100644 --- a/src/GraphDynamicsInterop/neuron_interop.jl +++ b/src/GraphDynamicsInterop/neuron_interop.jl @@ -24,6 +24,35 @@ function define_neuron(sys; mod=@__MODULE__()) system = structural_simplify(sys.system; fully_determined=false) params = get_ps(system) t = Symbol(get_iv(system)) +function define_neurons() + for (name, T) ∈ [(:hhne, HHNeuronExciBlox) + (:hhni, HHNeuronInhibBlox) + (:hhni_msn_adam, HHNeuronInhib_MSN_Adam_Blox) + (:hhni_fsi_adam, HHNeuronInhib_FSI_Adam_Blox) + (:hhne_stn_adam, HHNeuronExci_STN_Adam_Blox) + (:hhni_GPe_adam, HHNeuronInhib_GPe_Adam_Blox) + (:ngei, NextGenerationEIBlox) + (:wc, WilsonCowan) + (:ho, HarmonicOscillator) + (:jr, JansenRit) # Note! Regular JansenRit can support delays, and I have not yet implemented this! + (:if, IFNeuron) + (:lif, LIFNeuron) + (:qif, QIFNeuron) + (:izh, IzhikevichNeuron) + (:lif_exci, LIFExciNeuron) + (:lif_inh, LIFInhNeuron) + (:pexci, PINGNeuronExci) + (:pinhib, PINGNeuronInhib) + (:VdP, VanDerPol{NonNoisy}) + (:VdPN, VanDerPol{Noisy}) + (:ko, KuramotoOscillator{NonNoisy}) + (:kon, KuramotoOscillator{Noisy}) + (:adpyr, AdamPYR) + (:adinp, AdamINP)] + sys = T(;name) + system = structural_simplify(sys.system; fully_determined=false) + params = get_ps(system) + t = Symbol(get_iv(system)) states = [s for s ∈ unknowns(system) if !MTK.isinput(s)] inputs = [s for s ∈ unknowns(system) if MTK.isinput(s)] diff --git a/src/Neuroblox.jl b/src/Neuroblox.jl index 3979495d..598584a1 100644 --- a/src/Neuroblox.jl +++ b/src/Neuroblox.jl @@ -113,6 +113,7 @@ include("blox/subcortical_blox.jl") include("blox/stochastic.jl") include("blox/discrete.jl") include("blox/ping_neuron_examples.jl") +include("blox/adam_neurons.jl") include("blox/reinforcement_learning.jl") include("gui/GUI.jl") include("blox/connections.jl") @@ -260,4 +261,5 @@ export voltage_timeseries, meanfield_timeseries, state_timeseries, get_neurons, export AdjacencyMatrix, Connector, connection_rule, connection_equations, connection_spike_affects, connection_learning_rules, connection_callbacks export inputs, outputs, equations, unknowns, parameters, discrete_events export MetabolicHHNeuron +export AdamPYR, AdamINP end diff --git a/src/blox/adam_examples.jl b/src/blox/adam_examples.jl new file mode 100644 index 00000000..13f15c5b --- /dev/null +++ b/src/blox/adam_examples.jl @@ -0,0 +1,174 @@ +using Neuroblox +using OrdinaryDiffEq +using Random, Distributions +using Plots +import Neuroblox: AbstractNeuronBlox, paramscoping +using BenchmarkTools + +abstract type AbstractAdamNeuron <: AbstractNeuronBlox end + +struct AdamPYR <: AbstractAdamNeuron + params + system + namespace + + function AdamPYR(;name, + namespace=nothing, + C=1.0, + Eₙₐ=50, + ḡₙₐ=100, + Eₖ=-100, + ḡₖ=80, + Eₗ=-67, + ḡₗ=0.05, + τₑ=1.5, + Iₐₚₚ=-0.25, + Iₙₒᵢₛₑ=0.0) + p = paramscoping(C=C, Eₙₐ=Eₙₐ, ḡₙₐ=ḡₙₐ, Eₖ=Eₖ, ḡₖ=ḡₖ, Eₗ=Eₗ, ḡₗ=ḡₗ, Iₐₚₚ=Iₐₚₚ, Iₙₒᵢₛₑ=Iₙₒᵢₛₑ) + C, Eₙₐ, ḡₙₐ, Eₖ, ḡₖ, Eₗ, ḡₗ, Iₐₚₚ, Iₙₒᵢₛₑ = p + sts = @variables V(t)=0.0 m(t)=0.0 h(t)=0.0 n(t)=0.0 sₐₘₚₐ(t)=0.0 [output=true] jcn(t) [input=true] + + αₘ(v) = 0.32*(v+54.0)/(1.0 - exp(-(v+54.0)/4.0)) + βₘ(v) = 0.28*(v+27.0)/(exp((v+27.0)/5.0) - 1.0) + αₕ(v) = 0.128*exp((v+50.0)/18.0) + βₕ(v) = 4.0/(1.0 + exp(-(v+27.0)/5.0)) + αₙ(v) = 0.032*(v+52.0)/(1.0 - exp(-(v+52.0)/5.0)) + βₙ(v) = 0.5*exp(-(v+57.0)/40.0) + + m∞(v) = αₘ(v)/(αₘ(v) + βₘ(v)) + h∞(v) = αₕ(v)/(αₕ(v) + βₕ(v)) + n∞(v) = αₙ(v)/(αₙ(v) + βₙ(v)) + + τₘ(v) = 1.0/(αₘ(v) + βₘ(v)) + τₕ(v) = 1.0/(αₕ(v) + βₕ(v)) + τₙ(v) = 1.0/(αₙ(v) + βₙ(v)) + + gₐₘₚₐ(v) = 5*(1+tanh(v/4)) + + eqs = [D(V) ~ (Iₐₚₚ + Iₙₒᵢₛₑ - ḡₙₐ*m^3*h*(V - Eₙₐ) - ḡₖ*n^4*(V - Eₖ) - ḡₗ*(V - Eₗ) - jcn)/C, + D(m) ~ (m∞(V) - m)/τₘ(V), + D(h) ~ (h∞(V) - h)/τₕ(V), + D(n) ~ (n∞(V) - n)/τₙ(V), + D(sₐₘₚₐ) ~ gₐₘₚₐ(V)*(1-sₐₘₚₐ) - sₐₘₚₐ/τₑ + ] + + sys = ODESystem(eqs, t, sts, p; name=name) + + new(p, sys, namespace) + + end +end + +struct AdamINP <: AbstractAdamNeuron + params + system + namespace + + function AdamINP(;name, + namespace=nothing, + C=1.0, + Eₙₐ=50, + ḡₙₐ=100, + Eₖ=-100, + ḡₖ=80, + Eₗ=-67, + ḡₗ=0.05, + τᵢ=6, + Iₐₚₚ=0.1, + Iₙₒᵢₛₑ=0.0) + p = paramscoping(C=C, Eₙₐ=Eₙₐ, ḡₙₐ=ḡₙₐ, Eₖ=Eₖ, ḡₖ=ḡₖ, Eₗ=Eₗ, ḡₗ=ḡₗ, Iₐₚₚ=Iₐₚₚ, Iₙₒᵢₛₑ=Iₙₒᵢₛₑ) + C, Eₙₐ, ḡₙₐ, Eₖ, ḡₖ, Eₗ, ḡₗ, Iₐₚₚ, Iₙₒᵢₛₑ = p + sts = @variables V(t)=0.0 m(t)=0.0 h(t)=0.0 n(t)=0.0 sᵧ(t)=0.0 [output=true] jcn(t) [input=true] + + αₘ(v) = 0.32*(v+54.0)/(1.0 - exp(-(v+54.0)/4.0)) + βₘ(v) = 0.28*(v+27.0)/(exp((v+27.0)/5.0) - 1.0) + αₕ(v) = 0.128*exp((v+50.0)/18.0) + βₕ(v) = 4.0/(1.0 + exp(-(v+27.0)/5.0)) + αₙ(v) = 0.032*(v+52.0)/(1.0 - exp(-(v+52.0)/5.0)) + βₙ(v) = 0.5*exp(-(v+57.0)/40.0) + + m∞(v) = αₘ(v)/(αₘ(v) + βₘ(v)) + h∞(v) = αₕ(v)/(αₕ(v) + βₕ(v)) + n∞(v) = αₙ(v)/(αₙ(v) + βₙ(v)) + + τₘ(v) = 1.0/(αₘ(v) + βₘ(v)) + τₕ(v) = 1.0/(αₕ(v) + βₕ(v)) + τₙ(v) = 1.0/(αₙ(v) + βₙ(v)) + + gᵧ(v) = 2*(1+tanh(v/4)) + + eqs = [D(V) ~ (Iₐₚₚ + Iₙₒᵢₛₑ - ḡₙₐ*m^3*h*(V - Eₙₐ) - ḡₖ*n^4*(V - Eₖ) - ḡₗ*(V - Eₗ) - jcn)/C, + D(m) ~ (m∞(V) - m)/τₘ(V), + D(h) ~ (h∞(V) - h)/τₕ(V), + D(n) ~ (n∞(V) - n)/τₙ(V), + D(sᵧ) ~ gᵧ(V)*(1-sᵧ) - sᵧ/τᵢ + ] + + sys = ODESystem(eqs, t, sts, p; name=name) + + new(p, sys, namespace) + + end +end + +function Connector( + blox_src::AdamPYR, + blox_dest::AbstractAdamNeuron; + kwargs... +) + + sys_pre = blox_src.system + sys_post = blox_dest.system + w = generate_weight_param(blox_src, blox_dest; kwargs...) + E_Exc = haskey(kwargs, :E_Exc) ? kwargs[:E_Exc] : 0.0 + s = only(outputs(sys_pre; namespaced=true)) + + eq = sys_post.jcn ~ w*s*(sys_post.V - E_Exc) + + return Connector(nameof(sys_pre), nameof(sys_post); equation=eq, weight=w) + +end + +function Connector( + blox_src::AdamINP, + blox_dest::AbstractAdamNeuron; + kwargs... +) + + sys_pre = blox_src.system + sys_post = blox_dest.system + w = generate_weight_param(blox_src, blox_dest; kwargs...) + E_Inh = haskey(kwargs, :E_Inh) ? kwargs[:E_Inh] : -80.0 + s = only(outputs(sys_pre; namespaced=true)) + + eq = sys_post.jcn ~ w*s*(sys_post.V - E_Inh) + + return Connector(nameof(sys_pre), nameof(sys_post); equation=eq, weight=w) +end + +ḡᵢ = 0.5 +ḡₑ = 0.2 + +NI = 20 +NE = 80 + +exci = [AdamPYR(name=Symbol("PYR$i"), Iₐₚₚ=rand(Normal(0.25, 0.05))) for i in 1:NE] +inhi = [AdamINP(name=Symbol("INP$i"), Iₐₚₚ=rand(Normal(0.3, 0.05))) for i in 1:NI] # bump up to 0.3 + +g = MetaDiGraph() + +for ne ∈ exci + for ni ∈ inhi + add_edge!(g, ne => ni; weight=ḡₑ/NE) + add_edge!(g, ni => ne; weight=ḡᵢ/NI) + end +end + +tspan = (0.0, 2000.0) +begin + @btime @named sys = system_from_graph(g, graphdynamics=true) + @btime prob = ODEProblem(sys, [], tspan) + @btime sol = solve(prob, Tsit5(), saveat=0.5) +end + +plot(sol, idxs=1:5:(NE+NI)*5) \ No newline at end of file diff --git a/src/blox/adam_neurons.jl b/src/blox/adam_neurons.jl new file mode 100644 index 00000000..303c62a9 --- /dev/null +++ b/src/blox/adam_neurons.jl @@ -0,0 +1,105 @@ +abstract type AbstractAdamNeuron <: AbstractNeuronBlox end + +struct AdamPYR <: AbstractAdamNeuron + params + system + namespace + + function AdamPYR(;name, + namespace=nothing, + C=1.0, + Eₙₐ=50, + ḡₙₐ=100, + Eₖ=-100, + ḡₖ=80, + Eₗ=-67, + ḡₗ=0.05, + τₑ=1.5, + Iₐₚₚ=-0.25, + Iₙₒᵢₛₑ=0.0) + p = paramscoping(C=C, Eₙₐ=Eₙₐ, ḡₙₐ=ḡₙₐ, Eₖ=Eₖ, ḡₖ=ḡₖ, Eₗ=Eₗ, ḡₗ=ḡₗ, Iₐₚₚ=Iₐₚₚ, Iₙₒᵢₛₑ=Iₙₒᵢₛₑ) + C, Eₙₐ, ḡₙₐ, Eₖ, ḡₖ, Eₗ, ḡₗ, Iₐₚₚ, Iₙₒᵢₛₑ = p + sts = @variables V(t)=0.0 m(t)=0.0 h(t)=0.0 n(t)=0.0 sₐₘₚₐ(t)=0.0 [output=true] jcn(t) [input=true] + + αₘ(v) = 0.32*(v+54.0)/(1.0 - exp(-(v+54.0)/4.0)) + βₘ(v) = 0.28*(v+27.0)/(exp((v+27.0)/5.0) - 1.0) + αₕ(v) = 0.128*exp((v+50.0)/18.0) + βₕ(v) = 4.0/(1.0 + exp(-(v+27.0)/5.0)) + αₙ(v) = 0.032*(v+52.0)/(1.0 - exp(-(v+52.0)/5.0)) + βₙ(v) = 0.5*exp(-(v+57.0)/40.0) + + m∞(v) = αₘ(v)/(αₘ(v) + βₘ(v)) + h∞(v) = αₕ(v)/(αₕ(v) + βₕ(v)) + n∞(v) = αₙ(v)/(αₙ(v) + βₙ(v)) + + τₘ(v) = 1.0/(αₘ(v) + βₘ(v)) + τₕ(v) = 1.0/(αₕ(v) + βₕ(v)) + τₙ(v) = 1.0/(αₙ(v) + βₙ(v)) + + gₐₘₚₐ(v) = 5*(1+tanh(v/4)) + + eqs = [D(V) ~ (Iₐₚₚ + Iₙₒᵢₛₑ - ḡₙₐ*m^3*h*(V - Eₙₐ) - ḡₖ*n^4*(V - Eₖ) - ḡₗ*(V - Eₗ) - jcn)/C, + D(m) ~ (m∞(V) - m)/τₘ(V), + D(h) ~ (h∞(V) - h)/τₕ(V), + D(n) ~ (n∞(V) - n)/τₙ(V), + D(sₐₘₚₐ) ~ gₐₘₚₐ(V)*(1-sₐₘₚₐ) - sₐₘₚₐ/τₑ + ] + + sys = ODESystem(eqs, t, sts, p; name=name) + + new(p, sys, namespace) + + end +end + +struct AdamINP <: AbstractAdamNeuron + params + system + namespace + + function AdamINP(;name, + namespace=nothing, + C=1.0, + Eₙₐ=50, + ḡₙₐ=100, + Eₖ=-100, + ḡₖ=80, + Eₗ=-67, + ḡₗ=0.05, + τᵢ=6, + Iₐₚₚ=0.1, + Iₙₒᵢₛₑ=0.0) + p = paramscoping(C=C, Eₙₐ=Eₙₐ, ḡₙₐ=ḡₙₐ, Eₖ=Eₖ, ḡₖ=ḡₖ, Eₗ=Eₗ, ḡₗ=ḡₗ, Iₐₚₚ=Iₐₚₚ, Iₙₒᵢₛₑ=Iₙₒᵢₛₑ) + C, Eₙₐ, ḡₙₐ, Eₖ, ḡₖ, Eₗ, ḡₗ, Iₐₚₚ, Iₙₒᵢₛₑ = p + sts = @variables V(t)=0.0 m(t)=0.0 h(t)=0.0 n(t)=0.0 sᵧ(t)=0.0 [output=true] jcn(t) [input=true] + + αₘ(v) = 0.32*(v+54.0)/(1.0 - exp(-(v+54.0)/4.0)) + βₘ(v) = 0.28*(v+27.0)/(exp((v+27.0)/5.0) - 1.0) + αₕ(v) = 0.128*exp((v+50.0)/18.0) + βₕ(v) = 4.0/(1.0 + exp(-(v+27.0)/5.0)) + αₙ(v) = 0.032*(v+52.0)/(1.0 - exp(-(v+52.0)/5.0)) + βₙ(v) = 0.5*exp(-(v+57.0)/40.0) + + m∞(v) = αₘ(v)/(αₘ(v) + βₘ(v)) + h∞(v) = αₕ(v)/(αₕ(v) + βₕ(v)) + n∞(v) = αₙ(v)/(αₙ(v) + βₙ(v)) + + τₘ(v) = 1.0/(αₘ(v) + βₘ(v)) + τₕ(v) = 1.0/(αₕ(v) + βₕ(v)) + τₙ(v) = 1.0/(αₙ(v) + βₙ(v)) + + gᵧ(v) = 2*(1+tanh(v/4)) + + eqs = [D(V) ~ (Iₐₚₚ + Iₙₒᵢₛₑ - ḡₙₐ*m^3*h*(V - Eₙₐ) - ḡₖ*n^4*(V - Eₖ) - ḡₗ*(V - Eₗ) - jcn)/C, + D(m) ~ (m∞(V) - m)/τₘ(V), + D(h) ~ (h∞(V) - h)/τₕ(V), + D(n) ~ (n∞(V) - n)/τₙ(V), + D(sᵧ) ~ gᵧ(V)*(1-sᵧ) - sᵧ/τᵢ + ] + + sys = ODESystem(eqs, t, sts, p; name=name) + + new(p, sys, namespace) + + end +end \ No newline at end of file diff --git a/src/blox/connections.jl b/src/blox/connections.jl index 5f097db2..88e22922 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -1222,3 +1222,38 @@ function Connector( return Connector(nameof(sys_src), nameof(sys_dest); equation=eq, weight=w) end + +function Connector( + blox_src::AdamPYR, + blox_dest::AbstractAdamNeuron; + kwargs... +) + + sys_pre = blox_src.system + sys_post = blox_dest.system + w = generate_weight_param(blox_src, blox_dest; kwargs...) + E_Exc = haskey(kwargs, :E_Exc) ? kwargs[:E_Exc] : 0.0 + s = only(outputs(blox_src; namespaced=true)) + + eq = sys_post.jcn ~ w*s*(sys_post.V - E_Exc) + + return Connector(nameof(sys_pre), nameof(sys_post); equation=eq, weight=w) + +end + +function Connector( + blox_src::AdamINP, + blox_dest::AbstractAdamNeuron; + kwargs... +) + + sys_pre = blox_src.system + sys_post = blox_dest.system + w = generate_weight_param(blox_src, blox_dest; kwargs...) + E_Inh = haskey(kwargs, :E_Inh) ? kwargs[:E_Inh] : -80.0 + s = only(outputs(blox_src; namespaced=true)) + + eq = sys_post.jcn ~ w*s*(sys_post.V - E_Inh) + + return Connector(nameof(sys_pre), nameof(sys_post); equation=eq, weight=w) +end \ No newline at end of file From 51f02febe9ee8cb87c73396a55014c9e924dd869 Mon Sep 17 00:00:00 2001 From: Alexandra Bardon <52323995+abardon@users.noreply.github.com> Date: Wed, 26 Feb 2025 16:45:08 +0000 Subject: [PATCH 02/18] test commit --- src/blox/adam_neurons.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/blox/adam_neurons.jl b/src/blox/adam_neurons.jl index 303c62a9..16bc85a4 100644 --- a/src/blox/adam_neurons.jl +++ b/src/blox/adam_neurons.jl @@ -102,4 +102,4 @@ struct AdamINP <: AbstractAdamNeuron new(p, sys, namespace) end -end \ No newline at end of file +end From 9fb631434c4e40aac2803184d0e39f34e5bf06bf Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 28 Feb 2025 20:47:19 -0500 Subject: [PATCH 03/18] Abstract types and custom heavisdie --- src/blox/adam_examples.jl | 2 ++ src/blox/adam_neurons.jl | 11 +++++++++++ 2 files changed, 13 insertions(+) diff --git a/src/blox/adam_examples.jl b/src/blox/adam_examples.jl index 13f15c5b..6a0f4903 100644 --- a/src/blox/adam_examples.jl +++ b/src/blox/adam_examples.jl @@ -59,6 +59,8 @@ struct AdamPYR <: AbstractAdamNeuron end end +struct AdamGLU <: + struct AdamINP <: AbstractAdamNeuron params system diff --git a/src/blox/adam_neurons.jl b/src/blox/adam_neurons.jl index 16bc85a4..064ed667 100644 --- a/src/blox/adam_neurons.jl +++ b/src/blox/adam_neurons.jl @@ -1,5 +1,15 @@ +using IfElse + abstract type AbstractAdamNeuron <: AbstractNeuronBlox end +abstract type AbstractReceptor <: AbstractBlox end +abstract type AbstractNeurotransmitter <: AbstractBlox end + +# Custom IfElse function to ensure differentiability so the solvers don't complain +function heaviside(x) + IfElse(x > 0, 1.0, 0.0) +end + struct AdamPYR <: AbstractAdamNeuron params system @@ -103,3 +113,4 @@ struct AdamINP <: AbstractAdamNeuron end end + From c0b710d9c55dcaa83ebf34981a697604314b49f4 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 28 Feb 2025 20:49:00 -0500 Subject: [PATCH 04/18] Add Glu and fix Systems --- src/blox/adam_neurons.jl | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/src/blox/adam_neurons.jl b/src/blox/adam_neurons.jl index 064ed667..95ee73ed 100644 --- a/src/blox/adam_neurons.jl +++ b/src/blox/adam_neurons.jl @@ -55,7 +55,7 @@ struct AdamPYR <: AbstractAdamNeuron D(sₐₘₚₐ) ~ gₐₘₚₐ(V)*(1-sₐₘₚₐ) - sₐₘₚₐ/τₑ ] - sys = ODESystem(eqs, t, sts, p; name=name) + sys = System(eqs, t, sts, p; name=name) new(p, sys, namespace) @@ -107,10 +107,32 @@ struct AdamINP <: AbstractAdamNeuron D(sᵧ) ~ gᵧ(V)*(1-sᵧ) - sᵧ/τᵢ ] - sys = ODESystem(eqs, t, sts, p; name=name) + sys = System(eqs, t, sts, p; name=name) new(p, sys, namespace) end end +struct AdamGlu <: AbstractNeurotransmitter + params + system + namespace + + function AdamGlu(;name, + namespace=nothing, + Glu_max = 1.0, + τ_Glu=1.2, + θ=10.0) + + p = paramscoping(Glu_max=Glu_max, τ_Glu=τ_Glu, θ=θ) + Glu_max, τ_Glu, θ = p + sts = @variables Glu(t)=0.0 [output=true] jcn(t) [input=true] + + eqs = [D(Glu) ~ Glu_max*heaviside(jcn - θ) - Glu/τ_Glu] + + sys = System(eqs, t, sts, p; name=name) + + new(p, sys, namespace) + end +end \ No newline at end of file From f2c5347ee42a8f88ed2eea6fca4aeeecd544c192 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 28 Feb 2025 21:02:50 -0500 Subject: [PATCH 05/18] Cleanup because these are moved to other files Also adding IfElse to Project.toml --- Project.toml | 1 + src/blox/adam_examples.jl | 145 +------------------------------------- 2 files changed, 2 insertions(+), 144 deletions(-) diff --git a/Project.toml b/Project.toml index 716c21af..4ea5411d 100644 --- a/Project.toml +++ b/Project.toml @@ -16,6 +16,7 @@ ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" GraphDynamics = "bcd5d0fe-e6b7-4ef1-9848-780c183c7f4c" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" +IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688" diff --git a/src/blox/adam_examples.jl b/src/blox/adam_examples.jl index 6a0f4903..183bfaa2 100644 --- a/src/blox/adam_examples.jl +++ b/src/blox/adam_examples.jl @@ -5,149 +5,6 @@ using Plots import Neuroblox: AbstractNeuronBlox, paramscoping using BenchmarkTools -abstract type AbstractAdamNeuron <: AbstractNeuronBlox end - -struct AdamPYR <: AbstractAdamNeuron - params - system - namespace - - function AdamPYR(;name, - namespace=nothing, - C=1.0, - Eₙₐ=50, - ḡₙₐ=100, - Eₖ=-100, - ḡₖ=80, - Eₗ=-67, - ḡₗ=0.05, - τₑ=1.5, - Iₐₚₚ=-0.25, - Iₙₒᵢₛₑ=0.0) - p = paramscoping(C=C, Eₙₐ=Eₙₐ, ḡₙₐ=ḡₙₐ, Eₖ=Eₖ, ḡₖ=ḡₖ, Eₗ=Eₗ, ḡₗ=ḡₗ, Iₐₚₚ=Iₐₚₚ, Iₙₒᵢₛₑ=Iₙₒᵢₛₑ) - C, Eₙₐ, ḡₙₐ, Eₖ, ḡₖ, Eₗ, ḡₗ, Iₐₚₚ, Iₙₒᵢₛₑ = p - sts = @variables V(t)=0.0 m(t)=0.0 h(t)=0.0 n(t)=0.0 sₐₘₚₐ(t)=0.0 [output=true] jcn(t) [input=true] - - αₘ(v) = 0.32*(v+54.0)/(1.0 - exp(-(v+54.0)/4.0)) - βₘ(v) = 0.28*(v+27.0)/(exp((v+27.0)/5.0) - 1.0) - αₕ(v) = 0.128*exp((v+50.0)/18.0) - βₕ(v) = 4.0/(1.0 + exp(-(v+27.0)/5.0)) - αₙ(v) = 0.032*(v+52.0)/(1.0 - exp(-(v+52.0)/5.0)) - βₙ(v) = 0.5*exp(-(v+57.0)/40.0) - - m∞(v) = αₘ(v)/(αₘ(v) + βₘ(v)) - h∞(v) = αₕ(v)/(αₕ(v) + βₕ(v)) - n∞(v) = αₙ(v)/(αₙ(v) + βₙ(v)) - - τₘ(v) = 1.0/(αₘ(v) + βₘ(v)) - τₕ(v) = 1.0/(αₕ(v) + βₕ(v)) - τₙ(v) = 1.0/(αₙ(v) + βₙ(v)) - - gₐₘₚₐ(v) = 5*(1+tanh(v/4)) - - eqs = [D(V) ~ (Iₐₚₚ + Iₙₒᵢₛₑ - ḡₙₐ*m^3*h*(V - Eₙₐ) - ḡₖ*n^4*(V - Eₖ) - ḡₗ*(V - Eₗ) - jcn)/C, - D(m) ~ (m∞(V) - m)/τₘ(V), - D(h) ~ (h∞(V) - h)/τₕ(V), - D(n) ~ (n∞(V) - n)/τₙ(V), - D(sₐₘₚₐ) ~ gₐₘₚₐ(V)*(1-sₐₘₚₐ) - sₐₘₚₐ/τₑ - ] - - sys = ODESystem(eqs, t, sts, p; name=name) - - new(p, sys, namespace) - - end -end - -struct AdamGLU <: - -struct AdamINP <: AbstractAdamNeuron - params - system - namespace - - function AdamINP(;name, - namespace=nothing, - C=1.0, - Eₙₐ=50, - ḡₙₐ=100, - Eₖ=-100, - ḡₖ=80, - Eₗ=-67, - ḡₗ=0.05, - τᵢ=6, - Iₐₚₚ=0.1, - Iₙₒᵢₛₑ=0.0) - p = paramscoping(C=C, Eₙₐ=Eₙₐ, ḡₙₐ=ḡₙₐ, Eₖ=Eₖ, ḡₖ=ḡₖ, Eₗ=Eₗ, ḡₗ=ḡₗ, Iₐₚₚ=Iₐₚₚ, Iₙₒᵢₛₑ=Iₙₒᵢₛₑ) - C, Eₙₐ, ḡₙₐ, Eₖ, ḡₖ, Eₗ, ḡₗ, Iₐₚₚ, Iₙₒᵢₛₑ = p - sts = @variables V(t)=0.0 m(t)=0.0 h(t)=0.0 n(t)=0.0 sᵧ(t)=0.0 [output=true] jcn(t) [input=true] - - αₘ(v) = 0.32*(v+54.0)/(1.0 - exp(-(v+54.0)/4.0)) - βₘ(v) = 0.28*(v+27.0)/(exp((v+27.0)/5.0) - 1.0) - αₕ(v) = 0.128*exp((v+50.0)/18.0) - βₕ(v) = 4.0/(1.0 + exp(-(v+27.0)/5.0)) - αₙ(v) = 0.032*(v+52.0)/(1.0 - exp(-(v+52.0)/5.0)) - βₙ(v) = 0.5*exp(-(v+57.0)/40.0) - - m∞(v) = αₘ(v)/(αₘ(v) + βₘ(v)) - h∞(v) = αₕ(v)/(αₕ(v) + βₕ(v)) - n∞(v) = αₙ(v)/(αₙ(v) + βₙ(v)) - - τₘ(v) = 1.0/(αₘ(v) + βₘ(v)) - τₕ(v) = 1.0/(αₕ(v) + βₕ(v)) - τₙ(v) = 1.0/(αₙ(v) + βₙ(v)) - - gᵧ(v) = 2*(1+tanh(v/4)) - - eqs = [D(V) ~ (Iₐₚₚ + Iₙₒᵢₛₑ - ḡₙₐ*m^3*h*(V - Eₙₐ) - ḡₖ*n^4*(V - Eₖ) - ḡₗ*(V - Eₗ) - jcn)/C, - D(m) ~ (m∞(V) - m)/τₘ(V), - D(h) ~ (h∞(V) - h)/τₕ(V), - D(n) ~ (n∞(V) - n)/τₙ(V), - D(sᵧ) ~ gᵧ(V)*(1-sᵧ) - sᵧ/τᵢ - ] - - sys = ODESystem(eqs, t, sts, p; name=name) - - new(p, sys, namespace) - - end -end - -function Connector( - blox_src::AdamPYR, - blox_dest::AbstractAdamNeuron; - kwargs... -) - - sys_pre = blox_src.system - sys_post = blox_dest.system - w = generate_weight_param(blox_src, blox_dest; kwargs...) - E_Exc = haskey(kwargs, :E_Exc) ? kwargs[:E_Exc] : 0.0 - s = only(outputs(sys_pre; namespaced=true)) - - eq = sys_post.jcn ~ w*s*(sys_post.V - E_Exc) - - return Connector(nameof(sys_pre), nameof(sys_post); equation=eq, weight=w) - -end - -function Connector( - blox_src::AdamINP, - blox_dest::AbstractAdamNeuron; - kwargs... -) - - sys_pre = blox_src.system - sys_post = blox_dest.system - w = generate_weight_param(blox_src, blox_dest; kwargs...) - E_Inh = haskey(kwargs, :E_Inh) ? kwargs[:E_Inh] : -80.0 - s = only(outputs(sys_pre; namespaced=true)) - - eq = sys_post.jcn ~ w*s*(sys_post.V - E_Inh) - - return Connector(nameof(sys_pre), nameof(sys_post); equation=eq, weight=w) -end - ḡᵢ = 0.5 ḡₑ = 0.2 @@ -166,7 +23,7 @@ for ne ∈ exci end end -tspan = (0.0, 2000.0) +tspan = (0.0, 500.0) begin @btime @named sys = system_from_graph(g, graphdynamics=true) @btime prob = ODEProblem(sys, [], tspan) From 6852c2a0e58027c80f4f72590354481f4cf9c340 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 28 Feb 2025 21:08:18 -0500 Subject: [PATCH 06/18] Oops fix merge from rebase --- src/GraphDynamicsInterop/neuron_interop.jl | 35 +++------------------- 1 file changed, 4 insertions(+), 31 deletions(-) diff --git a/src/GraphDynamicsInterop/neuron_interop.jl b/src/GraphDynamicsInterop/neuron_interop.jl index 55d5ec98..6dddda81 100644 --- a/src/GraphDynamicsInterop/neuron_interop.jl +++ b/src/GraphDynamicsInterop/neuron_interop.jl @@ -24,35 +24,6 @@ function define_neuron(sys; mod=@__MODULE__()) system = structural_simplify(sys.system; fully_determined=false) params = get_ps(system) t = Symbol(get_iv(system)) -function define_neurons() - for (name, T) ∈ [(:hhne, HHNeuronExciBlox) - (:hhni, HHNeuronInhibBlox) - (:hhni_msn_adam, HHNeuronInhib_MSN_Adam_Blox) - (:hhni_fsi_adam, HHNeuronInhib_FSI_Adam_Blox) - (:hhne_stn_adam, HHNeuronExci_STN_Adam_Blox) - (:hhni_GPe_adam, HHNeuronInhib_GPe_Adam_Blox) - (:ngei, NextGenerationEIBlox) - (:wc, WilsonCowan) - (:ho, HarmonicOscillator) - (:jr, JansenRit) # Note! Regular JansenRit can support delays, and I have not yet implemented this! - (:if, IFNeuron) - (:lif, LIFNeuron) - (:qif, QIFNeuron) - (:izh, IzhikevichNeuron) - (:lif_exci, LIFExciNeuron) - (:lif_inh, LIFInhNeuron) - (:pexci, PINGNeuronExci) - (:pinhib, PINGNeuronInhib) - (:VdP, VanDerPol{NonNoisy}) - (:VdPN, VanDerPol{Noisy}) - (:ko, KuramotoOscillator{NonNoisy}) - (:kon, KuramotoOscillator{Noisy}) - (:adpyr, AdamPYR) - (:adinp, AdamINP)] - sys = T(;name) - system = structural_simplify(sys.system; fully_determined=false) - params = get_ps(system) - t = Symbol(get_iv(system)) states = [s for s ∈ unknowns(system) if !MTK.isinput(s)] inputs = [s for s ∈ unknowns(system) if MTK.isinput(s)] @@ -198,7 +169,9 @@ for sys ∈ [HHNeuronExciBlox(name=:hhne) VanDerPol{NonNoisy}(name=:VdP) VanDerPol{Noisy}(name=:VdPN) KuramotoOscillator{NonNoisy}(name=:ko) - KuramotoOscillator{Noisy}(name=:kon)] + KuramotoOscillator{Noisy}(name=:kon) + AdamPYR(name=:adam_pyr) + AdamINP(name=:adam_inp)] define_neuron(sys) end @@ -309,4 +282,4 @@ function GraphDynamics.apply_discrete_event!(integrator, _, vparams, s::Subsyste params = get_params(s) vparams[] = @set params.jcn_ = jcn nothing -end +end \ No newline at end of file From da7bcc5f864a18ee18c14de49e4e7b77116e6342 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 28 Feb 2025 21:18:33 -0500 Subject: [PATCH 07/18] Turn off benchmarking for now --- src/blox/adam_examples.jl | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/blox/adam_examples.jl b/src/blox/adam_examples.jl index 183bfaa2..aaf0b8ab 100644 --- a/src/blox/adam_examples.jl +++ b/src/blox/adam_examples.jl @@ -3,7 +3,7 @@ using OrdinaryDiffEq using Random, Distributions using Plots import Neuroblox: AbstractNeuronBlox, paramscoping -using BenchmarkTools +#using BenchmarkTools ḡᵢ = 0.5 ḡₑ = 0.2 @@ -24,10 +24,14 @@ for ne ∈ exci end tspan = (0.0, 500.0) -begin - @btime @named sys = system_from_graph(g, graphdynamics=true) - @btime prob = ODEProblem(sys, [], tspan) - @btime sol = solve(prob, Tsit5(), saveat=0.5) -end +# begin +# @btime @named sys = system_from_graph(g, graphdynamics=true) +# @btime prob = ODEProblem(sys, [], tspan) +# @btime sol = solve(prob, Tsit5(), saveat=0.5) +# end + +@named sys = system_from_graph(g, graphdynamics=true) +prob = ODEProblem(sys, [], tspan) +sol = solve(prob, Tsit5(), saveat=0.5) plot(sol, idxs=1:5:(NE+NI)*5) \ No newline at end of file From 35885db59f1515ed67312ee9041bda48a6aaa60a Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 28 Feb 2025 21:35:12 -0500 Subject: [PATCH 08/18] Oops typo --- src/blox/adam_neurons.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/blox/adam_neurons.jl b/src/blox/adam_neurons.jl index 95ee73ed..ce3abf22 100644 --- a/src/blox/adam_neurons.jl +++ b/src/blox/adam_neurons.jl @@ -7,7 +7,7 @@ abstract type AbstractNeurotransmitter <: AbstractBlox end # Custom IfElse function to ensure differentiability so the solvers don't complain function heaviside(x) - IfElse(x > 0, 1.0, 0.0) + IfElse.ifelse(x > 0, 1.0, 0.0) end struct AdamPYR <: AbstractAdamNeuron From 27fc2021a4e70dc234a9346d106f2c1f2064f7d1 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 28 Feb 2025 21:42:09 -0500 Subject: [PATCH 09/18] Add glutamate blox for [Glu] --- src/Neuroblox.jl | 1 + src/blox/adam_examples.jl | 18 +++++++++++++++++- src/blox/adam_neurons.jl | 3 ++- src/blox/connections.jl | 16 ++++++++++++++++ 4 files changed, 36 insertions(+), 2 deletions(-) diff --git a/src/Neuroblox.jl b/src/Neuroblox.jl index 598584a1..3b356cad 100644 --- a/src/Neuroblox.jl +++ b/src/Neuroblox.jl @@ -262,4 +262,5 @@ export AdjacencyMatrix, Connector, connection_rule, connection_equations, connec export inputs, outputs, equations, unknowns, parameters, discrete_events export MetabolicHHNeuron export AdamPYR, AdamINP +export AdamGlu end diff --git a/src/blox/adam_examples.jl b/src/blox/adam_examples.jl index aaf0b8ab..b19fe4a8 100644 --- a/src/blox/adam_examples.jl +++ b/src/blox/adam_examples.jl @@ -34,4 +34,20 @@ tspan = (0.0, 500.0) prob = ODEProblem(sys, [], tspan) sol = solve(prob, Tsit5(), saveat=0.5) -plot(sol, idxs=1:5:(NE+NI)*5) \ No newline at end of file +plot(sol, idxs=1:5:(NE+NI)*5) + + + +## Testing Glu for threshold setting +## Commented out for now but useful for tuning later so leaving in the file +# exci = AdamPYR(name=:PYR, Iₐₚₚ=0.25) +# glur = AdamGlu(name=:Glu, θ=-59.0) + +# g = MetaDiGraph() +# add_edge!(g, exci => glur; weight=1.0) + +# tspan = (0.0, 500.0) +# @named sys = system_from_graph(g, graphdynamics=false) +# prob = ODEProblem(sys, [], tspan) +# sol = solve(prob, Tsit5(), saveat=0.5) +# plot(sol, idxs=6) \ No newline at end of file diff --git a/src/blox/adam_neurons.jl b/src/blox/adam_neurons.jl index ce3abf22..538ae379 100644 --- a/src/blox/adam_neurons.jl +++ b/src/blox/adam_neurons.jl @@ -114,6 +114,7 @@ struct AdamINP <: AbstractAdamNeuron end end +# Threshold θ is set to -59 mV so that the total impulse of an average spike is about 1.0 struct AdamGlu <: AbstractNeurotransmitter params system @@ -123,7 +124,7 @@ struct AdamGlu <: AbstractNeurotransmitter namespace=nothing, Glu_max = 1.0, τ_Glu=1.2, - θ=10.0) + θ=-59.0) p = paramscoping(Glu_max=Glu_max, τ_Glu=τ_Glu, θ=θ) Glu_max, τ_Glu, θ = p diff --git a/src/blox/connections.jl b/src/blox/connections.jl index 88e22922..897a5e07 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -1256,4 +1256,20 @@ function Connector( eq = sys_post.jcn ~ w*s*(sys_post.V - E_Inh) return Connector(nameof(sys_pre), nameof(sys_post); equation=eq, weight=w) +end + +function Connector( + blox_src::AdamPYR, + blox_dest::AdamGlu; + kwargs... +) + + sys_pre = blox_src.system + sys_post = blox_dest.system + w = generate_weight_param(blox_src, blox_dest; kwargs...) + + eq = sys_post.jcn ~ sys_pre.V + + return Connector(nameof(sys_pre), nameof(sys_post); equation=eq, weight=w) + end \ No newline at end of file From 70d1641463add2447963cc0799c86d43fca33d4b Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 28 Feb 2025 22:32:24 -0500 Subject: [PATCH 10/18] Initial NMDAR Still missing post-synaptic connector --- src/Neuroblox.jl | 1 + src/blox/adam_examples.jl | 14 ++++++++++ src/blox/adam_neurons.jl | 56 +++++++++++++++++++++++++++++++++++++++ src/blox/connections.jl | 35 ++++++++++++++++++++++++ 4 files changed, 106 insertions(+) diff --git a/src/Neuroblox.jl b/src/Neuroblox.jl index 3b356cad..a85d2f17 100644 --- a/src/Neuroblox.jl +++ b/src/Neuroblox.jl @@ -263,4 +263,5 @@ export inputs, outputs, equations, unknowns, parameters, discrete_events export MetabolicHHNeuron export AdamPYR, AdamINP export AdamGlu +export AdamNMDAR end diff --git a/src/blox/adam_examples.jl b/src/blox/adam_examples.jl index b19fe4a8..0ebcb80f 100644 --- a/src/blox/adam_examples.jl +++ b/src/blox/adam_examples.jl @@ -36,7 +36,21 @@ sol = solve(prob, Tsit5(), saveat=0.5) plot(sol, idxs=1:5:(NE+NI)*5) +exci = AdamPYR(name=:PYR, Iₐₚₚ=0.25) +glur = AdamGlu(name=:Glu, θ=-59.0) +nmda = AdamNMDAR(name=:NMDA) +exci2 = AdamPYR(name=:PYR2, Iₐₚₚ=0.33) +g = MetaDiGraph() +add_edge!(g, exci => glur; weight=1.0) +add_edge!(g, glur => nmda; weight=1.0) +add_edge!(g, exci2 => nmda; weight=1.0) + +tspan = (0.0, 500.0) +@named sys = system_from_graph(g, graphdynamics=false) +prob = ODEProblem(sys, [], tspan) +sol = solve(prob, Tsit5(), saveat=0.5) +plot(sol, idxs=11) ## Testing Glu for threshold setting ## Commented out for now but useful for tuning later so leaving in the file diff --git a/src/blox/adam_neurons.jl b/src/blox/adam_neurons.jl index 538ae379..da2d7750 100644 --- a/src/blox/adam_neurons.jl +++ b/src/blox/adam_neurons.jl @@ -134,6 +134,62 @@ struct AdamGlu <: AbstractNeurotransmitter sys = System(eqs, t, sts, p; name=name) + new(p, sys, namespace) + end +end + +struct AdamNMDAR <: AbstractReceptor + params + system + namespace + + function AdamNMDAR(;name, + namespace=nothing, + k_on=5, + k_off=0.0055, + k_r=0.0018, + k_d=0.0084, + k_unblock=5.4, + k_block=0.61, + α=0.0916, + β=0.0465) + + p = paramscoping(k_on=k_on, k_off=k_off, k_r=k_r, k_d=k_d, k_unblock=k_unblock, k_block=k_block, α=α, β=β) + k_on, k_off, k_r, k_d, k_unblock, k_block, α, β = p + + sts = @variables begin + Glu(t) + [input=true] + V(t) + [input=true] + C(t)=0.5 + C_A(t)=0.0 + C_AA(t)=0.0 + D_AA(t)=0.0 + O_AA(t)=0.0 + [output = true] + O_AAB(t)=0.0 + C_AAB(t)=0.0 + D_AAB(t)=0.0 + C_AB(t)=0.0 + C_B(t)=0.5 + end + + eqs = [ + D(C) ~ k_off*C_A - 2*k_on*Glu*C, + D(C_A) ~ 2*k_off*C_AA + 2*k_on*Glu*C - (k_on*Glu + k_off)*C_A, + D(C_AA) ~ k_on*Glu*C_A + α*O_AA + k_r*D_AA - (2*k_off + β + k_d)*C_AA, + D(D_AA) ~ k_d*C_AA - k_r*D_AA, + D(O_AA) ~ β*C_AA + k_unblock*exp(V/47)*O_AAB - (α + k_block*exp(-V/17))*O_AA, + D(O_AAB) ~ k_block*exp(-V/17)*O_AA + β*C_AAB - (k_unblock*exp(V/47) + α)*O_AAB, + D(C_AAB) ~ α*O_AAB + k_on*Glu*C_AB + k_r*D_AAB - (β + 2*k_off + k_d)*C_AAB, + D(D_AAB) ~ k_d*C_AAB - k_r*D_AAB, + D(C_AB) ~ 2*k_off*C_AAB + 2*k_on*Glu*C_B - (k_on*Glu + k_off)*C_AB, + D(C_B) ~ k_off*C_AB - 2*k_on*Glu*C_B + ] + + sys = System( eqs, t, sts, p; name=name) + new(p, sys, namespace) end end \ No newline at end of file diff --git a/src/blox/connections.jl b/src/blox/connections.jl index 897a5e07..136b5037 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -1272,4 +1272,39 @@ function Connector( return Connector(nameof(sys_pre), nameof(sys_post); equation=eq, weight=w) +end + +# An alternative thought - could lower threshold to have larger impulses and set a max here +# to prevent [Glu] > [Glu]_max. Not doing that for now because don't want to mess with the +# potential softmax issue (probably negligible though). +function Connector( + blox_src::AdamGlu, + blox_dest::AdamNMDAR; + kwargs... +) + + sys_pre = blox_src.system + sys_post = blox_dest.system + w = generate_weight_param(blox_src, blox_dest; kwargs...) + + eq = sys_post.Glu ~ sys_pre.Glu + + return Connector(nameof(sys_pre), nameof(sys_post); equation=eq, weight=w) + +end + +function Connector( + blox_src::AbstractAdamNeuron, + blox_dest::AdamNMDAR; + kwargs... +) + + sys_pre = blox_src.system + sys_post = blox_dest.system + w = generate_weight_param(blox_src, blox_dest; kwargs...) + + eq = sys_post.V ~ sys_pre.V + + return Connector(nameof(sys_pre), nameof(sys_post); equation=eq, weight=w) + end \ No newline at end of file From a1e326ac0f9830fb705d421e4743dbc51db7c61e Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 1 Mar 2025 15:47:50 -0500 Subject: [PATCH 11/18] NDMA receptor connection updates --- src/blox/adam_neurons.jl | 6 ++++-- src/blox/connections.jl | 16 ++++++++++++++++ 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/src/blox/adam_neurons.jl b/src/blox/adam_neurons.jl index da2d7750..1a40e542 100644 --- a/src/blox/adam_neurons.jl +++ b/src/blox/adam_neurons.jl @@ -152,7 +152,9 @@ struct AdamNMDAR <: AbstractReceptor k_unblock=5.4, k_block=0.61, α=0.0916, - β=0.0465) + β=0.0465, + g_nmda=8.5, + E_nmda=0.0) p = paramscoping(k_on=k_on, k_off=k_off, k_r=k_r, k_d=k_d, k_unblock=k_unblock, k_block=k_block, α=α, β=β) k_on, k_off, k_r, k_d, k_unblock, k_block, α, β = p @@ -188,7 +190,7 @@ struct AdamNMDAR <: AbstractReceptor D(C_B) ~ k_off*C_AB - 2*k_on*Glu*C_B ] - sys = System( eqs, t, sts, p; name=name) + sys = System(eqs, t, sts, p; name=name) new(p, sys, namespace) end diff --git a/src/blox/connections.jl b/src/blox/connections.jl index 136b5037..47033520 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -1307,4 +1307,20 @@ function Connector( return Connector(nameof(sys_pre), nameof(sys_post); equation=eq, weight=w) +end + +function Connector( + blox_src::AdamNMDAR, + blox_dest::AbstractAdamNeuron; + kwargs... +) + + sys_pre = blox_src.system + sys_post = blox_dest.system + w = generate_weight_param(blox_src, blox_dest; kwargs...) + + eq = sys_post.jcn ~ w*sys_pre.g_NMDA*sys_pre.OAA*(sys_post.V - sys_pre.E_NMDA) + + return Connector(nameof(sys_pre), nameof(sys_post); equation=eq, weight=w) + end \ No newline at end of file From b6b9217e2364ff3ee059e041b8736ef87f3b1a92 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 1 Mar 2025 15:52:33 -0500 Subject: [PATCH 12/18] Oops need to add to parameter list --- src/blox/adam_examples.jl | 1 + src/blox/adam_neurons.jl | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/blox/adam_examples.jl b/src/blox/adam_examples.jl index 0ebcb80f..7d97bc1a 100644 --- a/src/blox/adam_examples.jl +++ b/src/blox/adam_examples.jl @@ -45,6 +45,7 @@ g = MetaDiGraph() add_edge!(g, exci => glur; weight=1.0) add_edge!(g, glur => nmda; weight=1.0) add_edge!(g, exci2 => nmda; weight=1.0) +add_edge!(g, nmda => exci2; weight=1.0) tspan = (0.0, 500.0) @named sys = system_from_graph(g, graphdynamics=false) diff --git a/src/blox/adam_neurons.jl b/src/blox/adam_neurons.jl index 1a40e542..4f745b6f 100644 --- a/src/blox/adam_neurons.jl +++ b/src/blox/adam_neurons.jl @@ -156,8 +156,8 @@ struct AdamNMDAR <: AbstractReceptor g_nmda=8.5, E_nmda=0.0) - p = paramscoping(k_on=k_on, k_off=k_off, k_r=k_r, k_d=k_d, k_unblock=k_unblock, k_block=k_block, α=α, β=β) - k_on, k_off, k_r, k_d, k_unblock, k_block, α, β = p + p = paramscoping(k_on=k_on, k_off=k_off, k_r=k_r, k_d=k_d, k_unblock=k_unblock, k_block=k_block, α=α, β=β, g_nmda=g_nmda, E_nmda=E_nmda) + k_on, k_off, k_r, k_d, k_unblock, k_block, α, β, g_nmda, E_nmda = p sts = @variables begin Glu(t) From dfd76f5ed0ff05ceb7a2742433e839aa16851227 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 1 Mar 2025 15:55:45 -0500 Subject: [PATCH 13/18] Annnd fix param calls --- src/blox/connections.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/blox/connections.jl b/src/blox/connections.jl index 47033520..985fed5c 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -1319,7 +1319,7 @@ function Connector( sys_post = blox_dest.system w = generate_weight_param(blox_src, blox_dest; kwargs...) - eq = sys_post.jcn ~ w*sys_pre.g_NMDA*sys_pre.OAA*(sys_post.V - sys_pre.E_NMDA) + eq = sys_post.jcn ~ w*blox_src.params.g_NMDA*sys_pre.OAA*(sys_post.V - blox_src.params.E_NMDA) return Connector(nameof(sys_pre), nameof(sys_post); equation=eq, weight=w) From 30142293f9d6c3385600bb8a4dfe3f0676599d50 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 1 Mar 2025 16:09:06 -0500 Subject: [PATCH 14/18] Ok last param fix --- src/blox/adam_examples.jl | 2 +- src/blox/adam_neurons.jl | 9 ++++----- src/blox/connections.jl | 5 +++-- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/blox/adam_examples.jl b/src/blox/adam_examples.jl index 7d97bc1a..0404454b 100644 --- a/src/blox/adam_examples.jl +++ b/src/blox/adam_examples.jl @@ -51,7 +51,7 @@ tspan = (0.0, 500.0) @named sys = system_from_graph(g, graphdynamics=false) prob = ODEProblem(sys, [], tspan) sol = solve(prob, Tsit5(), saveat=0.5) -plot(sol, idxs=11) +plot(sol) ## Testing Glu for threshold setting ## Commented out for now but useful for tuning later so leaving in the file diff --git a/src/blox/adam_neurons.jl b/src/blox/adam_neurons.jl index 4f745b6f..194d0273 100644 --- a/src/blox/adam_neurons.jl +++ b/src/blox/adam_neurons.jl @@ -152,12 +152,11 @@ struct AdamNMDAR <: AbstractReceptor k_unblock=5.4, k_block=0.61, α=0.0916, - β=0.0465, - g_nmda=8.5, - E_nmda=0.0) + β=0.0465 + ) - p = paramscoping(k_on=k_on, k_off=k_off, k_r=k_r, k_d=k_d, k_unblock=k_unblock, k_block=k_block, α=α, β=β, g_nmda=g_nmda, E_nmda=E_nmda) - k_on, k_off, k_r, k_d, k_unblock, k_block, α, β, g_nmda, E_nmda = p + p = paramscoping(k_on=k_on, k_off=k_off, k_r=k_r, k_d=k_d, k_unblock=k_unblock, k_block=k_block, α=α, β=β) + k_on, k_off, k_r, k_d, k_unblock, k_block, α, β = p sts = @variables begin Glu(t) diff --git a/src/blox/connections.jl b/src/blox/connections.jl index 985fed5c..b8df5d92 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -1318,8 +1318,9 @@ function Connector( sys_pre = blox_src.system sys_post = blox_dest.system w = generate_weight_param(blox_src, blox_dest; kwargs...) - - eq = sys_post.jcn ~ w*blox_src.params.g_NMDA*sys_pre.OAA*(sys_post.V - blox_src.params.E_NMDA) + g_NMDA = haskey(kwargs, :g_NMDA) ? kwargs[:g_NMDA] : 8.5 + E_NMDA = haskey(kwargs, :E_NMDA) ? kwargs[:E_NMDA] : 0.0 + eq = sys_post.jcn ~ w*g_NMDA*sys_pre.O_AA*(sys_post.V - E_NMDA) return Connector(nameof(sys_pre), nameof(sys_post); equation=eq, weight=w) From 0628806debb3131f53472aaac026ae6a481ffe77 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 1 Mar 2025 22:58:21 -0500 Subject: [PATCH 15/18] Very ugly but very performant NMDAR with GraphDynamics --- .../GraphDynamicsInterop.jl | 6 +- .../connection_interop.jl | 28 ++++- src/GraphDynamicsInterop/neuron_interop.jl | 30 ++++- src/blox/adam_examples.jl | 103 +++++++++++++++++- 4 files changed, 157 insertions(+), 10 deletions(-) diff --git a/src/GraphDynamicsInterop/GraphDynamicsInterop.jl b/src/GraphDynamicsInterop/GraphDynamicsInterop.jl index e590bd96..205a3dfd 100644 --- a/src/GraphDynamicsInterop/GraphDynamicsInterop.jl +++ b/src/GraphDynamicsInterop/GraphDynamicsInterop.jl @@ -54,7 +54,11 @@ using ..Neuroblox: VanDerPol, AdamPYR, AdamINP, - AbstractAdamNeuron + AbstractAdamNeuron, + AdamNMDAR, + AbstractReceptor, + AdamGlu, + AbstractNeurotransmitter using GraphDynamics: GraphDynamics, diff --git a/src/GraphDynamicsInterop/connection_interop.jl b/src/GraphDynamicsInterop/connection_interop.jl index 508955b5..a144c15b 100644 --- a/src/GraphDynamicsInterop/connection_interop.jl +++ b/src/GraphDynamicsInterop/connection_interop.jl @@ -1050,7 +1050,6 @@ end # #------------------------- # Adam Network # #------------------------- -# PING Network struct AdamConnection <: ConnectionRule w::Float64 V_E::Float64 @@ -1071,14 +1070,33 @@ end function (c::AdamConnection)(blox_src::Subsystem{AdamPYR}, blox_dst::Subsystem{<:AbstractAdamNeuron}) (; w, V_E) = c - (;sₐₘₚₐ) = blox_src - (;V) = blox_dst + (; sₐₘₚₐ) = blox_src + (; V) = blox_dst (; jcn = w * sₐₘₚₐ * (V - V_E)) end function (c::AdamConnection)(blox_src::Subsystem{AdamINP}, blox_dst::Subsystem{<:AbstractAdamNeuron}) (; w, V_I) = c - (;sᵧ) = blox_src - (;V) = blox_dst + (; sᵧ) = blox_src + (; V) = blox_dst (; jcn = w * sᵧ * (V - V_I)) end + +function (c::BasicConnection)(blox_src::Subsystem{AdamPYR}, blox_dst::Subsystem{AdamGlu}) + w = c.weight + (; V) = blox_src + (; jcn = V) +end + +# using Accessors +# function (c::AdamConnection)(sys_src::Subsystem{AdamGlu}, sys_dst::Subsystem{AdamNMDAR}) +# w = c.weight +# input = initialize_input(sys_dest) +# @set input.Glu = sys_src.Glu +# end + +# function (c::AdamConnection)(sys_src::Subsystem{<:AbstractAdamNeuron}, sys_dst::Subsystem{AdamNMDAR}) +# w = c.weight +# input = initialize_input(sys_dest) +# @set input.V = sys_src.V +# end diff --git a/src/GraphDynamicsInterop/neuron_interop.jl b/src/GraphDynamicsInterop/neuron_interop.jl index 6dddda81..cb202a05 100644 --- a/src/GraphDynamicsInterop/neuron_interop.jl +++ b/src/GraphDynamicsInterop/neuron_interop.jl @@ -171,7 +171,9 @@ for sys ∈ [HHNeuronExciBlox(name=:hhne) KuramotoOscillator{NonNoisy}(name=:ko) KuramotoOscillator{Noisy}(name=:kon) AdamPYR(name=:adam_pyr) - AdamINP(name=:adam_inp)] + AdamINP(name=:adam_inp) + AdamGlu(name=:adam_glu) + AdamNMDAR(name=:adam_nmdar)] define_neuron(sys) end @@ -282,4 +284,28 @@ function GraphDynamics.apply_discrete_event!(integrator, _, vparams, s::Subsyste params = get_params(s) vparams[] = @set params.jcn_ = jcn nothing -end \ No newline at end of file +end + +# GraphDynamics.initialize_input(s::Subsystem{AdamNMDAR}) = (; Glu = 0.0, V = 0.0) + +# function GraphDynamics.subsystem_differential(sys::Subsystem{AdamNMDAR}, inputs, t) +# # Unpack the system +# (; C, C_A, C_AA, D_AA, O_AA, O_AAB, C_AAB, D_AAB, C_AB, C_B) = sys +# (; k_on, k_off, k_r, k_d, k_unblock, k_block, α, β) = sys +# (; Glu, V) = inputs + +# # Differential equations +# SubsystemStates{AdamNMDAR}( +# #=d/dt=# C = k_off*C_A - 2*k_on*Glu*C, +# #=d/dt=# C_A = 2*k_off*C_AA + 2*k_on*Glu*C - (k_on*Glu + k_off)*C_A, +# #=d/dt=# C_AA = k_on*Glu*C_A + α*O_AA + k_r*D_AA - (2*k_off + β + k_d)*C_AA, +# #=d/dt=# D_AA = k_d*C_AA - k_r*D_AA, +# #=d/dt=# O_AA = k_r*D_AA - α*O_AA, +# #=d/dt=# O_AAB = k_unblock*C_AAB - k_block*O_AAB, +# #=d/dt=# C_AAB = k_block*O_AAB - k_unblock*C_AAB, +# #=d/dt=# D_AAB = k_d*C_AAB - k_r*D_AAB, +# #=d/dt=# C_AB = k_off*C_AAB - 2*k_on*Glu*C_AB, +# #=d/dt=# C_B = k_off*C_AB - 2*k_on*Glu*C_B +# ) + +# end \ No newline at end of file diff --git a/src/blox/adam_examples.jl b/src/blox/adam_examples.jl index 0404454b..9fe5ce2a 100644 --- a/src/blox/adam_examples.jl +++ b/src/blox/adam_examples.jl @@ -2,7 +2,8 @@ using Neuroblox using OrdinaryDiffEq using Random, Distributions using Plots -import Neuroblox: AbstractNeuronBlox, paramscoping +using GraphDynamics +#import Neuroblox: AbstractNeuronBlox, paramscoping #using BenchmarkTools ḡᵢ = 0.5 @@ -36,6 +37,81 @@ sol = solve(prob, Tsit5(), saveat=0.5) plot(sol, idxs=1:5:(NE+NI)*5) +using Neuroblox.GraphDynamicsInterop +GraphDynamicsInterop.issupported(::AdamNMDAR) = true +GraphDynamicsInterop.components(v::AdamNMDAR) = (v,) + +function GraphDynamicsInterop.to_subsystem(v::AdamNMDAR) + # Extract default parameter values + k_on = GraphDynamicsInterop.recursive_getdefault(v.system.k_on) + k_off = GraphDynamicsInterop.recursive_getdefault(v.system.k_off) + k_r = GraphDynamicsInterop.recursive_getdefault(v.system.k_r) + k_d = GraphDynamicsInterop.recursive_getdefault(v.system.k_d) + k_unblock = GraphDynamicsInterop.recursive_getdefault(v.system.k_unblock) + k_block = GraphDynamicsInterop.recursive_getdefault(v.system.k_block) + α = GraphDynamicsInterop.recursive_getdefault(v.system.α) + β = GraphDynamicsInterop.recursive_getdefault(v.system.β) + + params = SubsystemParams{AdamNMDAR}(; k_on, k_off, k_r, k_d, k_unblock, k_block, α, β) + + # Extract the default values of states + C = GraphDynamicsInterop.recursive_getdefault(v.system.C) + C_A = GraphDynamicsInterop.recursive_getdefault(v.system.C_A) + C_AA = GraphDynamicsInterop.recursive_getdefault(v.system.C_AA) + D_AA = GraphDynamicsInterop.recursive_getdefault(v.system.D_AA) + O_AA = GraphDynamicsInterop.recursive_getdefault(v.system.O_AA) + O_AAB = GraphDynamicsInterop.recursive_getdefault(v.system.O_AAB) + C_AAB = GraphDynamicsInterop.recursive_getdefault(v.system.C_AAB) + D_AAB = GraphDynamicsInterop.recursive_getdefault(v.system.D_AAB) + C_AB = GraphDynamicsInterop.recursive_getdefault(v.system.C_AB) + C_B = GraphDynamicsInterop.recursive_getdefault(v.system.C_B) + states = SubsystemStates{AdamNMDAR}(; C, C_A, C_AA, D_AA, O_AA, O_AAB, C_AAB, D_AAB, C_AB, C_B) + + Subsystem(states, params) +end + +GraphDynamics.initialize_input(s::Subsystem{AdamNMDAR}) = (; Glu = 0.0, V = 0.0) + +function GraphDynamics.subsystem_differential(s::Subsystem{AdamNMDAR}, inputs, t) + # Unpack + (; Glu, V) = inputs + (; C, C_A, C_AA, D_AA, O_AA, O_AAB, C_AAB, D_AAB, C_AB, C_B) = s + (; k_on, k_off, k_r, k_d, k_unblock, k_block, α, β) = s + return SubsystemStates{AdamNMDAR}( + #=d/dt=# C = k_off*C_A - 2*k_on*Glu*C, + #=d/dt=# C_A = 2*k_off*C_AA + 2*k_on*Glu*C - (k_on*Glu + k_off)*C_A, + #=d/dt=# C_AA = k_on*Glu*C_A + α*O_AA + k_r*D_AA - (2*k_off + β + k_d)*C_AA, + #=d/dt=# D_AA = k_d*C_AA - k_r*D_AA, + #=d/dt=# O_AA = k_r*D_AA - α*O_AA, + #=d/dt=# O_AAB = k_unblock*C_AAB - k_block*O_AAB, + #=d/dt=# C_AAB = k_block*O_AAB - k_unblock*C_AAB, + #=d/dt=# D_AAB = k_d*C_AAB - k_r*D_AAB, + #=d/dt=# C_AB = k_off*C_AAB - 2*k_on*Glu*C_AB, + #=d/dt=# C_B = k_off*C_AB - 2*k_on*Glu*C_B + ) +end + +function (c::GraphDynamicsInterop.BasicConnection)(sys_src::Subsystem{AdamGlu}, sys_dst::Subsystem{AdamNMDAR}, t) + w = c.weight + Glu = sys_src.Glu + V = 0.0 + (; Glu, V) +end + +function (c::GraphDynamicsInterop.BasicConnection)(sys_src::Subsystem{<:Neuroblox.AbstractAdamNeuron}, sys_dst::Subsystem{AdamNMDAR}, t) + w = c.weight + Glu = 0.0 + V = sys_src.V + (; Glu, V) +end + +function (c::GraphDynamicsInterop.BasicConnection)(sys_src::Subsystem{AdamNMDAR}, sys_dst::Subsystem{<:Neuroblox.AbstractAdamNeuron}, t) + w = c.weight + O_AA = sys_src.O_AA + jcn = w*O_AA*sys_dst.V + (; jcn) +end + exci = AdamPYR(name=:PYR, Iₐₚₚ=0.25) glur = AdamGlu(name=:Glu, θ=-59.0) nmda = AdamNMDAR(name=:NMDA) @@ -48,11 +124,34 @@ add_edge!(g, exci2 => nmda; weight=1.0) add_edge!(g, nmda => exci2; weight=1.0) tspan = (0.0, 500.0) -@named sys = system_from_graph(g, graphdynamics=false) +@named sys = system_from_graph(g, graphdynamics=true) prob = ODEProblem(sys, [], tspan) sol = solve(prob, Tsit5(), saveat=0.5) plot(sol) +NE = 80 +NI = 80 + +exci = [AdamPYR(name=Symbol("PYR$i"), Iₐₚₚ=rand(Normal(0.25, 0.05))) for i in 1:NE] +inhi = [AdamINP(name=Symbol("INP$i"), Iₐₚₚ=rand(Normal(0.3, 0.05))) for i in 1:NI] # bump up to 0.3 +nmdar = [AdamNMDAR(name=Symbol("NMDA$i")) for i in 1:NE] +glu = [AdamGlu(name=Symbol("Glu$i")) for i in 1:NI] + +g = MetaDiGraph() + +for i in axes(exci, 1) + add_edge!(g, exci[i] => glu[i]; weight=1.0) + add_edge!(g, glu[i] => nmdar[i]; weight=1.0) + add_edge!(g, inhi[i] => nmdar[i]; weight=1.0) + add_edge!(g, nmdar[i] => inhi[i]; weight=1.0) +end + +tspan = (0.0, 500.0) +@time @named sys = system_from_graph(g, graphdynamics=true) +@time prob = ODEProblem(sys, [], tspan) +@time sol = solve(prob, Tsit5(), saveat=0.5) + + ## Testing Glu for threshold setting ## Commented out for now but useful for tuning later so leaving in the file # exci = AdamPYR(name=:PYR, Iₐₚₚ=0.25) From 2632706813fae7c5ddf4c9441c18855a13966b04 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 1 Mar 2025 23:00:56 -0500 Subject: [PATCH 16/18] Update N to be more realistic --- src/blox/adam_examples.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/blox/adam_examples.jl b/src/blox/adam_examples.jl index 9fe5ce2a..99fee4a2 100644 --- a/src/blox/adam_examples.jl +++ b/src/blox/adam_examples.jl @@ -129,8 +129,8 @@ prob = ODEProblem(sys, [], tspan) sol = solve(prob, Tsit5(), saveat=0.5) plot(sol) -NE = 80 -NI = 80 +NE = 800 +NI = 800 exci = [AdamPYR(name=Symbol("PYR$i"), Iₐₚₚ=rand(Normal(0.25, 0.05))) for i in 1:NE] inhi = [AdamINP(name=Symbol("INP$i"), Iₐₚₚ=rand(Normal(0.3, 0.05))) for i in 1:NI] # bump up to 0.3 From fbf4cab6e77ede13d850d2d51dcabbe7277e1e72 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sun, 2 Mar 2025 15:24:13 -0500 Subject: [PATCH 17/18] Helper function and balancing E/I a bit --- src/blox/adam_examples.jl | 39 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/src/blox/adam_examples.jl b/src/blox/adam_examples.jl index 99fee4a2..0724adfd 100644 --- a/src/blox/adam_examples.jl +++ b/src/blox/adam_examples.jl @@ -151,6 +151,45 @@ tspan = (0.0, 500.0) @time prob = ODEProblem(sys, [], tspan) @time sol = solve(prob, Tsit5(), saveat=0.5) +function make_nmda_edge!(g, prenrn, postnrn) + glu = AdamGlu(name=Symbol("Glu$(prenrn.name)_$(postnrn.name)")) + nmda = AdamNMDAR(name=Symbol("NMDA$(prenrn.name)_$(postnrn.name)")) + add_edge!(g, prenrn => glu; weight=1.0) + add_edge!(g, glu => nmda; weight=1.0) + add_edge!(g, postnrn => nmda; weight=1.0) + add_edge!(g, nmda => postnrn; weight=1.0) +end + +NE = 80 +NI = 20 + +exci = [AdamPYR(name=Symbol("PYR$i"), Iₐₚₚ=rand(Normal(1.5, 0.05))) for i in 1:NE] +inhi = [AdamINP(name=Symbol("INP$i"), Iₐₚₚ=rand(Normal(0.1, 0.05))) for i in 1:NI] # bump up to 0.3 + +g = MetaDiGraph() + +for ne ∈ exci + for ni ∈ inhi + make_nmda_edge!(g, ne, ni) + end +end + +for ne ∈ exci + for ne ∈ exci + add_edge!(g, ne => ne; weight=1.0) + end +end + +for ni ∈ inhi + for ne ∈ exci[1:20] + add_edge!(g, ni => ne; weight=ḡᵢ/NI) + end +end + +tspan = (0.0, 500.0) +@time @named sys = system_from_graph(g, graphdynamics=true) +@time prob = ODEProblem(sys, [], tspan) +@time sol = solve(prob, Tsit5(), saveat=0.5) ## Testing Glu for threshold setting ## Commented out for now but useful for tuning later so leaving in the file From 22da6e338c3a39d6562ccbe146662fba01285e10 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Mon, 3 Mar 2025 12:40:19 -0500 Subject: [PATCH 18/18] Tidy up GraphDynamics functions --- src/blox/adam_examples.jl | 155 +++++++------------------- src/blox/adams_graphdynamics_draft.jl | 88 +++++++++++++++ 2 files changed, 129 insertions(+), 114 deletions(-) create mode 100644 src/blox/adams_graphdynamics_draft.jl diff --git a/src/blox/adam_examples.jl b/src/blox/adam_examples.jl index 0724adfd..e318bcc8 100644 --- a/src/blox/adam_examples.jl +++ b/src/blox/adam_examples.jl @@ -2,13 +2,51 @@ using Neuroblox using OrdinaryDiffEq using Random, Distributions using Plots -using GraphDynamics #import Neuroblox: AbstractNeuronBlox, paramscoping #using BenchmarkTools +include("adams_graphdynamics_draft.jl") + ḡᵢ = 0.5 ḡₑ = 0.2 +NE = 80 +NI = 20 + +exci = [AdamPYR(name=Symbol("PYR$i"), Iₐₚₚ=rand(Normal(1.5, 0.05))) for i in 1:NE] +inhi = [AdamINP(name=Symbol("INP$i"), Iₐₚₚ=rand(Normal(0.1, 0.05))) for i in 1:NI] # bump up to 0.3 + +g = MetaDiGraph() + +for ne ∈ exci + for ni ∈ inhi + make_nmda_edge!(g, ne, ni) + end +end + +for ne ∈ exci + for ne ∈ exci + add_edge!(g, ne => ne; weight=1.0) + end +end + +for ni ∈ inhi + for ne ∈ exci[1:20] + add_edge!(g, ni => ne; weight=ḡᵢ/NI) + end +end + +begin + tspan = (0.0, 1000.0) + @time sys = system_from_graph(g, graphdynamics=true) + @time prob = ODEProblem(sys, [], tspan) + @time sol = solve(prob, Tsit5(), saveat=0.5) +end + + +### Older tests + +## Test network without NMDAR connections NI = 20 NE = 80 @@ -37,80 +75,8 @@ sol = solve(prob, Tsit5(), saveat=0.5) plot(sol, idxs=1:5:(NE+NI)*5) -using Neuroblox.GraphDynamicsInterop -GraphDynamicsInterop.issupported(::AdamNMDAR) = true -GraphDynamicsInterop.components(v::AdamNMDAR) = (v,) - -function GraphDynamicsInterop.to_subsystem(v::AdamNMDAR) - # Extract default parameter values - k_on = GraphDynamicsInterop.recursive_getdefault(v.system.k_on) - k_off = GraphDynamicsInterop.recursive_getdefault(v.system.k_off) - k_r = GraphDynamicsInterop.recursive_getdefault(v.system.k_r) - k_d = GraphDynamicsInterop.recursive_getdefault(v.system.k_d) - k_unblock = GraphDynamicsInterop.recursive_getdefault(v.system.k_unblock) - k_block = GraphDynamicsInterop.recursive_getdefault(v.system.k_block) - α = GraphDynamicsInterop.recursive_getdefault(v.system.α) - β = GraphDynamicsInterop.recursive_getdefault(v.system.β) - - params = SubsystemParams{AdamNMDAR}(; k_on, k_off, k_r, k_d, k_unblock, k_block, α, β) - - # Extract the default values of states - C = GraphDynamicsInterop.recursive_getdefault(v.system.C) - C_A = GraphDynamicsInterop.recursive_getdefault(v.system.C_A) - C_AA = GraphDynamicsInterop.recursive_getdefault(v.system.C_AA) - D_AA = GraphDynamicsInterop.recursive_getdefault(v.system.D_AA) - O_AA = GraphDynamicsInterop.recursive_getdefault(v.system.O_AA) - O_AAB = GraphDynamicsInterop.recursive_getdefault(v.system.O_AAB) - C_AAB = GraphDynamicsInterop.recursive_getdefault(v.system.C_AAB) - D_AAB = GraphDynamicsInterop.recursive_getdefault(v.system.D_AAB) - C_AB = GraphDynamicsInterop.recursive_getdefault(v.system.C_AB) - C_B = GraphDynamicsInterop.recursive_getdefault(v.system.C_B) - states = SubsystemStates{AdamNMDAR}(; C, C_A, C_AA, D_AA, O_AA, O_AAB, C_AAB, D_AAB, C_AB, C_B) - - Subsystem(states, params) -end - -GraphDynamics.initialize_input(s::Subsystem{AdamNMDAR}) = (; Glu = 0.0, V = 0.0) - -function GraphDynamics.subsystem_differential(s::Subsystem{AdamNMDAR}, inputs, t) - # Unpack - (; Glu, V) = inputs - (; C, C_A, C_AA, D_AA, O_AA, O_AAB, C_AAB, D_AAB, C_AB, C_B) = s - (; k_on, k_off, k_r, k_d, k_unblock, k_block, α, β) = s - return SubsystemStates{AdamNMDAR}( - #=d/dt=# C = k_off*C_A - 2*k_on*Glu*C, - #=d/dt=# C_A = 2*k_off*C_AA + 2*k_on*Glu*C - (k_on*Glu + k_off)*C_A, - #=d/dt=# C_AA = k_on*Glu*C_A + α*O_AA + k_r*D_AA - (2*k_off + β + k_d)*C_AA, - #=d/dt=# D_AA = k_d*C_AA - k_r*D_AA, - #=d/dt=# O_AA = k_r*D_AA - α*O_AA, - #=d/dt=# O_AAB = k_unblock*C_AAB - k_block*O_AAB, - #=d/dt=# C_AAB = k_block*O_AAB - k_unblock*C_AAB, - #=d/dt=# D_AAB = k_d*C_AAB - k_r*D_AAB, - #=d/dt=# C_AB = k_off*C_AAB - 2*k_on*Glu*C_AB, - #=d/dt=# C_B = k_off*C_AB - 2*k_on*Glu*C_B - ) -end - -function (c::GraphDynamicsInterop.BasicConnection)(sys_src::Subsystem{AdamGlu}, sys_dst::Subsystem{AdamNMDAR}, t) - w = c.weight - Glu = sys_src.Glu - V = 0.0 - (; Glu, V) -end - -function (c::GraphDynamicsInterop.BasicConnection)(sys_src::Subsystem{<:Neuroblox.AbstractAdamNeuron}, sys_dst::Subsystem{AdamNMDAR}, t) - w = c.weight - Glu = 0.0 - V = sys_src.V - (; Glu, V) -end -function (c::GraphDynamicsInterop.BasicConnection)(sys_src::Subsystem{AdamNMDAR}, sys_dst::Subsystem{<:Neuroblox.AbstractAdamNeuron}, t) - w = c.weight - O_AA = sys_src.O_AA - jcn = w*O_AA*sys_dst.V - (; jcn) -end +### Testing single neuron connections exci = AdamPYR(name=:PYR, Iₐₚₚ=0.25) glur = AdamGlu(name=:Glu, θ=-59.0) @@ -129,6 +95,7 @@ prob = ODEProblem(sys, [], tspan) sol = solve(prob, Tsit5(), saveat=0.5) plot(sol) +### Testing multiple neuron connections NE = 800 NI = 800 @@ -151,46 +118,6 @@ tspan = (0.0, 500.0) @time prob = ODEProblem(sys, [], tspan) @time sol = solve(prob, Tsit5(), saveat=0.5) -function make_nmda_edge!(g, prenrn, postnrn) - glu = AdamGlu(name=Symbol("Glu$(prenrn.name)_$(postnrn.name)")) - nmda = AdamNMDAR(name=Symbol("NMDA$(prenrn.name)_$(postnrn.name)")) - add_edge!(g, prenrn => glu; weight=1.0) - add_edge!(g, glu => nmda; weight=1.0) - add_edge!(g, postnrn => nmda; weight=1.0) - add_edge!(g, nmda => postnrn; weight=1.0) -end - -NE = 80 -NI = 20 - -exci = [AdamPYR(name=Symbol("PYR$i"), Iₐₚₚ=rand(Normal(1.5, 0.05))) for i in 1:NE] -inhi = [AdamINP(name=Symbol("INP$i"), Iₐₚₚ=rand(Normal(0.1, 0.05))) for i in 1:NI] # bump up to 0.3 - -g = MetaDiGraph() - -for ne ∈ exci - for ni ∈ inhi - make_nmda_edge!(g, ne, ni) - end -end - -for ne ∈ exci - for ne ∈ exci - add_edge!(g, ne => ne; weight=1.0) - end -end - -for ni ∈ inhi - for ne ∈ exci[1:20] - add_edge!(g, ni => ne; weight=ḡᵢ/NI) - end -end - -tspan = (0.0, 500.0) -@time @named sys = system_from_graph(g, graphdynamics=true) -@time prob = ODEProblem(sys, [], tspan) -@time sol = solve(prob, Tsit5(), saveat=0.5) - ## Testing Glu for threshold setting ## Commented out for now but useful for tuning later so leaving in the file # exci = AdamPYR(name=:PYR, Iₐₚₚ=0.25) diff --git a/src/blox/adams_graphdynamics_draft.jl b/src/blox/adams_graphdynamics_draft.jl new file mode 100644 index 00000000..e29ce2b2 --- /dev/null +++ b/src/blox/adams_graphdynamics_draft.jl @@ -0,0 +1,88 @@ +using Neuroblox +using GraphDynamics + +using Neuroblox.GraphDynamicsInterop +GraphDynamicsInterop.issupported(::AdamNMDAR) = true +GraphDynamicsInterop.components(v::AdamNMDAR) = (v,) + +function GraphDynamicsInterop.to_subsystem(v::AdamNMDAR) + # Extract default parameter values + k_on = GraphDynamicsInterop.recursive_getdefault(v.system.k_on) + k_off = GraphDynamicsInterop.recursive_getdefault(v.system.k_off) + k_r = GraphDynamicsInterop.recursive_getdefault(v.system.k_r) + k_d = GraphDynamicsInterop.recursive_getdefault(v.system.k_d) + k_unblock = GraphDynamicsInterop.recursive_getdefault(v.system.k_unblock) + k_block = GraphDynamicsInterop.recursive_getdefault(v.system.k_block) + α = GraphDynamicsInterop.recursive_getdefault(v.system.α) + β = GraphDynamicsInterop.recursive_getdefault(v.system.β) + + params = SubsystemParams{AdamNMDAR}(; k_on, k_off, k_r, k_d, k_unblock, k_block, α, β) + + # Extract the default values of states + C = GraphDynamicsInterop.recursive_getdefault(v.system.C) + C_A = GraphDynamicsInterop.recursive_getdefault(v.system.C_A) + C_AA = GraphDynamicsInterop.recursive_getdefault(v.system.C_AA) + D_AA = GraphDynamicsInterop.recursive_getdefault(v.system.D_AA) + O_AA = GraphDynamicsInterop.recursive_getdefault(v.system.O_AA) + O_AAB = GraphDynamicsInterop.recursive_getdefault(v.system.O_AAB) + C_AAB = GraphDynamicsInterop.recursive_getdefault(v.system.C_AAB) + D_AAB = GraphDynamicsInterop.recursive_getdefault(v.system.D_AAB) + C_AB = GraphDynamicsInterop.recursive_getdefault(v.system.C_AB) + C_B = GraphDynamicsInterop.recursive_getdefault(v.system.C_B) + states = SubsystemStates{AdamNMDAR}(; C, C_A, C_AA, D_AA, O_AA, O_AAB, C_AAB, D_AAB, C_AB, C_B) + + Subsystem(states, params) +end + +GraphDynamics.initialize_input(s::Subsystem{AdamNMDAR}) = (; Glu = 0.0, V = 0.0) + +function GraphDynamics.subsystem_differential(s::Subsystem{AdamNMDAR}, inputs, t) + # Unpack + (; Glu, V) = inputs + (; C, C_A, C_AA, D_AA, O_AA, O_AAB, C_AAB, D_AAB, C_AB, C_B) = s + (; k_on, k_off, k_r, k_d, k_unblock, k_block, α, β) = s + return SubsystemStates{AdamNMDAR}( + #=d/dt=# C = k_off*C_A - 2*k_on*Glu*C, + #=d/dt=# C_A = 2*k_off*C_AA + 2*k_on*Glu*C - (k_on*Glu + k_off)*C_A, + #=d/dt=# C_AA = k_on*Glu*C_A + α*O_AA + k_r*D_AA - (2*k_off + β + k_d)*C_AA, + #=d/dt=# D_AA = k_d*C_AA - k_r*D_AA, + #=d/dt=# O_AA = k_r*D_AA - α*O_AA, + #=d/dt=# O_AAB = k_unblock*C_AAB - k_block*O_AAB, + #=d/dt=# C_AAB = k_block*O_AAB - k_unblock*C_AAB, + #=d/dt=# D_AAB = k_d*C_AAB - k_r*D_AAB, + #=d/dt=# C_AB = k_off*C_AAB - 2*k_on*Glu*C_AB, + #=d/dt=# C_B = k_off*C_AB - 2*k_on*Glu*C_B + ) +end + +function (c::GraphDynamicsInterop.BasicConnection)(sys_src::Subsystem{AdamGlu}, sys_dst::Subsystem{AdamNMDAR}, t) + w = c.weight + Glu = sys_src.Glu + V = 0.0 + (; Glu, V) +end + +function (c::GraphDynamicsInterop.BasicConnection)(sys_src::Subsystem{<:Neuroblox.AbstractAdamNeuron}, sys_dst::Subsystem{AdamNMDAR}, t) + w = c.weight + Glu = 0.0 + V = sys_src.V + (; Glu, V) +end + +function (c::GraphDynamicsInterop.BasicConnection)(sys_src::Subsystem{AdamNMDAR}, sys_dst::Subsystem{<:Neuroblox.AbstractAdamNeuron}, t) + w = c.weight + O_AA = sys_src.O_AA + jcn = w*O_AA*sys_dst.V + (; jcn) +end + +## Not GraphDynamics but work on connections +# Helper function for connections +function make_nmda_edge!(g, prenrn, postnrn) + glu = AdamGlu(name=Symbol("Glu$(prenrn.name)_$(postnrn.name)")) + nmda = AdamNMDAR(name=Symbol("NMDA$(prenrn.name)_$(postnrn.name)")) + add_edge!(g, prenrn => glu; weight=1.0) + add_edge!(g, glu => nmda; weight=1.0) + add_edge!(g, postnrn => nmda; weight=1.0) + add_edge!(g, nmda => postnrn; weight=1.0) +end \ No newline at end of file