From f458f6e8f306bef6604c8e9c1f8a5b45b7f9a188 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 30 Aug 2024 16:36:23 -0400 Subject: [PATCH 01/30] Add Izhikevich NGNMM Still need connection rule and will create others with different defaults for the other populations in the Adams paper --- src/Neuroblox.jl | 1 + src/blox/neural_mass.jl | 34 ++++++++++++++++++++++++++++++++++ 2 files changed, 35 insertions(+) diff --git a/src/Neuroblox.jl b/src/Neuroblox.jl index 0f0495bf..1d744b91 100644 --- a/src/Neuroblox.jl +++ b/src/Neuroblox.jl @@ -219,5 +219,6 @@ export run_experiment!, run_trial! export addnontunableparams export get_weights, get_dynamic_states, get_idx_tagged_vars, get_eqidx_tagged_vars export BalloonModel,LeadField, boldsignal_endo_balloon +export PYR_Izh end diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index 5de98ab9..b1b71ebf 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -502,3 +502,37 @@ struct KuramotoOscillator <: NeuralMassBlox end end end + +struct PYR_Izh <: NeuralMassBlox + params + output + jcn + odesystem + namespace + function PYR_Izh(; + name, + namespace=nothing, + Δ, + α=0.6215, + gₛ=1.2308, + η̄=0.12, + I_ext=0.0, + eᵣ=1.0, + a=0.0077, + b=-0.0062, + wⱼ=0.0189, + sⱼ=1.2308, + τₛ2.6, + κ=1.0) + p = paramscoping(Δ=Δ, α=α, gₛ=gₛ, η̄=η̄, I_ext=I_ext, eᵣ=eᵣ, a=a, b=b, wⱼ=wⱼ, sⱼ=sⱼ, κ=κ) + Δ, α, gₛ, η̄, I_ext, eᵣ, a, b, wⱼ, sⱼ, κ = p + sts = @variables r(t)=0.0 v(t)=0.0 w(t)=0.0 s(t)=0.0 [output=true] jcn(t)=0.0 [input=true] + eqs = [ D(r) ~ Δ/π + 2*r*v - (α+gₛ*s)*r, + D(v) ~ v^2 - α*v - w + η̄ + I_ext + gₛ*s*κ*(eᵣ - v) + jcn - (π*r)^2, + D(w) ~ a*(b*v - w) + wⱼ*r, + D(s) ~ -s/τₛ + sⱼ*r + ] + sys = System(eqs, t, sts, p; name=name) + new(p, sts[4], sts[5], sys, namespace) + end +end \ No newline at end of file From 354e228d1c1163747552a77938553ad41fcda5e7 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 30 Aug 2024 16:37:15 -0400 Subject: [PATCH 02/30] Oops typo --- src/blox/neural_mass.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index b1b71ebf..b6e2cdac 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -522,7 +522,7 @@ struct PYR_Izh <: NeuralMassBlox b=-0.0062, wⱼ=0.0189, sⱼ=1.2308, - τₛ2.6, + τₛ=2.6, κ=1.0) p = paramscoping(Δ=Δ, α=α, gₛ=gₛ, η̄=η̄, I_ext=I_ext, eᵣ=eᵣ, a=a, b=b, wⱼ=wⱼ, sⱼ=sⱼ, κ=κ) Δ, α, gₛ, η̄, I_ext, eᵣ, a, b, wⱼ, sⱼ, κ = p From af0b51820139bf44e8be3c740a0f27d927a03550 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 30 Aug 2024 16:42:58 -0400 Subject: [PATCH 03/30] Need access to voltage for connections --- examples/adams_example_for_brain_r01.jl | 0 src/blox/neural_mass.jl | 3 ++- 2 files changed, 2 insertions(+), 1 deletion(-) create mode 100644 examples/adams_example_for_brain_r01.jl diff --git a/examples/adams_example_for_brain_r01.jl b/examples/adams_example_for_brain_r01.jl new file mode 100644 index 00000000..e69de29b diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index b6e2cdac..86a9db71 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -506,6 +506,7 @@ end struct PYR_Izh <: NeuralMassBlox params output + voltage jcn odesystem namespace @@ -533,6 +534,6 @@ struct PYR_Izh <: NeuralMassBlox D(s) ~ -s/τₛ + sⱼ*r ] sys = System(eqs, t, sts, p; name=name) - new(p, sts[4], sts[5], sys, namespace) + new(p, sts[4], sts[2], sts[5], sys, namespace) end end \ No newline at end of file From 40c6f77904a1b9f0b372726b3c103e48022af266 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 30 Aug 2024 16:46:10 -0400 Subject: [PATCH 04/30] NGNMM connections Done via firing rate -> voltage coupling --- src/blox/connections.jl | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/blox/connections.jl b/src/blox/connections.jl index c4be906c..b80018a0 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -965,3 +965,21 @@ function (bc::BloxConnector)( bc(stim, neuron; kwargs...) end end + +function (bc::BloxConnector)( + bloxout::PYR_Izh, + bloxin::PYR_Izh; + kwargs... +) + sys_out = get_namespaced_sys(bloxout) + sys_in = get_namespaced_sys(bloxin) + + w = generate_weight_param(bloxout, bloxin; kwargs...) + push!(bc.weights, w) + + s_presyn = namespace_expr(bloxout.output, sys_out) + v_postsyn = namespace_expr(bloxin.voltage, sys_in) + eq = sys_in.jcn ~ (1-sys_in.κ)*sys_out.gₛ*s_presyn*(sys_in.eᵣ-v_postsyn) + + accumulate_equation!(bc, eq) +end \ No newline at end of file From cda1724632d952c030a69c9b88ef6e7ab7e8c930 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 30 Aug 2024 17:03:20 -0400 Subject: [PATCH 05/30] Oops --- src/blox/neural_mass.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index 86a9db71..b86b8aa2 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -513,7 +513,7 @@ struct PYR_Izh <: NeuralMassBlox function PYR_Izh(; name, namespace=nothing, - Δ, + Δ=0.02, α=0.6215, gₛ=1.2308, η̄=0.12, From e52e1286ddae13958f6a6ea303b93f5d8f70c118 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 30 Aug 2024 17:07:56 -0400 Subject: [PATCH 06/30] Example to reproduce interacting neural masses Reproduces Figure 3 (right half) of the Chen/Campbell paper --- examples/adams_example_for_brain_r01.jl | 31 +++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/examples/adams_example_for_brain_r01.jl b/examples/adams_example_for_brain_r01.jl index e69de29b..d4592b65 100644 --- a/examples/adams_example_for_brain_r01.jl +++ b/examples/adams_example_for_brain_r01.jl @@ -0,0 +1,31 @@ +using Neuroblox +using DifferentialEquations +using DataFrames +using Test +using Distributions +using Statistics +using LinearAlgebra +using Graphs +using MetaGraphs +using Random + +@named popP = PYR_Izh(η̄=0.08, κ=0.8) +@named popQ = PYR_Izh(η̄=0.08, κ=0.2, wⱼ=0.0095, a=0.077) + +g = MetaDiGraph() +add_blox!.(Ref(g), [popP, popQ]) +add_edge!(g, popP => popQ; weight=1.0) #weight is acutaally meaningless here +add_edge!(g, popQ => popP; weight=1.0) #weight is acutaally meaningless here + +@named sys = system_from_graph(g) +sys = structural_simplify(sys) + +sim_dur = 800.0 +prob = ODEProblem(sys, [], (0.0, sim_dur)) +sol = solve(prob, Tsit5(), saveat=1.0) + +# Reproduce Chen/Campbell figure 3 panels on the right, especially 3f +using Plots +plot(sol, idxs=[1, 5]) +plot(sol, idxs=[2, 6]) +plot(sol, idxs=[3, 7]) \ No newline at end of file From 58d4f8b7f3c8f4f667dafa03a51d64eb2488a99d Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 31 Aug 2024 17:42:16 -0400 Subject: [PATCH 07/30] Weight shouldn't be meaningless --- 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 b80018a0..940cd6eb 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -979,7 +979,7 @@ function (bc::BloxConnector)( s_presyn = namespace_expr(bloxout.output, sys_out) v_postsyn = namespace_expr(bloxin.voltage, sys_in) - eq = sys_in.jcn ~ (1-sys_in.κ)*sys_out.gₛ*s_presyn*(sys_in.eᵣ-v_postsyn) + eq = sys_in.jcn ~ w*(1-sys_in.κ)*sys_out.gₛ*s_presyn*(sys_in.eᵣ-v_postsyn) accumulate_equation!(bc, eq) end \ No newline at end of file From 6ac7091b0dda3fd44ad21a2e409a85a5a0dac2d2 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 31 Aug 2024 20:43:02 -0400 Subject: [PATCH 08/30] Ugh PING I guess --- examples/adams_example_for_brain_r01.jl | 63 ++++++++++++++++++++++++- src/Neuroblox.jl | 2 +- src/blox/connections.jl | 3 +- src/blox/neural_mass.jl | 32 ++++++++++++- 4 files changed, 95 insertions(+), 5 deletions(-) diff --git a/examples/adams_example_for_brain_r01.jl b/examples/adams_example_for_brain_r01.jl index d4592b65..79cdacac 100644 --- a/examples/adams_example_for_brain_r01.jl +++ b/examples/adams_example_for_brain_r01.jl @@ -28,4 +28,65 @@ sol = solve(prob, Tsit5(), saveat=1.0) using Plots plot(sol, idxs=[1, 5]) plot(sol, idxs=[2, 6]) -plot(sol, idxs=[3, 7]) \ No newline at end of file +plot(sol, idxs=[3, 7]) + + + + +## DO NOT TOUCH LARGELY FAILURES THAT RUN sol_dde_with_delays +# I want to salvage a bit of this so leaving in for now +# -AGC + +# # Reproducing ketamine dynamics +# @named PYR = PYR_Izh(η̄=0.08, κ=0.5, I_ext=0.25) +# @named INP = PYR_Izh(η̄=0.08, κ=0.5, wⱼ=0.0095, a=0.077, I_ext=0.5) + +# g = MetaDiGraph() +# add_blox!.(Ref(g), [PYR, INP]) +# add_edge!(g, PYR => INP; weight=1.0) #weight is acutaally meaningless here +# add_edge!(g, INP => PYR; weight=1.0) #weight is acutaally meaningless here + +# @named sys = system_from_graph(g) +# sys = structural_simplify(sys) + +# sim_dur = 3000.0 +# prob = ODEProblem(sys, [], (0.0, sim_dur)) +# sol = solve(prob, Tsit5(), saveat=1.0) + +# # Reproduce Chen/Campbell figure 3 panels on the right, especially 3f +# using Plots +# plot(sol, idxs=[1, 5]) +# plot(sol, idxs=[2, 6]) +# plot(sol, idxs=[3, 7]) + +# # Reproducing ketamine dynamics +# kappa_mod = 0.9 +# I_adj = 0.1 +# ω=4.0 +# @named PYR = PYR_Izh(η̄=0.08, κ=kappa_mod, I_ext=I_adj, ω=ω*2*π/1000) +# @named INP = PYR_Izh(η̄=0.08, κ=1-kappa_mod, wⱼ=0.0095*15, a=0.077, τₛ=5.0, I_ext=I_adj, ω=ω*2*π/1000) + +# g = MetaDiGraph() +# add_blox!.(Ref(g), [PYR, INP]) +# add_edge!(g, PYR => INP; weight=1.0) +# add_edge!(g, INP => PYR; weight=1.0) + +# @named sys = system_from_graph(g) +# sys = structural_simplify(sys) + +# sim_dur = 1000000.0 +# prob = ODEProblem(sys, [], (0.0, sim_dur)) +# sol = solve(prob, Tsit5(), saveat=0.1) + +# # Reproduce Chen/Campbell figure 3 panels on the right, especially 3f +# #using Plots, DSP +# plot(sol, idxs=[1, 5]) +# plot(sol, idxs=[2, 6]) + +# data = Array(sol) +# data = data .+ rand(Normal(0, 0.1), size(data)) +# wp = welch_pgram(data[6, :]; fs=10000) +# plot(wp.freq, pow2db.(wp.power), xlim=(0, 50)) + +# spec = spectrogram(data[6, :], 200; fs=10000) +# heatmap(spec.time, spec.freq, pow2db.(spec.power), ylim=(0, 50), xlim=(0, 100)) \ No newline at end of file diff --git a/src/Neuroblox.jl b/src/Neuroblox.jl index 1d744b91..19950910 100644 --- a/src/Neuroblox.jl +++ b/src/Neuroblox.jl @@ -219,6 +219,6 @@ export run_experiment!, run_trial! export addnontunableparams export get_weights, get_dynamic_states, get_idx_tagged_vars, get_eqidx_tagged_vars export BalloonModel,LeadField, boldsignal_endo_balloon -export PYR_Izh +export PYR_Izh, QIF_PING_NGNMM end diff --git a/src/blox/connections.jl b/src/blox/connections.jl index 940cd6eb..5e1ef44f 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -982,4 +982,5 @@ function (bc::BloxConnector)( eq = sys_in.jcn ~ w*(1-sys_in.κ)*sys_out.gₛ*s_presyn*(sys_in.eᵣ-v_postsyn) accumulate_equation!(bc, eq) -end \ No newline at end of file +end + diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index b86b8aa2..db34afe3 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -524,16 +524,44 @@ struct PYR_Izh <: NeuralMassBlox wⱼ=0.0189, sⱼ=1.2308, τₛ=2.6, - κ=1.0) + κ=1.0, + ω=0.0) p = paramscoping(Δ=Δ, α=α, gₛ=gₛ, η̄=η̄, I_ext=I_ext, eᵣ=eᵣ, a=a, b=b, wⱼ=wⱼ, sⱼ=sⱼ, κ=κ) Δ, α, gₛ, η̄, I_ext, eᵣ, a, b, wⱼ, sⱼ, κ = p sts = @variables r(t)=0.0 v(t)=0.0 w(t)=0.0 s(t)=0.0 [output=true] jcn(t)=0.0 [input=true] eqs = [ D(r) ~ Δ/π + 2*r*v - (α+gₛ*s)*r, - D(v) ~ v^2 - α*v - w + η̄ + I_ext + gₛ*s*κ*(eᵣ - v) + jcn - (π*r)^2, + D(v) ~ v^2 - α*v - w + η̄ + I_ext*sin(ω*t) + gₛ*s*κ*(eᵣ - v) + jcn - (π*r)^2, D(w) ~ a*(b*v - w) + wⱼ*r, D(s) ~ -s/τₛ + sⱼ*r ] sys = System(eqs, t, sts, p; name=name) new(p, sts[4], sts[2], sts[5], sys, namespace) end +end + +struct QIF_PING_NGNMM <: NeuralMassBlox + params + output + voltage + jcn + odesystem + namespace + function QIF_PING_NGNMM(; + name, + namespace=noathing + Δ=0.02, + τₘ, + H, + I_ext, + J_internal, + A) + p = paramscoping(Δ=Δ, τₘ=τₘ, H=H, I_ext=I_ext, J_internal=J_internal, A=A) + Δ, τₘ, H, I_ext, J_internal, A = p + sts = @variables r(t)=0.0 [output=true] v(t)=0.0 jcn(t)=0.0 [input=true] + @brownian ξ + eqs = [D(r) ~ Δ/(π*τₘ^2) + 2*r*v/τₘ, + D(v) ~ (v^2 + H + I_ext)/τₘ - τₘ*(π*r)^2 + J_internal*r + A*ξ + jcn] + sys = System(eqs, t, sts, p; name=name) + new(p, sts[1], sts[2], sts[3], sys, namespace) + end end \ No newline at end of file From 5aa2e053cd99c9d44be2ad5c8df830f2756f143a Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 31 Aug 2024 20:43:03 -0400 Subject: [PATCH 09/30] Create saving_voltage_utils.jl --- saving_voltage_utils.jl | 46 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 saving_voltage_utils.jl diff --git a/saving_voltage_utils.jl b/saving_voltage_utils.jl new file mode 100644 index 00000000..af7ce44e --- /dev/null +++ b/saving_voltage_utils.jl @@ -0,0 +1,46 @@ +using FileIO +using DataFrames + +# Loading single file - will loop to load all eventually +# save as csv +data = load("/Users/achesebro/Downloads/sim0003.jld2") +df = data["df"] +all_voltages = select(df, r"₊V") +time_vec = select(df, :timestamp)[:, 1] + +# Get all region names + +# Make into utility function +# Create test (eventually) +all_names = names(all_voltages) +unique_names = Vector{String}() + +for i in eachindex(all_names) + name_temp = all_names[i] + name = split(name_temp, "₊")[1] #check for innernamespaceof for global namespace removal + push!(unique_names, name) +end + +unique_names = unique(unique_names) + +# Get unique time indices (not originally unique because of the callbacks) +time_indices = Vector{Int}() +for i in eachindex(time_vec) + if i == 1 || time_vec[i] != time_vec[i-1] + push!(time_indices, i) + end +end + +# Extract the voltages at the unique time indices +regions = zeros(length(time_indices), length(unique_names)) + +for i in eachindex(unique_names) + name = unique_names[i] + region = select(all_voltages, Regex(name)) + temp = reduce(+, eachcol(region)) ./ ncol(region) + println(size(region)) + regions[:, i] = temp[time_indices] +end + +using Plots +plot(time_vec[time_indices], regions, xlabel="Time", ylabel="Voltage", title="Average Voltage in Each Region") \ No newline at end of file From c2dc922ef76f07105de8447abce9928a56286ead Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 31 Aug 2024 20:44:58 -0400 Subject: [PATCH 10/30] QIF PING connectors --- src/blox/connections.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/blox/connections.jl b/src/blox/connections.jl index 5e1ef44f..5b23f569 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -984,3 +984,19 @@ function (bc::BloxConnector)( accumulate_equation!(bc, eq) end +function (bc::BloxConnector)( + bloxout::QIF_PING_NGNMM, + bloxin::QIF_PING_NGNMM; + kwargs... +) + sys_out = get_namespaced_sys(bloxout) + sys_in = get_namespaced_sys(bloxin) + + w = generate_weight_param(bloxout, bloxin; kwargs...) + push!(bc.weights, w) + + x = namespace_expr(bloxout.output, sys_out) + eq = sys_in.jcn ~ w*x + + accumulate_equation!(bc, eq) +end \ No newline at end of file From adce0262b18aa1561b2f168434b2f913e40e7474 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 31 Aug 2024 20:51:33 -0400 Subject: [PATCH 11/30] Remove noise for now Doesn't look like it's necessary so remove to speed up sims --- examples/adams_example_for_brain_r01.jl | 6 +++++- src/blox/neural_mass.jl | 20 ++++++++++---------- 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/examples/adams_example_for_brain_r01.jl b/examples/adams_example_for_brain_r01.jl index 79cdacac..16799ead 100644 --- a/examples/adams_example_for_brain_r01.jl +++ b/examples/adams_example_for_brain_r01.jl @@ -9,6 +9,11 @@ using Graphs using MetaGraphs using Random +# PING QIF theta-nested gamma oscillations + + + +# Chen/Campbell populations - limited utility for now @named popP = PYR_Izh(η̄=0.08, κ=0.8) @named popQ = PYR_Izh(η̄=0.08, κ=0.2, wⱼ=0.0095, a=0.077) @@ -32,7 +37,6 @@ plot(sol, idxs=[3, 7]) - ## DO NOT TOUCH LARGELY FAILURES THAT RUN sol_dde_with_delays # I want to salvage a bit of this so leaving in for now # -AGC diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index db34afe3..5c332e06 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -549,18 +549,18 @@ struct QIF_PING_NGNMM <: NeuralMassBlox function QIF_PING_NGNMM(; name, namespace=noathing - Δ=0.02, - τₘ, - H, - I_ext, - J_internal, - A) - p = paramscoping(Δ=Δ, τₘ=τₘ, H=H, I_ext=I_ext, J_internal=J_internal, A=A) - Δ, τₘ, H, I_ext, J_internal, A = p + Δ=1.0, + τₘ=20.0, + H=1.3, + I_ext=0.0, + ω=0.0, + J_internal=8.0) + p = paramscoping(Δ=Δ, τₘ=τₘ, H=H, I_ext=I_ext, J_internal=J_internal) + Δ, τₘ, H, I_ext, J_internal = p sts = @variables r(t)=0.0 [output=true] v(t)=0.0 jcn(t)=0.0 [input=true] - @brownian ξ + #@brownian ξ eqs = [D(r) ~ Δ/(π*τₘ^2) + 2*r*v/τₘ, - D(v) ~ (v^2 + H + I_ext)/τₘ - τₘ*(π*r)^2 + J_internal*r + A*ξ + jcn] + D(v) ~ (v^2 + H + I_ext*sin(ω*t))/τₘ - τₘ*(π*r)^2 + J_internal*r + jcn] sys = System(eqs, t, sts, p; name=name) new(p, sts[1], sts[2], sts[3], sys, namespace) end From 484ff72a74596cb22888671be4d99945a8825673 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 31 Aug 2024 20:51:59 -0400 Subject: [PATCH 12/30] Oops typo --- src/blox/neural_mass.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index 5c332e06..0acbcf41 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -548,7 +548,7 @@ struct QIF_PING_NGNMM <: NeuralMassBlox namespace function QIF_PING_NGNMM(; name, - namespace=noathing + namespace=nothing, Δ=1.0, τₘ=20.0, H=1.3, From b6b70fc889962f43b58cae15484f7cf13b08354f Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 31 Aug 2024 22:04:30 -0400 Subject: [PATCH 13/30] Ok it kinda works --- examples/adams_example_for_brain_r01.jl | 23 +++++++++++++++++++++++ src/blox/neural_mass.jl | 7 ++++--- 2 files changed, 27 insertions(+), 3 deletions(-) diff --git a/examples/adams_example_for_brain_r01.jl b/examples/adams_example_for_brain_r01.jl index 16799ead..a7170373 100644 --- a/examples/adams_example_for_brain_r01.jl +++ b/examples/adams_example_for_brain_r01.jl @@ -10,6 +10,29 @@ using MetaGraphs using Random # PING QIF theta-nested gamma oscillations +@named exci_PING = QIF_PING_NGNMM(I_ext=10.0, ω=5*2*π/1000, J_internal=8.0, H=1.3, Δ=1.0, τₘ=20.0, A=0.2) +@named inhi_PING = QIF_PING_NGNMM(I_ext=5.0, ω=5*2*π/1000, J_internal=0.0, H=-5.0, Δ=1.0, τₘ=10.0, A=0.0) + +g = MetaDiGraph() +add_blox!.(Ref(g), [exci_PING, inhi_PING]) +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) + +sim_dur = 1000.0 +prob = SDEProblem(sys, [], (0.0, sim_dur)) +sol = solve(prob, RKMil(), saveat=0.1) + +using Plots, DSP, CSV +plot(sol, idxs=[1, 3]) +plot(sol, idxs=[2, 4]) + +data = Array(sol) +#data = data .+ rand(Normal(0, 10), size(data)) +wp = welch_pgram(data[2, :]; fs=2000) +plot(wp.freq, pow2db.(wp.power), xlim=(0, 50), xlabel="Frequency (Hz)" , ylabel="Power (dB)", title="Power Spectrum of Excitatory Population") diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index 0acbcf41..ab988f49 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -554,13 +554,14 @@ struct QIF_PING_NGNMM <: NeuralMassBlox H=1.3, I_ext=0.0, ω=0.0, - J_internal=8.0) + J_internal=8.0, + A=0.0) p = paramscoping(Δ=Δ, τₘ=τₘ, H=H, I_ext=I_ext, J_internal=J_internal) Δ, τₘ, H, I_ext, J_internal = p sts = @variables r(t)=0.0 [output=true] v(t)=0.0 jcn(t)=0.0 [input=true] - #@brownian ξ + @brownian ξ eqs = [D(r) ~ Δ/(π*τₘ^2) + 2*r*v/τₘ, - D(v) ~ (v^2 + H + I_ext*sin(ω*t))/τₘ - τₘ*(π*r)^2 + J_internal*r + jcn] + D(v) ~ (v^2 + H + I_ext*sin(ω*t))/τₘ - τₘ*(π*r)^2 + J_internal*r + A*ξ + jcn] sys = System(eqs, t, sts, p; name=name) new(p, sts[1], sts[2], sts[3], sys, namespace) end From b5bc169a176c1e1ec24caacebe78a886be0f2605 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 30 Aug 2024 16:36:23 -0400 Subject: [PATCH 14/30] Add Izhikevich NGNMM Still need connection rule and will create others with different defaults for the other populations in the Adams paper --- src/Neuroblox.jl | 1 + src/blox/neural_mass.jl | 34 ++++++++++++++++++++++++++++++++++ 2 files changed, 35 insertions(+) diff --git a/src/Neuroblox.jl b/src/Neuroblox.jl index f8457f8e..d3789d11 100644 --- a/src/Neuroblox.jl +++ b/src/Neuroblox.jl @@ -237,6 +237,7 @@ export run_experiment!, run_trial! export addnontunableparams export get_weights, get_dynamic_states, get_idx_tagged_vars, get_eqidx_tagged_vars export BalloonModel,LeadField, boldsignal_endo_balloon +export PYR_Izh export meanfield, meanfield!, rasterplot, rasterplot!, stackplot, stackplot!, voltage_stack end diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index 3834c5e4..c6735c2b 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -502,3 +502,37 @@ struct KuramotoOscillator <: NeuralMassBlox end end end + +struct PYR_Izh <: NeuralMassBlox + params + output + jcn + odesystem + namespace + function PYR_Izh(; + name, + namespace=nothing, + Δ, + α=0.6215, + gₛ=1.2308, + η̄=0.12, + I_ext=0.0, + eᵣ=1.0, + a=0.0077, + b=-0.0062, + wⱼ=0.0189, + sⱼ=1.2308, + τₛ2.6, + κ=1.0) + p = paramscoping(Δ=Δ, α=α, gₛ=gₛ, η̄=η̄, I_ext=I_ext, eᵣ=eᵣ, a=a, b=b, wⱼ=wⱼ, sⱼ=sⱼ, κ=κ) + Δ, α, gₛ, η̄, I_ext, eᵣ, a, b, wⱼ, sⱼ, κ = p + sts = @variables r(t)=0.0 v(t)=0.0 w(t)=0.0 s(t)=0.0 [output=true] jcn(t)=0.0 [input=true] + eqs = [ D(r) ~ Δ/π + 2*r*v - (α+gₛ*s)*r, + D(v) ~ v^2 - α*v - w + η̄ + I_ext + gₛ*s*κ*(eᵣ - v) + jcn - (π*r)^2, + D(w) ~ a*(b*v - w) + wⱼ*r, + D(s) ~ -s/τₛ + sⱼ*r + ] + sys = System(eqs, t, sts, p; name=name) + new(p, sts[4], sts[5], sys, namespace) + end +end \ No newline at end of file From ca99a5f86fb662e1a09e7d606937f524a3ee797f Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 30 Aug 2024 16:37:15 -0400 Subject: [PATCH 15/30] Oops typo --- src/blox/neural_mass.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index c6735c2b..103e6d2f 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -522,7 +522,7 @@ struct PYR_Izh <: NeuralMassBlox b=-0.0062, wⱼ=0.0189, sⱼ=1.2308, - τₛ2.6, + τₛ=2.6, κ=1.0) p = paramscoping(Δ=Δ, α=α, gₛ=gₛ, η̄=η̄, I_ext=I_ext, eᵣ=eᵣ, a=a, b=b, wⱼ=wⱼ, sⱼ=sⱼ, κ=κ) Δ, α, gₛ, η̄, I_ext, eᵣ, a, b, wⱼ, sⱼ, κ = p From 4ffe1680c18d5560812d800380eb6dd02856b3c6 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 30 Aug 2024 16:42:58 -0400 Subject: [PATCH 16/30] Need access to voltage for connections --- examples/adams_example_for_brain_r01.jl | 0 src/blox/neural_mass.jl | 3 ++- 2 files changed, 2 insertions(+), 1 deletion(-) create mode 100644 examples/adams_example_for_brain_r01.jl diff --git a/examples/adams_example_for_brain_r01.jl b/examples/adams_example_for_brain_r01.jl new file mode 100644 index 00000000..e69de29b diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index 103e6d2f..eabf40be 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -506,6 +506,7 @@ end struct PYR_Izh <: NeuralMassBlox params output + voltage jcn odesystem namespace @@ -533,6 +534,6 @@ struct PYR_Izh <: NeuralMassBlox D(s) ~ -s/τₛ + sⱼ*r ] sys = System(eqs, t, sts, p; name=name) - new(p, sts[4], sts[5], sys, namespace) + new(p, sts[4], sts[2], sts[5], sys, namespace) end end \ No newline at end of file From 45143e1c35a59763cdb3e2f518e8030d0d923027 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 30 Aug 2024 16:46:10 -0400 Subject: [PATCH 17/30] NGNMM connections Done via firing rate -> voltage coupling --- src/blox/connections.jl | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/blox/connections.jl b/src/blox/connections.jl index 6a99f0ef..67c66b62 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -959,3 +959,21 @@ function (bc::BloxConnector)( bc(stim, neuron; kwargs...) end end + +function (bc::BloxConnector)( + bloxout::PYR_Izh, + bloxin::PYR_Izh; + kwargs... +) + sys_out = get_namespaced_sys(bloxout) + sys_in = get_namespaced_sys(bloxin) + + w = generate_weight_param(bloxout, bloxin; kwargs...) + push!(bc.weights, w) + + s_presyn = namespace_expr(bloxout.output, sys_out) + v_postsyn = namespace_expr(bloxin.voltage, sys_in) + eq = sys_in.jcn ~ (1-sys_in.κ)*sys_out.gₛ*s_presyn*(sys_in.eᵣ-v_postsyn) + + accumulate_equation!(bc, eq) +end \ No newline at end of file From a0a749e2491d193798df22cde8b401cca1133e8b Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 30 Aug 2024 17:03:20 -0400 Subject: [PATCH 18/30] Oops --- src/blox/neural_mass.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index eabf40be..f097e596 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -513,7 +513,7 @@ struct PYR_Izh <: NeuralMassBlox function PYR_Izh(; name, namespace=nothing, - Δ, + Δ=0.02, α=0.6215, gₛ=1.2308, η̄=0.12, From 68f79fb3753b02d06350d86f1235e2ba38108a27 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 30 Aug 2024 17:07:56 -0400 Subject: [PATCH 19/30] Example to reproduce interacting neural masses Reproduces Figure 3 (right half) of the Chen/Campbell paper --- examples/adams_example_for_brain_r01.jl | 31 +++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/examples/adams_example_for_brain_r01.jl b/examples/adams_example_for_brain_r01.jl index e69de29b..d4592b65 100644 --- a/examples/adams_example_for_brain_r01.jl +++ b/examples/adams_example_for_brain_r01.jl @@ -0,0 +1,31 @@ +using Neuroblox +using DifferentialEquations +using DataFrames +using Test +using Distributions +using Statistics +using LinearAlgebra +using Graphs +using MetaGraphs +using Random + +@named popP = PYR_Izh(η̄=0.08, κ=0.8) +@named popQ = PYR_Izh(η̄=0.08, κ=0.2, wⱼ=0.0095, a=0.077) + +g = MetaDiGraph() +add_blox!.(Ref(g), [popP, popQ]) +add_edge!(g, popP => popQ; weight=1.0) #weight is acutaally meaningless here +add_edge!(g, popQ => popP; weight=1.0) #weight is acutaally meaningless here + +@named sys = system_from_graph(g) +sys = structural_simplify(sys) + +sim_dur = 800.0 +prob = ODEProblem(sys, [], (0.0, sim_dur)) +sol = solve(prob, Tsit5(), saveat=1.0) + +# Reproduce Chen/Campbell figure 3 panels on the right, especially 3f +using Plots +plot(sol, idxs=[1, 5]) +plot(sol, idxs=[2, 6]) +plot(sol, idxs=[3, 7]) \ No newline at end of file From 6ee74eba0e4a72c56eb5176cd6641a7b97467da8 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 31 Aug 2024 17:42:16 -0400 Subject: [PATCH 20/30] Weight shouldn't be meaningless --- 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 67c66b62..1dc1fc80 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -973,7 +973,7 @@ function (bc::BloxConnector)( s_presyn = namespace_expr(bloxout.output, sys_out) v_postsyn = namespace_expr(bloxin.voltage, sys_in) - eq = sys_in.jcn ~ (1-sys_in.κ)*sys_out.gₛ*s_presyn*(sys_in.eᵣ-v_postsyn) + eq = sys_in.jcn ~ w*(1-sys_in.κ)*sys_out.gₛ*s_presyn*(sys_in.eᵣ-v_postsyn) accumulate_equation!(bc, eq) end \ No newline at end of file From 7ab33ab6630cc8dfd6b2af8b68c190279f424b54 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 31 Aug 2024 20:43:02 -0400 Subject: [PATCH 21/30] Ugh PING I guess --- examples/adams_example_for_brain_r01.jl | 63 ++++++++++++++++++++++++- src/Neuroblox.jl | 2 +- src/blox/connections.jl | 3 +- src/blox/neural_mass.jl | 32 ++++++++++++- 4 files changed, 95 insertions(+), 5 deletions(-) diff --git a/examples/adams_example_for_brain_r01.jl b/examples/adams_example_for_brain_r01.jl index d4592b65..79cdacac 100644 --- a/examples/adams_example_for_brain_r01.jl +++ b/examples/adams_example_for_brain_r01.jl @@ -28,4 +28,65 @@ sol = solve(prob, Tsit5(), saveat=1.0) using Plots plot(sol, idxs=[1, 5]) plot(sol, idxs=[2, 6]) -plot(sol, idxs=[3, 7]) \ No newline at end of file +plot(sol, idxs=[3, 7]) + + + + +## DO NOT TOUCH LARGELY FAILURES THAT RUN sol_dde_with_delays +# I want to salvage a bit of this so leaving in for now +# -AGC + +# # Reproducing ketamine dynamics +# @named PYR = PYR_Izh(η̄=0.08, κ=0.5, I_ext=0.25) +# @named INP = PYR_Izh(η̄=0.08, κ=0.5, wⱼ=0.0095, a=0.077, I_ext=0.5) + +# g = MetaDiGraph() +# add_blox!.(Ref(g), [PYR, INP]) +# add_edge!(g, PYR => INP; weight=1.0) #weight is acutaally meaningless here +# add_edge!(g, INP => PYR; weight=1.0) #weight is acutaally meaningless here + +# @named sys = system_from_graph(g) +# sys = structural_simplify(sys) + +# sim_dur = 3000.0 +# prob = ODEProblem(sys, [], (0.0, sim_dur)) +# sol = solve(prob, Tsit5(), saveat=1.0) + +# # Reproduce Chen/Campbell figure 3 panels on the right, especially 3f +# using Plots +# plot(sol, idxs=[1, 5]) +# plot(sol, idxs=[2, 6]) +# plot(sol, idxs=[3, 7]) + +# # Reproducing ketamine dynamics +# kappa_mod = 0.9 +# I_adj = 0.1 +# ω=4.0 +# @named PYR = PYR_Izh(η̄=0.08, κ=kappa_mod, I_ext=I_adj, ω=ω*2*π/1000) +# @named INP = PYR_Izh(η̄=0.08, κ=1-kappa_mod, wⱼ=0.0095*15, a=0.077, τₛ=5.0, I_ext=I_adj, ω=ω*2*π/1000) + +# g = MetaDiGraph() +# add_blox!.(Ref(g), [PYR, INP]) +# add_edge!(g, PYR => INP; weight=1.0) +# add_edge!(g, INP => PYR; weight=1.0) + +# @named sys = system_from_graph(g) +# sys = structural_simplify(sys) + +# sim_dur = 1000000.0 +# prob = ODEProblem(sys, [], (0.0, sim_dur)) +# sol = solve(prob, Tsit5(), saveat=0.1) + +# # Reproduce Chen/Campbell figure 3 panels on the right, especially 3f +# #using Plots, DSP +# plot(sol, idxs=[1, 5]) +# plot(sol, idxs=[2, 6]) + +# data = Array(sol) +# data = data .+ rand(Normal(0, 0.1), size(data)) +# wp = welch_pgram(data[6, :]; fs=10000) +# plot(wp.freq, pow2db.(wp.power), xlim=(0, 50)) + +# spec = spectrogram(data[6, :], 200; fs=10000) +# heatmap(spec.time, spec.freq, pow2db.(spec.power), ylim=(0, 50), xlim=(0, 100)) \ No newline at end of file diff --git a/src/Neuroblox.jl b/src/Neuroblox.jl index d3789d11..fae53a72 100644 --- a/src/Neuroblox.jl +++ b/src/Neuroblox.jl @@ -237,7 +237,7 @@ export run_experiment!, run_trial! export addnontunableparams export get_weights, get_dynamic_states, get_idx_tagged_vars, get_eqidx_tagged_vars export BalloonModel,LeadField, boldsignal_endo_balloon -export PYR_Izh +export PYR_Izh, QIF_PING_NGNMM export meanfield, meanfield!, rasterplot, rasterplot!, stackplot, stackplot!, voltage_stack end diff --git a/src/blox/connections.jl b/src/blox/connections.jl index 1dc1fc80..cbc0ceb2 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -976,4 +976,5 @@ function (bc::BloxConnector)( eq = sys_in.jcn ~ w*(1-sys_in.κ)*sys_out.gₛ*s_presyn*(sys_in.eᵣ-v_postsyn) accumulate_equation!(bc, eq) -end \ No newline at end of file +end + diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index f097e596..e81e5717 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -524,16 +524,44 @@ struct PYR_Izh <: NeuralMassBlox wⱼ=0.0189, sⱼ=1.2308, τₛ=2.6, - κ=1.0) + κ=1.0, + ω=0.0) p = paramscoping(Δ=Δ, α=α, gₛ=gₛ, η̄=η̄, I_ext=I_ext, eᵣ=eᵣ, a=a, b=b, wⱼ=wⱼ, sⱼ=sⱼ, κ=κ) Δ, α, gₛ, η̄, I_ext, eᵣ, a, b, wⱼ, sⱼ, κ = p sts = @variables r(t)=0.0 v(t)=0.0 w(t)=0.0 s(t)=0.0 [output=true] jcn(t)=0.0 [input=true] eqs = [ D(r) ~ Δ/π + 2*r*v - (α+gₛ*s)*r, - D(v) ~ v^2 - α*v - w + η̄ + I_ext + gₛ*s*κ*(eᵣ - v) + jcn - (π*r)^2, + D(v) ~ v^2 - α*v - w + η̄ + I_ext*sin(ω*t) + gₛ*s*κ*(eᵣ - v) + jcn - (π*r)^2, D(w) ~ a*(b*v - w) + wⱼ*r, D(s) ~ -s/τₛ + sⱼ*r ] sys = System(eqs, t, sts, p; name=name) new(p, sts[4], sts[2], sts[5], sys, namespace) end +end + +struct QIF_PING_NGNMM <: NeuralMassBlox + params + output + voltage + jcn + odesystem + namespace + function QIF_PING_NGNMM(; + name, + namespace=noathing + Δ=0.02, + τₘ, + H, + I_ext, + J_internal, + A) + p = paramscoping(Δ=Δ, τₘ=τₘ, H=H, I_ext=I_ext, J_internal=J_internal, A=A) + Δ, τₘ, H, I_ext, J_internal, A = p + sts = @variables r(t)=0.0 [output=true] v(t)=0.0 jcn(t)=0.0 [input=true] + @brownian ξ + eqs = [D(r) ~ Δ/(π*τₘ^2) + 2*r*v/τₘ, + D(v) ~ (v^2 + H + I_ext)/τₘ - τₘ*(π*r)^2 + J_internal*r + A*ξ + jcn] + sys = System(eqs, t, sts, p; name=name) + new(p, sts[1], sts[2], sts[3], sys, namespace) + end end \ No newline at end of file From 82ad5daf38a9ae71434ee6bdd4c3811040c3c193 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 31 Aug 2024 20:43:03 -0400 Subject: [PATCH 22/30] Create saving_voltage_utils.jl --- saving_voltage_utils.jl | 46 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 saving_voltage_utils.jl diff --git a/saving_voltage_utils.jl b/saving_voltage_utils.jl new file mode 100644 index 00000000..af7ce44e --- /dev/null +++ b/saving_voltage_utils.jl @@ -0,0 +1,46 @@ +using FileIO +using DataFrames + +# Loading single file - will loop to load all eventually +# save as csv +data = load("/Users/achesebro/Downloads/sim0003.jld2") +df = data["df"] +all_voltages = select(df, r"₊V") +time_vec = select(df, :timestamp)[:, 1] + +# Get all region names + +# Make into utility function +# Create test (eventually) +all_names = names(all_voltages) +unique_names = Vector{String}() + +for i in eachindex(all_names) + name_temp = all_names[i] + name = split(name_temp, "₊")[1] #check for innernamespaceof for global namespace removal + push!(unique_names, name) +end + +unique_names = unique(unique_names) + +# Get unique time indices (not originally unique because of the callbacks) +time_indices = Vector{Int}() +for i in eachindex(time_vec) + if i == 1 || time_vec[i] != time_vec[i-1] + push!(time_indices, i) + end +end + +# Extract the voltages at the unique time indices +regions = zeros(length(time_indices), length(unique_names)) + +for i in eachindex(unique_names) + name = unique_names[i] + region = select(all_voltages, Regex(name)) + temp = reduce(+, eachcol(region)) ./ ncol(region) + println(size(region)) + regions[:, i] = temp[time_indices] +end + +using Plots +plot(time_vec[time_indices], regions, xlabel="Time", ylabel="Voltage", title="Average Voltage in Each Region") \ No newline at end of file From 12e9a381fcfc6e953cf83697793be14e7ee384b0 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 31 Aug 2024 20:44:58 -0400 Subject: [PATCH 23/30] QIF PING connectors --- src/blox/connections.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/blox/connections.jl b/src/blox/connections.jl index cbc0ceb2..c3eb4c86 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -978,3 +978,19 @@ function (bc::BloxConnector)( accumulate_equation!(bc, eq) end +function (bc::BloxConnector)( + bloxout::QIF_PING_NGNMM, + bloxin::QIF_PING_NGNMM; + kwargs... +) + sys_out = get_namespaced_sys(bloxout) + sys_in = get_namespaced_sys(bloxin) + + w = generate_weight_param(bloxout, bloxin; kwargs...) + push!(bc.weights, w) + + x = namespace_expr(bloxout.output, sys_out) + eq = sys_in.jcn ~ w*x + + accumulate_equation!(bc, eq) +end \ No newline at end of file From d2d56a9c92ecf509efe0847e355bcfe941f81e4b Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 31 Aug 2024 20:51:33 -0400 Subject: [PATCH 24/30] Remove noise for now Doesn't look like it's necessary so remove to speed up sims --- examples/adams_example_for_brain_r01.jl | 6 +++++- src/blox/neural_mass.jl | 20 ++++++++++---------- 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/examples/adams_example_for_brain_r01.jl b/examples/adams_example_for_brain_r01.jl index 79cdacac..16799ead 100644 --- a/examples/adams_example_for_brain_r01.jl +++ b/examples/adams_example_for_brain_r01.jl @@ -9,6 +9,11 @@ using Graphs using MetaGraphs using Random +# PING QIF theta-nested gamma oscillations + + + +# Chen/Campbell populations - limited utility for now @named popP = PYR_Izh(η̄=0.08, κ=0.8) @named popQ = PYR_Izh(η̄=0.08, κ=0.2, wⱼ=0.0095, a=0.077) @@ -32,7 +37,6 @@ plot(sol, idxs=[3, 7]) - ## DO NOT TOUCH LARGELY FAILURES THAT RUN sol_dde_with_delays # I want to salvage a bit of this so leaving in for now # -AGC diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index e81e5717..a8ed1a9d 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -549,18 +549,18 @@ struct QIF_PING_NGNMM <: NeuralMassBlox function QIF_PING_NGNMM(; name, namespace=noathing - Δ=0.02, - τₘ, - H, - I_ext, - J_internal, - A) - p = paramscoping(Δ=Δ, τₘ=τₘ, H=H, I_ext=I_ext, J_internal=J_internal, A=A) - Δ, τₘ, H, I_ext, J_internal, A = p + Δ=1.0, + τₘ=20.0, + H=1.3, + I_ext=0.0, + ω=0.0, + J_internal=8.0) + p = paramscoping(Δ=Δ, τₘ=τₘ, H=H, I_ext=I_ext, J_internal=J_internal) + Δ, τₘ, H, I_ext, J_internal = p sts = @variables r(t)=0.0 [output=true] v(t)=0.0 jcn(t)=0.0 [input=true] - @brownian ξ + #@brownian ξ eqs = [D(r) ~ Δ/(π*τₘ^2) + 2*r*v/τₘ, - D(v) ~ (v^2 + H + I_ext)/τₘ - τₘ*(π*r)^2 + J_internal*r + A*ξ + jcn] + D(v) ~ (v^2 + H + I_ext*sin(ω*t))/τₘ - τₘ*(π*r)^2 + J_internal*r + jcn] sys = System(eqs, t, sts, p; name=name) new(p, sts[1], sts[2], sts[3], sys, namespace) end From 82d1e4bdedb2b43a05ae61ec3c5f66072a10e5a9 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 31 Aug 2024 20:51:59 -0400 Subject: [PATCH 25/30] Oops typo --- src/blox/neural_mass.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index a8ed1a9d..c36c0f56 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -548,7 +548,7 @@ struct QIF_PING_NGNMM <: NeuralMassBlox namespace function QIF_PING_NGNMM(; name, - namespace=noathing + namespace=nothing, Δ=1.0, τₘ=20.0, H=1.3, From 1a29ffe59fde797755b4170a0fe40829f5725ffc Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Sat, 31 Aug 2024 22:04:30 -0400 Subject: [PATCH 26/30] Ok it kinda works --- examples/adams_example_for_brain_r01.jl | 23 +++++++++++++++++++++++ src/blox/neural_mass.jl | 7 ++++--- 2 files changed, 27 insertions(+), 3 deletions(-) diff --git a/examples/adams_example_for_brain_r01.jl b/examples/adams_example_for_brain_r01.jl index 16799ead..a7170373 100644 --- a/examples/adams_example_for_brain_r01.jl +++ b/examples/adams_example_for_brain_r01.jl @@ -10,6 +10,29 @@ using MetaGraphs using Random # PING QIF theta-nested gamma oscillations +@named exci_PING = QIF_PING_NGNMM(I_ext=10.0, ω=5*2*π/1000, J_internal=8.0, H=1.3, Δ=1.0, τₘ=20.0, A=0.2) +@named inhi_PING = QIF_PING_NGNMM(I_ext=5.0, ω=5*2*π/1000, J_internal=0.0, H=-5.0, Δ=1.0, τₘ=10.0, A=0.0) + +g = MetaDiGraph() +add_blox!.(Ref(g), [exci_PING, inhi_PING]) +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) + +sim_dur = 1000.0 +prob = SDEProblem(sys, [], (0.0, sim_dur)) +sol = solve(prob, RKMil(), saveat=0.1) + +using Plots, DSP, CSV +plot(sol, idxs=[1, 3]) +plot(sol, idxs=[2, 4]) + +data = Array(sol) +#data = data .+ rand(Normal(0, 10), size(data)) +wp = welch_pgram(data[2, :]; fs=2000) +plot(wp.freq, pow2db.(wp.power), xlim=(0, 50), xlabel="Frequency (Hz)" , ylabel="Power (dB)", title="Power Spectrum of Excitatory Population") diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index c36c0f56..f3fbb795 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -554,13 +554,14 @@ struct QIF_PING_NGNMM <: NeuralMassBlox H=1.3, I_ext=0.0, ω=0.0, - J_internal=8.0) + J_internal=8.0, + A=0.0) p = paramscoping(Δ=Δ, τₘ=τₘ, H=H, I_ext=I_ext, J_internal=J_internal) Δ, τₘ, H, I_ext, J_internal = p sts = @variables r(t)=0.0 [output=true] v(t)=0.0 jcn(t)=0.0 [input=true] - #@brownian ξ + @brownian ξ eqs = [D(r) ~ Δ/(π*τₘ^2) + 2*r*v/τₘ, - D(v) ~ (v^2 + H + I_ext*sin(ω*t))/τₘ - τₘ*(π*r)^2 + J_internal*r + jcn] + D(v) ~ (v^2 + H + I_ext*sin(ω*t))/τₘ - τₘ*(π*r)^2 + J_internal*r + A*ξ + jcn] sys = System(eqs, t, sts, p; name=name) new(p, sts[1], sts[2], sts[3], sys, namespace) end From 920700f46971052a485e2e3f34fdb0ab389df1d1 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 6 Sep 2024 20:53:45 -0400 Subject: [PATCH 27/30] Save voltage version Removing unused utilities and cleaning up the RF learning with saving. --- .../RF_learning_using_BLOX_save_voltages.jl | 306 ++++++++++++++++++ saving_voltage_utils.jl | 46 --- 2 files changed, 306 insertions(+), 46 deletions(-) create mode 100644 examples/RF_learning_using_BLOX_save_voltages.jl delete mode 100644 saving_voltage_utils.jl diff --git a/examples/RF_learning_using_BLOX_save_voltages.jl b/examples/RF_learning_using_BLOX_save_voltages.jl new file mode 100644 index 00000000..af51f215 --- /dev/null +++ b/examples/RF_learning_using_BLOX_save_voltages.jl @@ -0,0 +1,306 @@ +### A Pluto.jl notebook ### +# v0.19.36 + +using Markdown +using InteractiveUtils + +# ╔═╡ 7d500a22-c5ed-48d1-b009-4d4a0f03fc2a +using Pkg + +# ╔═╡ 98f74502-74e6-11ee-34ff-7b16a35b2860 +Pkg.activate(".") + +# ╔═╡ 0e839efb-a15b-410b-93e7-a8668f515423 +Pkg.add("CSV") + +# ╔═╡ 7067bfc7-0e4f-4474-b79e-403e67bfeed8 +Pkg.add("DataFrames") + +# ╔═╡ 038809e1-7512-4f5d-8fbf-feef2cdce117 +Pkg.add("Plots") + +# ╔═╡ b4e8ba6e-e91a-474c-9314-8dd75fa2e82e +Pkg.add("DifferentialEquations") + +# ╔═╡ b6a6209f-97f3-4844-928a-948f3d53b4d7 +using Statistics + +# ╔═╡ 1f184999-9dd8-4958-8158-157b683b6f5c +using CSV + +# ╔═╡ 536b8983-8548-4144-be00-4d516fe7b75f +using DataFrames + +# ╔═╡ e6a453f3-e2a8-44ec-9be2-1f6d9293c911 +using Neuroblox + +# ╔═╡ e8deef43-7ee9-484f-a8a8-2df6e7006e22 +using Plots + +# ╔═╡ ee0c04b9-57c9-44bc-b092-218019a545ec +using DifferentialEquations + +# ╔═╡ d37ba3dd-78df-4e18-ac76-3864a2e9216e +using MetaGraphs + +# ╔═╡ 8abc4a4b-4acf-4c25-a10c-33278329c9f4 +begin + + + time_block_dur = 90 # ms (size of discrete time blocks) + N_trials = 700 #number of trials + + global_ns = :g + + fn = joinpath(@__DIR__, "image_example.csv") #stimulus image file + data = CSV.read(fn, DataFrame) + + #define the circuit blox + @named stim = ImageStimulus(data[1:N_trials,:]; namespace=global_ns, t_stimulus=600, t_pause=1000) + + @named LC = NextGenerationEIBlox(;namespace=global_ns, Cₑ=2*26,Cᵢ=1*26, Δₑ=0.5, Δᵢ=0.5, η_0ₑ=10.0, η_0ᵢ=0.0, v_synₑₑ=10.0, v_synₑᵢ=-10.0, v_synᵢₑ=10.0, v_synᵢᵢ=-10.0, alpha_invₑₑ=10.0/26, alpha_invₑᵢ=0.8/26, alpha_invᵢₑ=10.0/26, alpha_invᵢᵢ=0.8/26, kₑₑ=0.0*26, kₑᵢ=0.6*26, kᵢₑ=0.6*26, kᵢᵢ=0*26) + + @named ITN = NextGenerationEIBlox(;namespace=global_ns, Cₑ=2*36,Cᵢ=1*36, Δₑ=0.5, Δᵢ=0.5, η_0ₑ=10.0, η_0ᵢ=0.0, v_synₑₑ=10.0, v_synₑᵢ=-10.0, v_synᵢₑ=10.0, v_synᵢᵢ=-10.0, alpha_invₑₑ=10.0/36, alpha_invₑᵢ=0.8/36, alpha_invᵢₑ=10.0/36, alpha_invᵢᵢ=0.8/36, kₑₑ=0.0*36, kₑᵢ=0.6*36, kᵢₑ=0.6*36, kᵢᵢ=0*36) + + @named VC = CorticalBlox(N_wta=45, N_exci=5, density=0.01, weight=1,I_bg_ar=0;namespace=global_ns) + + @named PFC = CorticalBlox(N_wta=20, N_exci=5, density=0.01, weight=1,I_bg_ar=0;namespace=global_ns) + + @named STR1 = Striatum(N_inhib=25;namespace=global_ns) + @named STR2 = Striatum(N_inhib=25;namespace=global_ns) + + @named tan_nrn = HHNeuronExciBlox(;namespace=global_ns) + + @named gpi1 = GPi(N_inhib=25;namespace=global_ns) + @named gpi2 = GPi(N_inhib=25;namespace=global_ns) + + @named gpe1 = GPe(N_inhib=15;namespace=global_ns) + @named gpe2 = GPe(N_inhib=15;namespace=global_ns) + + @named STN1 = STN(N_exci=15,I_bg=3*ones(25);namespace=global_ns) + @named STN2 = STN(N_exci=15,I_bg=3*ones(25);namespace=global_ns) + + @named Thal1 = Thalamus(N_exci=25;namespace=global_ns) + @named Thal2 = Thalamus(N_exci=25;namespace=global_ns) + + + @named tan_pop1 = TAN(;namespace=global_ns) + @named tan_pop2 = TAN(;namespace=global_ns) + + @named AS = GreedyPolicy(namespace=global_ns, t_decision=180.5) + @named SNcb = SNc(namespace=global_ns) + + + assembly = [LC, ITN, VC, PFC, STR1, STR2, tan_nrn, gpi1, gpi2, gpe1, gpe2, STN1, STN2, Thal1, Thal2, stim, tan_pop1, tan_pop2, AS, SNcb] + d = Dict(b => i for (i,b) in enumerate(assembly)) + + + #define learning rules + hebbian_mod = HebbianModulationPlasticity(K=0.04, decay=0.01, α=2.5, θₘ=1, modulator=SNcb, t_pre=1600-eps(), t_post=1600-eps(), t_mod=90) + + hebbian_cort = HebbianPlasticity(K=5e-4, W_lim=5, t_pre=1600-eps(), t_post=1600-eps()) + + hebbian_thal_cort = HebbianPlasticity(K=1.7e-5, W_lim=6, t_pre=1600-eps(), t_post=1600-eps()) + + + g = MetaDiGraph() + add_blox!.(Ref(g), assembly) + + #connect the bloxs + + add_edge!(g, d[LC], d[VC], Dict(:weight => 44)) #LC->VC + add_edge!(g, d[LC], d[PFC], Dict(:weight => 44)) #LC->pfc + add_edge!(g, d[ITN], d[tan_nrn], Dict(:weight => 100)) #ITN->tan + add_edge!(g, d[VC], d[PFC], Dict(:weight => 1, :density => 0.08, :learning_rule => hebbian_cort)) #VC->pfc + add_edge!(g, d[PFC], d[STR1], Dict(:weight => 0.075, :density => 0.04, :learning_rule => hebbian_mod)) #pfc->str1 + add_edge!(g, d[PFC], d[STR2], Dict(:weight => 0.075, :density => 0.04, :learning_rule => hebbian_mod)) #pfc->str2 + add_edge!(g, d[tan_nrn], d[STR1], Dict(:weight => 0.17)) #tan->str1 + add_edge!(g, d[tan_nrn], d[STR2], Dict(:weight => 0.17)) #tan->str2 + add_edge!(g, d[STR1], d[gpi1], Dict(:weight => 4, :density => 0.04)) #str1->gpi1 + add_edge!(g, d[STR2], d[gpi2], Dict(:weight => 4, :density => 0.04)) #str2->gpi2 + add_edge!(g, d[gpi1], d[Thal1], Dict(:weight => 0.16, :density => 0.04)) #gpi1->thal1 + add_edge!(g, d[gpi2], d[Thal2], Dict(:weight => 0.16, :density => 0.04)) #gpi2->thal2 + add_edge!(g, d[Thal1], d[PFC], Dict(:weight => 0.2, :density => 0.32, :learning_rule => hebbian_thal_cort, :sta => true)) #thal1->pfc + add_edge!(g, d[Thal2], d[PFC], Dict(:weight => 0.2, :density => 0.32, :learning_rule => hebbian_thal_cort, :sta => true)) #thal2->pfc + add_edge!(g, d[STR1], d[gpe1], Dict(:weight => 4, :density => 0.04)) #str1->gpe1 + add_edge!(g, d[STR2], d[gpe2], Dict(:weight => 4, :density => 0.04)) #str2->gpe2 + add_edge!(g, d[gpe1], d[gpi1], Dict(:weight => 0.2, :density => 0.04)) #gpe1->gpi1 + add_edge!(g, d[gpe2], d[gpi2], Dict(:weight => 0.2, :density => 0.04)) #gpe2->gpi2 + add_edge!(g, d[gpe1], d[STN1], Dict(:weight => 3.5, :density => 0.04)) #gpe1->stn1 + add_edge!(g, d[gpe2], d[STN2], Dict(:weight => 3.5, :density => 0.04)) #gpe2->stn2 + add_edge!(g, d[STN1], d[gpi1], Dict(:weight => 0.1, :density => 0.04)) #stn1->gpi1 + add_edge!(g, d[STN2], d[gpi2], Dict(:weight => 0.1, :density => 0.04)) #stn2->gpi2 + add_edge!(g, d[stim], d[VC], :weight, 14) #stim->VC + add_edge!(g, d[tan_pop1], d[STR1], Dict(:weight => 1, :t_event => 90.0)) #TAN pop1 -> str1 + add_edge!(g, d[tan_pop2], d[STR2], Dict(:weight => 1, :t_event => 90.0)) #TAN pop2 -> str2 + add_edge!(g, d[STR1], d[tan_pop1], Dict(:weight => 1)) #str1 -> TAN pop1 + add_edge!(g, d[STR2], d[tan_pop1], Dict(:weight => 1)) #str2 -> TAN pop1 + add_edge!(g, d[STR1], d[tan_pop2], Dict(:weight => 1)) #str1 -> TAN pop2 + add_edge!(g, d[STR2], d[tan_pop2], Dict(:weight => 1)) #str2 -> TAN pop2 + add_edge!(g, d[STR1], d[STR2], Dict(:weight => 1, :t_event => 181.0)) #str1 -> str2 + add_edge!(g, d[STR2], d[STR1], Dict(:weight => 1, :t_event => 181.0)) #str2 -> str1 + add_edge!(g, d[STR1], d[AS])# str1->AS + add_edge!(g, d[STR2], d[AS])# str2->AS + add_edge!(g, d[STR1], d[SNcb], Dict(:weight => 1)) # str1->Snc + add_edge!(g, d[STR2], d[SNcb], Dict(:weight => 1)) # str2->Snc + + +end + +# ╔═╡ 5ddd418f-23d2-4151-82c0-e8c9fae7e4f2 +#define the circuit as an Agent + agent = Agent(g; name=:ag, t_block = 90); + +# ╔═╡ 7576c26d-8617-4c66-ab5e-db0bd8bb9e17 + prob=agent.problem + +# ╔═╡ d80b3bf7-24d8-4395-8903-de974e6445f9 +begin + #extract membrane voltages of every neuron + getsys=agent.odesystem; + st=unknowns(getsys) + vlist=Int64[] + for ii = 1:length(st) + if contains(string(st[ii]), "V(t)") + push!(vlist,ii) + end + end +end + +# ╔═╡ bf1f09d9-7f9e-48f7-b3de-9fecae82a5a0 +begin +#define environment : contains stimuli and feedback + env = ClassificationEnvironment(stim; name=:env, namespace=global_ns) + +end + +# ╔═╡ 14e2806e-5cdd-4b0b-b8c7-8c6d2c7e349c +begin + + t_warmup=800 + t_trial = env.t_trial + tspan = (0, t_trial) + stim_params_in = Neuroblox.get_trial_stimulus(env) + sys = Neuroblox.get_sys(agent) + prob2 = remake(prob; p = merge(stim_params_in),tspan=(0,1600)) + + if t_warmup > 0 + prob2 = remake(prob2; tspan=(0,t_warmup)) + sol = solve(prob2, Vern7()) + u0 = sol[1:end,end] # last value of state vector + prob2 = remake(prob2; p = merge(stim_params_in),tspan=tspan, u0=u0) + else + prob2 = remake(prob3; p = merge(stim_params_in),tspan=tspan) + u0 = [] + end + + action_selection = agent.action_selection + learning_rules = agent.learning_rules + + defs = ModelingToolkit.get_defaults(sys) + weights = Dict{Num, Float64}() + for w in keys(learning_rules) + weights[w] = defs[w] + end + + perf=zeros(N_trials) + act=zeros(N_trials) +end + + +# ╔═╡ b7e84b20-0b80-478c-bc88-2883f80bcbb4 +begin + # sys = Neuroblox.get_sys(agent) + # learning_rules = agent.learning_rules + # weights = Dict{Num, Float64}() + # for w in keys(learning_rules) + # weights[w] = defs[w] + # end + + for ii = 1:N_trials + prob3 = agent.problem + stim_params = Neuroblox.get_trial_stimulus(env) + new_params = ModelingToolkit.MTKParameters(sys, merge(defs, weights, stim_params)) + prob3 = remake(prob3; p = new_params, u0=u0, tspan=(0,env.t_trial)) + @info env.current_trial + sol2 = solve(prob3, Vern7()) + agent.problem = prob3 + + if isnothing(action_selection) + feedback = 1 + else + action = action_selection(sol2) + feedback = env(action) + @info action + @info env.category[env.current_trial] + @info feedback + perf[ii]=feedback + act[ii] = action + end + + for (w, rule) in learning_rules + w_val = weights[w] + Δw = Neuroblox.weight_gradient(rule, sol2, w_val, feedback) + weights[w] += Δw + end + Neuroblox.increment_trial!(env) + + df = DataFrame() + df.t = sol2.t + df.VC = Neuroblox.meanfield_timeseries(VC, sol2) + df.PFC = Neuroblox.meanfield_timeseries(PFC, sol2) + df.STR1 = Neuroblox.meanfield_timeseries(STR1, sol2) + df.STR2 = Neuroblox.meanfield_timeseries(STR2, sol2) + df.gpi1 = Neuroblox.meanfield_timeseries(gpi1, sol2) + df.gpi2 = Neuroblox.meanfield_timeseries(gpi2, sol2) + df.gpe1 = Neuroblox.meanfield_timeseries(gpe1, sol2) + df.gpe2 = Neuroblox.meanfield_timeseries(gpe2, sol2) + df.STN1 = Neuroblox.meanfield_timeseries(STN1, sol2) + df.STN2 = Neuroblox.meanfield_timeseries(STN2, sol2) + df.Thal1 = Neuroblox.meanfield_timeseries(Thal1, sol2) + df.Thal2 = Neuroblox.meanfield_timeseries(Thal2, sol2) + CSV.write("sol_trial_$ii.csv", df, writeheader=true) + end + +end + +# ╔═╡ 9703f279-a311-4aa2-85e4-75d8898b21b6 +Gray.(act/2) + +# ╔═╡ a973a05c-e686-48c5-8242-688e6475624b +Gray.(perf) + +# ╔═╡ 35891714-fafc-4465-a2e5-bcd7aaf1a59c +begin + perf_smth=zeros(length(perf)-19) + for ii=20:length(perf) + perf_smth[ii-19] = mean(perf[ii-19:ii]) + + end + plot(collect(20:500),perf_smth[1:481],ylims=(0,1),xlims=(0,500)) +end + +# ╔═╡ Cell order: +# ╠═7d500a22-c5ed-48d1-b009-4d4a0f03fc2a +# ╠═98f74502-74e6-11ee-34ff-7b16a35b2860 +# ╠═b6a6209f-97f3-4844-928a-948f3d53b4d7 +# ╠═0e839efb-a15b-410b-93e7-a8668f515423 +# ╠═1f184999-9dd8-4958-8158-157b683b6f5c +# ╠═7067bfc7-0e4f-4474-b79e-403e67bfeed8 +# ╠═536b8983-8548-4144-be00-4d516fe7b75f +# ╠═e6a453f3-e2a8-44ec-9be2-1f6d9293c911 +# ╠═038809e1-7512-4f5d-8fbf-feef2cdce117 +# ╠═e8deef43-7ee9-484f-a8a8-2df6e7006e22 +# ╠═b4e8ba6e-e91a-474c-9314-8dd75fa2e82e +# ╠═ee0c04b9-57c9-44bc-b092-218019a545ec +# ╠═d37ba3dd-78df-4e18-ac76-3864a2e9216e +# ╠═8abc4a4b-4acf-4c25-a10c-33278329c9f4 +# ╠═5ddd418f-23d2-4151-82c0-e8c9fae7e4f2 +# ╠═7576c26d-8617-4c66-ab5e-db0bd8bb9e17 +# ╠═d80b3bf7-24d8-4395-8903-de974e6445f9 +# ╠═bf1f09d9-7f9e-48f7-b3de-9fecae82a5a0 +# ╠═14e2806e-5cdd-4b0b-b8c7-8c6d2c7e349c +# ╠═b7e84b20-0b80-478c-bc88-2883f80bcbb4 +# ╠═9703f279-a311-4aa2-85e4-75d8898b21b6 +# ╠═a973a05c-e686-48c5-8242-688e6475624b +# ╠═35891714-fafc-4465-a2e5-bcd7aaf1a59c diff --git a/saving_voltage_utils.jl b/saving_voltage_utils.jl deleted file mode 100644 index af7ce44e..00000000 --- a/saving_voltage_utils.jl +++ /dev/null @@ -1,46 +0,0 @@ -using FileIO -using DataFrames - -# Loading single file - will loop to load all eventually -# save as csv -data = load("/Users/achesebro/Downloads/sim0003.jld2") -df = data["df"] -all_voltages = select(df, r"₊V") -time_vec = select(df, :timestamp)[:, 1] - -# Get all region names - -# Make into utility function -# Create test (eventually) -all_names = names(all_voltages) -unique_names = Vector{String}() - -for i in eachindex(all_names) - name_temp = all_names[i] - name = split(name_temp, "₊")[1] #check for innernamespaceof for global namespace removal - push!(unique_names, name) -end - -unique_names = unique(unique_names) - -# Get unique time indices (not originally unique because of the callbacks) -time_indices = Vector{Int}() -for i in eachindex(time_vec) - if i == 1 || time_vec[i] != time_vec[i-1] - push!(time_indices, i) - end -end - -# Extract the voltages at the unique time indices -regions = zeros(length(time_indices), length(unique_names)) - -for i in eachindex(unique_names) - name = unique_names[i] - region = select(all_voltages, Regex(name)) - temp = reduce(+, eachcol(region)) ./ ncol(region) - println(size(region)) - regions[:, i] = temp[time_indices] -end - -using Plots -plot(time_vec[time_indices], regions, xlabel="Time", ylabel="Voltage", title="Average Voltage in Each Region") \ No newline at end of file From 75d3ea482ad2988369642d37c0d01794f8f81785 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 13 Sep 2024 11:58:57 -0400 Subject: [PATCH 28/30] Better processing of saved voltages --- examples/RF_processing_saved_voltages.jl | 26 ++++++++++++++++++++++++ 1 file changed, 26 insertions(+) create mode 100644 examples/RF_processing_saved_voltages.jl diff --git a/examples/RF_processing_saved_voltages.jl b/examples/RF_processing_saved_voltages.jl new file mode 100644 index 00000000..3e0f9913 --- /dev/null +++ b/examples/RF_processing_saved_voltages.jl @@ -0,0 +1,26 @@ +using CSV, DataFrames +using Interpolations + +function collapse_timesteps(data) + new_data = DataFrame() + for i ∈ eachindex(data.t) + if i > 1 && data.t[i] != data.t[i-1] + push!(new_data, data[i, :]) + end + end + return new_data +end + +function resample_timeseries(data, dt=1) + new_data = DataFrame() + new_t = 0:dt:data.t[end] + new_data.t = new_t + for i ∈ eachindex(Array(data[1, 2:end])) + temp_interp = linear_interpolation(data.t, data[:, i+1], extrapolation_bc=Line()) + new_data[!, names(data)[i+1]] = temp_interp(new_t) + end + return new_data +end + + +data = DataFrame(CSV.File("/Users/achesebro/Downloads/all_sims/sol_trial_1.csv")) From f802bfd8c86956576b077753aea575d2a98a1a20 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 13 Sep 2024 13:45:34 -0400 Subject: [PATCH 29/30] Update tests for new blocks --- examples/adams_example_for_brain_r01.jl | 4 +-- test/components.jl | 36 +++++++++++++++++++++++++ 2 files changed, 38 insertions(+), 2 deletions(-) diff --git a/examples/adams_example_for_brain_r01.jl b/examples/adams_example_for_brain_r01.jl index cada5590..29158a7b 100644 --- a/examples/adams_example_for_brain_r01.jl +++ b/examples/adams_example_for_brain_r01.jl @@ -71,8 +71,8 @@ f g = MetaDiGraph() add_blox!.(Ref(g), [popP, popQ]) -add_edge!(g, popP => popQ; weight=1.0) #weight is acutaally meaningless here -add_edge!(g, popQ => popP; weight=1.0) #weight is acutaally meaningless here +add_edge!(g, popP => popQ; weight=1.0) #weight is acutally meaningless here +add_edge!(g, popQ => popP; weight=1.0) #weight is acutally meaningless here @named sys = system_from_graph(g) sys = structural_simplify(sys) diff --git a/test/components.jl b/test/components.jl index e278a66c..d16f0967 100644 --- a/test/components.jl +++ b/test/components.jl @@ -772,3 +772,39 @@ end sol = solve(prob, Tsit5()) @test sol.retcode == ReturnCode.Success end + +@testset "QIF_PING" begin + @named exci_PING = QIF_PING_NGNMM(I_ext=10.0, ω=5*2*π/1000, J_internal=8.0, H=1.3, Δ=1.0, τₘ=20.0, A=0.2) + @named inhi_PING = QIF_PING_NGNMM(I_ext=5.0, ω=5*2*π/1000, J_internal=0.0, H=-5.0, Δ=1.0, τₘ=10.0, A=0.0) + + g = MetaDiGraph() + add_blox!.(Ref(g), [exci_PING, inhi_PING]) + 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) + + sim_dur = 100.0 + prob = SDEProblem(sys, [], (0.0, sim_dur)) + sol = solve(prob, RKMil(), saveat=0.1) + @test sol.retcode == ReturnCode.Success +end + +@testset "PYR_Izh" begin + @named popP = PYR_Izh(η̄=0.08, κ=0.8) + @named popQ = PYR_Izh(η̄=0.08, κ=0.2, wⱼ=0.0095, a=0.077) + + g = MetaDiGraph() + add_blox!.(Ref(g), [popP, popQ]) + add_edge!(g, popP => popQ; weight=1.0) + add_edge!(g, popQ => popP; weight=1.0) + + @named sys = system_from_graph(g) + sys = structural_simplify(sys) + + sim_dur = 200.0 + prob = ODEProblem(sys, [], (0.0, sim_dur)) + sol = solve(prob, Tsit5(), saveat=1.0) + @test sol.retcode == ReturnCode.Success +end \ No newline at end of file From d9348ff124c6faaa801290de69236c608233e1e6 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Fri, 13 Sep 2024 13:46:39 -0400 Subject: [PATCH 30/30] Organization --- .../RF_learning_using_BLOX_save_voltages.jl | 0 examples/{ => flattening_circuit}/RF_processing_saved_voltages.jl | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename examples/{ => flattening_circuit}/RF_learning_using_BLOX_save_voltages.jl (100%) rename examples/{ => flattening_circuit}/RF_processing_saved_voltages.jl (100%) diff --git a/examples/RF_learning_using_BLOX_save_voltages.jl b/examples/flattening_circuit/RF_learning_using_BLOX_save_voltages.jl similarity index 100% rename from examples/RF_learning_using_BLOX_save_voltages.jl rename to examples/flattening_circuit/RF_learning_using_BLOX_save_voltages.jl diff --git a/examples/RF_processing_saved_voltages.jl b/examples/flattening_circuit/RF_processing_saved_voltages.jl similarity index 100% rename from examples/RF_processing_saved_voltages.jl rename to examples/flattening_circuit/RF_processing_saved_voltages.jl