Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

initial version of the basal ganglia and Parkinson tutorial #22

Merged
merged 1 commit into from
Sep 27, 2024
Merged

initial version of the basal ganglia and Parkinson tutorial #22

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
198 changes: 198 additions & 0 deletions docs/src/tutorials/basal_ganglia_parkinson.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,198 @@
using Neuroblox
using DifferentialEquations
using Graphs
using MetaGraphs
using CairoMakie
# using Statistics
# using Plots
# using DSP
# using Peaks
# using SparseArrays
# using Random


# utility function to include in Neuroblox
get_system(blox) = blox.odesystem
# get_system(blox::AbstractBlox) = blox.odesystem


# Isolated MSN network in baseline condition
@named msn = Striatum_MSN_Adam();
sys = get_system(msn);
sys = structural_simplify(sys)

# check all the 100 neurons, each with its associated associated currents
unknowns(sys)
prob = SDEProblem(sys, [], (0.0, 5500.0), [])
sol = solve(prob, RKMil(); dt=0.05, saveat=0.05)

# plot voltage of a single neuron
plot(sol, idxs=1, axis = (xlabel = "time (ms)", ylabel = "membrane potential (mV)"))

# plot mean field
meanfield(msn, sol, axis = (xlabel = "time (ms)", ylabel = "membrane potential (mV)", title = "Mean Field"))

# get mean firing rate
spikes = detect_spikes(msn, sol; threshold=-55)
t, fr = mean_firing_rate(spikes, sol)

# raster plot
rasterplot(fig[1,3], msn, sol, threshold = -55.0, axis = (; title = "Neuron's Spikes - Mean Firing Rate: $(fr[1]) spikes/s"))

# power spectrum of the GABAa current


fig = Figure(resolution = (1500, 600))
powerspectrumplot(fig[1,1], msn, sol; state = "G",
ylims=(1e-5, 10),
alpha_start = 5,
alpha_label_position = (8.5, 4),
beta_label_position = (22, 4),
gamma_label_position = (60, 4),
axis = (; title = "Periodogram with no window"))


powerspectrumplot(fig[1,2], msn, sol; state = "G",
method=welch_pgram, window=hanning,
ylims=(1e-5, 10),
alpha_start = 5,
alpha_label_position = (8.5, 4),
beta_label_position = (22, 4),
gamma_label_position = (60, 4),
axis = (; title = "Welch's method with no Hanning window"))


# ens_prob = EnsembleProblem(prob)
# ens_sol = solve(ens_prob, RKMil(); dt=0.05, saveat=0.05, trajectories=5)

global_ns = :g
# Basal ganglia model in baseline condition
@named msn = Striatum_MSN_Adam(namespace=global_ns)
@named fsi = Striatum_FSI_Adam(namespace=global_ns)
@named gpe = GPe_Adam(namespace=global_ns)
@named stn = STN_Adam(namespace=global_ns)


assembly = [msn, fsi, gpe, stn]
g = MetaDiGraph()
add_blox!.(Ref(g), assembly)
add_edge!(g, 1, 3, Dict(:weight=> 2.5/33, :density=>0.33))
add_edge!(g, 2, 1, Dict(:weight=> 0.6/7.5, :density=>0.15))
add_edge!(g, 3, 4, Dict(:weight=> 0.3/4, :density=>0.05))
add_edge!(g, 4, 2, Dict(:weight=> 0.165/4, :density=>0.1))
# the fractions above represent ḡ_inh/number of presynaptic neurons

@named neuron_net = system_from_graph(g)
sys = structural_simplify(neuron_net)
prob = SDEProblem(sys, [], (0.0, 5500.0), [])
sol = solve(prob, RKMil(); dt=0.05, saveat=0.05)



fig = Figure(resolution = (1500, 600))
powerspectrumplot(fig[1,1], msn, sol; state = "G",
method=welch_pgram, window=hanning,
ylims=(1e-5, 10),
alpha_start = 5,
alpha_label_position = (8.5, 3),
beta_label_position = (22, 3),
gamma_label_position = (60, 3),
axis = (; title = "MSN"))


powerspectrumplot(fig[1,2], fsi, sol; state = "G",
method=welch_pgram, window=hanning,
ylims=(1e-5, 10),
alpha_start = 5,
alpha_label_position = (8.5, 3),
beta_label_position = (22, 3),
gamma_label_position = (60, 3),
axis = (; title = "FSI"))

powerspectrumplot(fig[1,3], gpe, sol; state = "G",
method=welch_pgram, window=hanning,
ylims=(1e-5, 10),
alpha_start = 5,
alpha_label_position = (8.5, 3),
beta_label_position = (22, 3),
gamma_label_position = (60, 3),
axis = (; title = "GPe"))


powerspectrumplot(fig[1,4], stn, sol; state = "G",
method=welch_pgram, window=hanning,
ylims=(1e-5, 10),
alpha_start = 5,
alpha_label_position = (8.5, 3),
beta_label_position = (22, 3),
gamma_label_position = (60, 3),
axis = (; title = "STN"))


fig




# Parkinson's condition

@named msn = Striatum_MSN_Adam(namespace=global_ns, I_bg = 1.2519*ones(100), G_M = 1.2)
@named fsi = Striatum_FSI_Adam(namespace=global_ns, I_bg = 4.511*ones(50), weight = 0.2, g_weight = 0.075)

assembly = [msn, fsi, gpe, stn]
g = MetaDiGraph()
add_blox!.(Ref(g), assembly)
ḡ_inh = 0.48
add_edge!(g, 2, 1, Dict(:weight=> ḡ_inh/7.5, :density=>0.15))
add_edge!(g, 1, 3, Dict(:weight=> 2.5/33, :density=>0.33))
add_edge!(g, 3, 4, Dict(:weight=> 0.3/4, :density=>0.05))
add_edge!(g, 4, 2, Dict(:weight=> 0.165/4, :density=>0.1))
# the fractions above represent ḡ_inh/number of presynaptic neurons

@named neuron_net = system_from_graph(g)
sys = structural_simplify(neuron_net)
prob = SDEProblem(sys, [], (0.0, 5500.0), [])
sol = solve(prob, RKMil(); dt=0.05, saveat=0.05)



powerspectrumplot(fig[2,1], msn, sol; state = "G",
method=welch_pgram, window=hanning,
ylims=(1e-5, 10),
alpha_start = 5,
alpha_label_position = (8.5, 3),
beta_label_position = (22, 3),
gamma_label_position = (60, 3),
axis = (; title = "MSN"))


powerspectrumplot(fig[2,2], fsi, sol; state = "G",
method=welch_pgram, window=hanning,
ylims=(1e-5, 10),
alpha_start = 5,
alpha_label_position = (8.5, 3),
beta_label_position = (22, 3),
gamma_label_position = (60, 3),
axis = (; title = "FSI"))

powerspectrumplot(fig[2,3], gpe, sol; state = "G",
method=welch_pgram, window=hanning,
ylims=(1e-5, 10),
alpha_start = 5,
alpha_label_position = (8.5, 3),
beta_label_position = (22, 3),
gamma_label_position = (60, 3),
axis = (; title = "GPe"))


powerspectrumplot(fig[2,4], stn, sol; state = "G",
method=welch_pgram, window=hanning,
ylims=(1e-5, 10),
alpha_start = 5,
alpha_label_position = (8.5, 3),
beta_label_position = (22, 3),
gamma_label_position = (60, 3),
axis = (; title = "STN"))


fig
Loading