From 1ab07ca11b5b0f96e4cbf66d106a6f5ce684fb8d Mon Sep 17 00:00:00 2001
From: agchesebro <76024790+agchesebro@users.noreply.github.com>
Date: Mon, 29 Jul 2024 18:17:08 -0400
Subject: [PATCH 01/10] Initial commit with to-do list

---
 src/Neuroblox.jl  | 1 +
 src/sermon_dbs.jl | 7 +++++++
 2 files changed, 8 insertions(+)
 create mode 100644 src/sermon_dbs.jl

diff --git a/src/Neuroblox.jl b/src/Neuroblox.jl
index a1e3138a..bf77f7c4 100644
--- a/src/Neuroblox.jl
+++ b/src/Neuroblox.jl
@@ -112,6 +112,7 @@ include("gui/GUI.jl")
 include("blox/connections.jl")
 include("blox/blox_utilities.jl")
 include("Neurographs.jl")
+include("blox/sermon_dbs.jl")
 
 function simulate(sys::ODESystem, u0, timespan, p, solver = AutoVern7(Rodas4()); kwargs...)
     prob = ODEProblem(sys, u0, timespan, p)
diff --git a/src/sermon_dbs.jl b/src/sermon_dbs.jl
new file mode 100644
index 00000000..9eba09b8
--- /dev/null
+++ b/src/sermon_dbs.jl
@@ -0,0 +1,7 @@
+# TODO: List of missing components needed for this tutorial
+# 1. SuperBlox of Kuramoto oscillators to collect dynamics
+# 2. DBS stimulator block as implemented in the paper
+# 3. Neurotransmitter pool blocks as discussed in the paper
+# 4. BloxConnector functionality for neurotransmitter pools
+# 5. BloxConnector functionality for DBS stimulator
+

From 45877e7c6a28b83b7ef7a57a7627aac719d74cad Mon Sep 17 00:00:00 2001
From: agchesebro <76024790+agchesebro@users.noreply.github.com>
Date: Tue, 30 Jul 2024 09:01:08 -0400
Subject: [PATCH 02/10] Create neurotransmitter pool block

---
 src/blox/sermon_dbs.jl | 28 ++++++++++++++++++++++++++++
 src/sermon_dbs.jl      |  7 -------
 2 files changed, 28 insertions(+), 7 deletions(-)
 create mode 100644 src/blox/sermon_dbs.jl
 delete mode 100644 src/sermon_dbs.jl

diff --git a/src/blox/sermon_dbs.jl b/src/blox/sermon_dbs.jl
new file mode 100644
index 00000000..81621cba
--- /dev/null
+++ b/src/blox/sermon_dbs.jl
@@ -0,0 +1,28 @@
+# TODO: List of missing components needed for this tutorial
+# 1. SuperBlox of Kuramoto oscillators to collect dynamics
+# 2. DBS stimulator block as implemented in the paper
+# 4. BloxConnector functionality for neurotransmitter pools
+# 5. BloxConnector functionality for DBS stimulator
+
+# Create neurotransmitter pools 
+struct SermonNPool <: NeuralMassBlox
+    params
+    output
+    jcn
+    odesystem
+    namespace
+    function SermonNPool(;
+                        name,
+                        namespace=nothing,
+                        τ_pool=0.1,
+                        p_pool=0.3
+            )
+        p = paramscoping(τ_pool=τ_pool, p_pool=p_pool)
+        τ_pool, p_pool = p
+        
+        sts = @variables n_pool(t)=0.0 [output = true] jcn(t)=0.0 [input=true]
+        eqs = [D(n_pool) ~ (1-n_pool)/τ_pool - p_pool*n_pool*jcn]
+        sys = System(eqs, t, sts, p; name=name)
+        new(p, sts[1], sts[2], sys, namespace)
+    end
+end
diff --git a/src/sermon_dbs.jl b/src/sermon_dbs.jl
deleted file mode 100644
index 9eba09b8..00000000
--- a/src/sermon_dbs.jl
+++ /dev/null
@@ -1,7 +0,0 @@
-# TODO: List of missing components needed for this tutorial
-# 1. SuperBlox of Kuramoto oscillators to collect dynamics
-# 2. DBS stimulator block as implemented in the paper
-# 3. Neurotransmitter pool blocks as discussed in the paper
-# 4. BloxConnector functionality for neurotransmitter pools
-# 5. BloxConnector functionality for DBS stimulator
-

From 110180d1508a0da5541734c69b7f542a91dcef32 Mon Sep 17 00:00:00 2001
From: agchesebro <76024790+agchesebro@users.noreply.github.com>
Date: Tue, 30 Jul 2024 14:15:55 -0400
Subject: [PATCH 03/10] Add DBS block

---
 src/Neuroblox.jl       |  1 +
 src/blox/sermon_dbs.jl | 30 ++++++++++++++++++++++++++++++
 2 files changed, 31 insertions(+)

diff --git a/src/Neuroblox.jl b/src/Neuroblox.jl
index bf77f7c4..bc7c212d 100644
--- a/src/Neuroblox.jl
+++ b/src/Neuroblox.jl
@@ -212,5 +212,6 @@ export run_experiment!, run_trial!
 export addnontunableparams
 export get_weights, get_dynamic_states, get_idx_tagged_vars, get_eqidx_tagged_vars
 export BalloonModel, boldsignal_endo_balloon
+export SermonNPool, SermonDBS
 
 end
diff --git a/src/blox/sermon_dbs.jl b/src/blox/sermon_dbs.jl
index 81621cba..c7307ba5 100644
--- a/src/blox/sermon_dbs.jl
+++ b/src/blox/sermon_dbs.jl
@@ -1,3 +1,6 @@
+using ModelingToolkitStandardLibrary.Blocks: smooth_square, smooth_step
+using IfElse
+
 # TODO: List of missing components needed for this tutorial
 # 1. SuperBlox of Kuramoto oscillators to collect dynamics
 # 2. DBS stimulator block as implemented in the paper
@@ -26,3 +29,30 @@ struct SermonNPool <: NeuralMassBlox
         new(p, sts[1], sts[2], sys, namespace)
     end
 end
+
+
+# Create a new DBS stimulator block
+# Largely based on MTK Standard Library smooth_square function in src/Blocks/sources.jl
+struct SermonDBS <: NeuralMassBlox
+    params
+    output
+    jcn
+    odesystem
+    namespace
+    function SermonDBS(;
+                        name,
+                        namespace=nothing,
+                        f_stim=130,
+                        amplitude=1.0,
+                        pulse_width=0.0001,
+                        start_time=0.0      
+            )
+        p = paramscoping(f_stim=f_stim, amplitude=amplitude, pulse_width=pulse_width, start_time=start_time)
+        f_stim, amplitude, pulse_width, start_time = p
+        
+        sts = @variables u(t)=0.0 [output = true]
+        eqs = [u ~ IfElse.ifelse(t > start_time, IfElse.ifelse(t % (1.0/f_stim) < pulse_width, amplitude, 0), 0)]
+        sys = System(eqs, t, sts, p; name=name)
+        new(p, sts[1], nothing, sys, namespace)
+    end
+end

From 61a9ab94c48f11d6870a2385823bad28b25f3db5 Mon Sep 17 00:00:00 2001
From: agchesebro <76024790+agchesebro@users.noreply.github.com>
Date: Tue, 30 Jul 2024 14:38:16 -0400
Subject: [PATCH 04/10] Create connector for DBS stimulator to neurotransmitter
 pool

---
 src/blox/sermon_dbs.jl | 20 +++++++++++++++++++-
 1 file changed, 19 insertions(+), 1 deletion(-)

diff --git a/src/blox/sermon_dbs.jl b/src/blox/sermon_dbs.jl
index c7307ba5..4894e97e 100644
--- a/src/blox/sermon_dbs.jl
+++ b/src/blox/sermon_dbs.jl
@@ -23,7 +23,7 @@ struct SermonNPool <: NeuralMassBlox
         p = paramscoping(τ_pool=τ_pool, p_pool=p_pool)
         τ_pool, p_pool = p
         
-        sts = @variables n_pool(t)=0.0 [output = true] jcn(t)=0.0 [input=true]
+        sts = @variables n_pool(t)=1.0 [output = true] jcn(t)=0.0 [input=true]
         eqs = [D(n_pool) ~ (1-n_pool)/τ_pool - p_pool*n_pool*jcn]
         sys = System(eqs, t, sts, p; name=name)
         new(p, sts[1], sts[2], sys, namespace)
@@ -56,3 +56,21 @@ struct SermonDBS <: NeuralMassBlox
         new(p, sts[1], nothing, sys, namespace)
     end
 end
+
+function (bc::BloxConnector)(
+    bloxout::SermonDBS, 
+    bloxin::SermonNPool; 
+    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
+

From 47a550c094fa78e37461a7fa69467c025dc1d1d7 Mon Sep 17 00:00:00 2001
From: agchesebro <76024790+agchesebro@users.noreply.github.com>
Date: Tue, 30 Jul 2024 16:19:34 -0400
Subject: [PATCH 05/10] Custom connector for Sermon Kuramoto connections

---
 src/blox/sermon_dbs.jl | 29 +++++++++++++++++++++++++++++
 1 file changed, 29 insertions(+)

diff --git a/src/blox/sermon_dbs.jl b/src/blox/sermon_dbs.jl
index 4894e97e..61eee082 100644
--- a/src/blox/sermon_dbs.jl
+++ b/src/blox/sermon_dbs.jl
@@ -74,3 +74,32 @@ function (bc::BloxConnector)(
     accumulate_equation!(bc, eq)
 end
 
+# Redfine Kuramoto oscillator connections with arbitrary connection rule
+function (bc::BloxConnector)(
+    bloxout::KuramotoOscillator, 
+    bloxin::KuramotoOscillator; 
+    kwargs...
+)
+    sys_out = get_namespaced_sys(bloxout)
+    sys_in = get_namespaced_sys(bloxin)
+
+    #technically these two lines aren't needed, but useful to have if weighting occurs here
+    w = generate_weight_param(bloxout, bloxin; kwargs...)
+    push!(bc.weights, w)
+
+    if haskey(kwargs, :sermon_rule)
+        xₒ = namespace_expr(bloxout.output, sys_out)
+        xᵢ = namespace_expr(bloxin.output, sys_in) #needed because this is also the θ term of the block receiving the connection
+        
+        # Custom values for the connection rule
+        f₀ = -0.780
+        f₁ = 0.198
+        f₂ = 0.302
+        f₃ = 0.851
+        f₄ = 0.998
+        x = xₒ - xᵢ
+        eq = sys_in.jcn ~ f₀ + f₁*cos(x) + f₂*sin(x) + f₃*cos(2*x) + f₄*sin(2*x)
+        accumulate_equation!(bc, eq)
+    end
+end
+

From 893a4d72b20eab52db729f8af1c9ea7fcc737fca Mon Sep 17 00:00:00 2001
From: agchesebro <76024790+agchesebro@users.noreply.github.com>
Date: Tue, 30 Jul 2024 16:28:13 -0400
Subject: [PATCH 06/10] Connector for DBS stimulator to Kuramoto oscillator

---
 src/blox/sermon_dbs.jl | 21 +++++++++++++++++----
 1 file changed, 17 insertions(+), 4 deletions(-)

diff --git a/src/blox/sermon_dbs.jl b/src/blox/sermon_dbs.jl
index 61eee082..d60ea8d1 100644
--- a/src/blox/sermon_dbs.jl
+++ b/src/blox/sermon_dbs.jl
@@ -3,9 +3,6 @@ using IfElse
 
 # TODO: List of missing components needed for this tutorial
 # 1. SuperBlox of Kuramoto oscillators to collect dynamics
-# 2. DBS stimulator block as implemented in the paper
-# 4. BloxConnector functionality for neurotransmitter pools
-# 5. BloxConnector functionality for DBS stimulator
 
 # Create neurotransmitter pools 
 struct SermonNPool <: NeuralMassBlox
@@ -30,7 +27,6 @@ struct SermonNPool <: NeuralMassBlox
     end
 end
 
-
 # Create a new DBS stimulator block
 # Largely based on MTK Standard Library smooth_square function in src/Blocks/sources.jl
 struct SermonDBS <: NeuralMassBlox
@@ -103,3 +99,20 @@ function (bc::BloxConnector)(
     end
 end
 
+function (bc::BloxConnector)(
+    bloxout::SermonDBS, 
+    bloxin::KuramotoOscillator; 
+    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)
+    xᵢ = namespace_expr(bloxin.output, sys_in) #needed because this is also the θ term of the block receiving the connection
+
+    eq = sys_in.jcn ~ w*xₒ*sin(xᵢ)
+    accumulate_equation!(bc, eq)
+end

From 159b33f94b13d77f91f21a95d12bdbf28744a86b Mon Sep 17 00:00:00 2001
From: agchesebro <76024790+agchesebro@users.noreply.github.com>
Date: Wed, 31 Jul 2024 09:40:57 -0400
Subject: [PATCH 07/10] Remove unused dependencies

---
 src/blox/sermon_dbs.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/blox/sermon_dbs.jl b/src/blox/sermon_dbs.jl
index d60ea8d1..25fd22bf 100644
--- a/src/blox/sermon_dbs.jl
+++ b/src/blox/sermon_dbs.jl
@@ -1,4 +1,3 @@
-using ModelingToolkitStandardLibrary.Blocks: smooth_square, smooth_step
 using IfElse
 
 # TODO: List of missing components needed for this tutorial

From e5394f1e1da2bf7a7c05ee124df98a5ce378a143 Mon Sep 17 00:00:00 2001
From: agchesebro <76024790+agchesebro@users.noreply.github.com>
Date: Wed, 31 Jul 2024 15:42:31 -0400
Subject: [PATCH 08/10] Tutorial initial draft - still missing global
 modulation

---
 src/blox/sermon_dbs_tutorial.jl | 51 +++++++++++++++++++++++++++++++++
 1 file changed, 51 insertions(+)
 create mode 100644 src/blox/sermon_dbs_tutorial.jl

diff --git a/src/blox/sermon_dbs_tutorial.jl b/src/blox/sermon_dbs_tutorial.jl
new file mode 100644
index 00000000..7d37e278
--- /dev/null
+++ b/src/blox/sermon_dbs_tutorial.jl
@@ -0,0 +1,51 @@
+using Neuroblox
+using DifferentialEquations
+using DataFrames
+using Distributions
+using Statistics
+using LinearAlgebra
+using Graphs
+using MetaGraphs
+using Random
+
+# First, create the three neurotransmitter pools. As they have the same structure,
+# we can use the same constructor with different parameters to create the systems.
+@named RRP = SermonNPool(τ_pool=1.0, p_pool=0.3)
+@named RP  = SermonNPool(τ_pool=10.0, p_pool=0.02)
+@named RtP = SermonNPool(τ_pool=50.0, p_pool=0.00035)
+
+# Next, create the DBS stimulator block. This block is a simple pulse with a given width
+# and period. 
+@named DBS = SermonDBS()
+
+# Next, create a collection of 50 Kuramotor oscillators with noise.
+# Set the number of oscillators
+N = 50
+
+# Define the natural distribution of oscillator frequencies
+Ω = 249
+σ = 26.317
+
+ks_blocks = [KuramotoOscillator(name=Symbol("KO$i"), 
+                                ω=rand(Normal(Ω, σ)),
+                                ζ=5.920,
+                                include_noise=true) for i in 1:N]
+
+# Create a graph and add all the oscillators to it
+g = MetaDiGraph()
+add_blox!.(Ref(g), ks_blocks)
+
+# Connect all oscillators to each other
+for i in 1:N
+    for j in 1:N
+        add_edge!(g, i, j, Dict(:weight => 1.0, :sermon_rule => true))
+    end
+end
+
+@named sys = system_from_graph(g)
+@time sys = structural_simplify(sys)
+
+
+sim_len = 150.0
+@time prob = SDEProblem(sys, [], (0.0, sim_len))
+@time sol = solve(prob, SRA1(), dt=0.0001, saveat=0.001)
\ No newline at end of file

From 34b6db4c1b196114a1413acdbb55e0590e491678 Mon Sep 17 00:00:00 2001
From: agchesebro <76024790+agchesebro@users.noreply.github.com>
Date: Wed, 7 Aug 2024 08:36:47 -0400
Subject: [PATCH 09/10] Updating tutorial to illustrate the issue with
 parameters in weight expressions

---
 src/blox/sermon_dbs.jl          | 14 +++++++++++--
 src/blox/sermon_dbs_tutorial.jl | 37 ++++++++++++++++++++++++++-------
 2 files changed, 41 insertions(+), 10 deletions(-)

diff --git a/src/blox/sermon_dbs.jl b/src/blox/sermon_dbs.jl
index 25fd22bf..d989bb75 100644
--- a/src/blox/sermon_dbs.jl
+++ b/src/blox/sermon_dbs.jl
@@ -37,7 +37,7 @@ struct SermonDBS <: NeuralMassBlox
     function SermonDBS(;
                         name,
                         namespace=nothing,
-                        f_stim=130,
+                        f_stim=130.0,
                         amplitude=1.0,
                         pulse_width=0.0001,
                         start_time=0.0      
@@ -83,6 +83,16 @@ function (bc::BloxConnector)(
     push!(bc.weights, w)
 
     if haskey(kwargs, :sermon_rule)
+        if haskey(kwargs, :extra_bloxs) && haskey(kwargs, :extra_params)
+            RRP, RP, RtP = kwargs[:extra_bloxs]
+            out_RRP = namespace_expr(RRP.output, get_namespaced_sys(RRP))
+            out_RP = namespace_expr(RP.output, get_namespaced_sys(RP))
+            out_RtP = namespace_expr(RtP.output, get_namespaced_sys(RtP))
+
+            M_RRP, M_RP, M_RtP, k_μ = kwargs[:extra_params]
+        else
+            error("Extra states and parameters need to be specified to use the Sermon rule")
+        end
         xₒ = namespace_expr(bloxout.output, sys_out)
         xᵢ = namespace_expr(bloxin.output, sys_in) #needed because this is also the θ term of the block receiving the connection
         
@@ -93,7 +103,7 @@ function (bc::BloxConnector)(
         f₃ = 0.851
         f₄ = 0.998
         x = xₒ - xᵢ
-        eq = sys_in.jcn ~ f₀ + f₁*cos(x) + f₂*sin(x) + f₃*cos(2*x) + f₄*sin(2*x)
+        eq = sys_in.jcn ~ w * max(M_RRP * out_RRP, M_RP * out_RP, M_RtP * out_RtP)*(f₀ + f₁*cos(x) + f₂*sin(x) + f₃*cos(2*x) + f₄*sin(2*x))
         accumulate_equation!(bc, eq)
     end
 end
diff --git a/src/blox/sermon_dbs_tutorial.jl b/src/blox/sermon_dbs_tutorial.jl
index 7d37e278..b681f0c3 100644
--- a/src/blox/sermon_dbs_tutorial.jl
+++ b/src/blox/sermon_dbs_tutorial.jl
@@ -7,6 +7,8 @@ using LinearAlgebra
 using Graphs
 using MetaGraphs
 using Random
+using DSP
+using ModelingToolkit: namespace_expr
 
 # First, create the three neurotransmitter pools. As they have the same structure,
 # we can use the same constructor with different parameters to create the systems.
@@ -16,7 +18,7 @@ using Random
 
 # Next, create the DBS stimulator block. This block is a simple pulse with a given width
 # and period. 
-@named DBS = SermonDBS()
+@named DBS = SermonDBS(start_time=50.0)
 
 # Next, create a collection of 50 Kuramotor oscillators with noise.
 # Set the number of oscillators
@@ -33,19 +35,38 @@ ks_blocks = [KuramotoOscillator(name=Symbol("KO$i"),
 
 # Create a graph and add all the oscillators to it
 g = MetaDiGraph()
-add_blox!.(Ref(g), ks_blocks)
+add_blox!.(Ref(g), vcat(ks_blocks, DBS, RRP, RP, RtP))
+
+# p = @parameters M_RRP=1.0 M_RP=0.5 M_RtP = 0.3 k_μ = 11976.0 
+
+# out_RRP = namespace_expr(RRP.output, get_namespaced_sys(RRP))
+# out_RP = namespace_expr(RP.output, get_namespaced_sys(RP))
+# out_RtP = namespace_expr(RtP.output, get_namespaced_sys(RtP))
+
+# # Connect all oscillators to each other
+# for i in 1:N
+#     for j in 1:N
+#         add_edge!(g, i, j, Dict(:weight => k_μ * max(M_RRP * out_RRP, M_RP * out_RP, M_RtP * out_RtP)/N, :sermon_rule => true))
+#     end
+# end
+
+p = @parameters M_RRP=1.0 M_RP=0.5 M_RtP = 0.3 k_μ = 11976.0 
 
 # Connect all oscillators to each other
 for i in 1:N
     for j in 1:N
-        add_edge!(g, i, j, Dict(:weight => 1.0, :sermon_rule => true))
+        add_edge!(g, i, j, Dict(:weight => k_μ/N, :sermon_rule => true, :extra_bloxs => [RRP, RP, RtP], :extra_params => p))
     end
 end
 
-@named sys = system_from_graph(g)
-@time sys = structural_simplify(sys)
+for i in 1:N
+    add_edge!(g, N+1, i, Dict(:weight => 72053.0))
+end
+
+add_edge!(g, N+1, N+2, Dict(:weight => 1.0))
+add_edge!(g, N+1, N+3, Dict(:weight => 1.0))
+add_edge!(g, N+1, N+4, Dict(:weight => 1.0))
 
+@named sys = system_from_graph(g, p)
+sys = structural_simplify(sys)
 
-sim_len = 150.0
-@time prob = SDEProblem(sys, [], (0.0, sim_len))
-@time sol = solve(prob, SRA1(), dt=0.0001, saveat=0.001)
\ No newline at end of file

From beb0acf1d46260af4d0c48790dbf574f414acbd5 Mon Sep 17 00:00:00 2001
From: agchesebro <76024790+agchesebro@users.noreply.github.com>
Date: Wed, 7 Aug 2024 19:37:16 -0400
Subject: [PATCH 10/10] What should be working tutorial code (but isn't)

---
 src/blox/sermon_dbs.jl          |  45 ++++++++++--
 src/blox/sermon_dbs_tutorial.jl | 124 ++++++++++++++++++++++++--------
 2 files changed, 135 insertions(+), 34 deletions(-)

diff --git a/src/blox/sermon_dbs.jl b/src/blox/sermon_dbs.jl
index d989bb75..4725536a 100644
--- a/src/blox/sermon_dbs.jl
+++ b/src/blox/sermon_dbs.jl
@@ -1,9 +1,7 @@
 using IfElse
 
-# TODO: List of missing components needed for this tutorial
-# 1. SuperBlox of Kuramoto oscillators to collect dynamics
-
-# Create neurotransmitter pools 
+# Neurotransmitter pool block 
+# Implements Equation 4
 struct SermonNPool <: NeuralMassBlox
     params
     output
@@ -27,7 +25,8 @@ struct SermonNPool <: NeuralMassBlox
 end
 
 # Create a new DBS stimulator block
-# Largely based on MTK Standard Library smooth_square function in src/Blocks/sources.jl
+# Largely based on MTK Standard Library pulse block from the digital electronics
+# Implements unnumbered equation below equation 2 for the pulse train, with default parameters from the paper
 struct SermonDBS <: NeuralMassBlox
     params
     output
@@ -52,6 +51,7 @@ struct SermonDBS <: NeuralMassBlox
     end
 end
 
+# Connects the DBS stimulator to the neurotransmitter pool
 function (bc::BloxConnector)(
     bloxout::SermonDBS, 
     bloxin::SermonNPool; 
@@ -89,7 +89,7 @@ function (bc::BloxConnector)(
             out_RP = namespace_expr(RP.output, get_namespaced_sys(RP))
             out_RtP = namespace_expr(RtP.output, get_namespaced_sys(RtP))
 
-            M_RRP, M_RP, M_RtP, k_μ = kwargs[:extra_params]
+            M_RRP, M_RP, M_RtP, k_μ, I = kwargs[:extra_params]
         else
             error("Extra states and parameters need to be specified to use the Sermon rule")
         end
@@ -97,17 +97,20 @@ function (bc::BloxConnector)(
         xᵢ = namespace_expr(bloxin.output, sys_in) #needed because this is also the θ term of the block receiving the connection
         
         # Custom values for the connection rule
+        # Values from supplementary material Table A
         f₀ = -0.780
         f₁ = 0.198
         f₂ = 0.302
         f₃ = 0.851
         f₄ = 0.998
         x = xₒ - xᵢ
+        # Equation 6 from the paper
         eq = sys_in.jcn ~ w * max(M_RRP * out_RRP, M_RP * out_RP, M_RtP * out_RtP)*(f₀ + f₁*cos(x) + f₂*sin(x) + f₃*cos(2*x) + f₄*sin(2*x))
         accumulate_equation!(bc, eq)
     end
 end
 
+# Connects the DBS stimulator to the Kuramoto oscillators (Iₜ term in equation 1)
 function (bc::BloxConnector)(
     bloxout::SermonDBS, 
     bloxin::KuramotoOscillator; 
@@ -125,3 +128,33 @@ function (bc::BloxConnector)(
     eq = sys_in.jcn ~ w*xₒ*sin(xᵢ)
     accumulate_equation!(bc, eq)
 end
+
+# Debugging purposes only
+# Connects the Kuramoto oscillators without neurotransmitter modulation
+function (bc::BloxConnector)(
+    bloxout::KuramotoOscillator, 
+    bloxin::KuramotoOscillator; 
+    kwargs...
+)
+    sys_out = get_namespaced_sys(bloxout)
+    sys_in = get_namespaced_sys(bloxin)
+
+    #technically these two lines aren't needed, but useful to have if weighting occurs here
+    w = generate_weight_param(bloxout, bloxin; kwargs...)
+    push!(bc.weights, w)
+
+    if haskey(kwargs, :sermon_rule_no_neurotransmitters)
+        xₒ = namespace_expr(bloxout.output, sys_out)
+        xᵢ = namespace_expr(bloxin.output, sys_in) #needed because this is also the θ term of the block receiving the connection
+        
+        # Custom values for the connection rule
+        f₀ = -0.780
+        f₁ = 0.198
+        f₂ = 0.302
+        f₃ = 0.851
+        f₄ = 0.998
+        x = xₒ - xᵢ
+        eq = sys_in.jcn ~ w *(f₀ + f₁*cos(x) + f₂*sin(x) + f₃*cos(2*x) + f₄*sin(2*x))
+        accumulate_equation!(bc, eq)
+    end
+end
\ No newline at end of file
diff --git a/src/blox/sermon_dbs_tutorial.jl b/src/blox/sermon_dbs_tutorial.jl
index b681f0c3..ed17bcc8 100644
--- a/src/blox/sermon_dbs_tutorial.jl
+++ b/src/blox/sermon_dbs_tutorial.jl
@@ -1,30 +1,29 @@
 using Neuroblox
 using DifferentialEquations
-using DataFrames
-using Distributions
-using Statistics
-using LinearAlgebra
-using Graphs
-using MetaGraphs
-using Random
+#using DataFrames
+using Distributions, Random
+#using Statistics
+#using LinearAlgebra
+using Graphs, MetaGraphs
 using DSP
-using ModelingToolkit: namespace_expr
 
 # First, create the three neurotransmitter pools. As they have the same structure,
 # we can use the same constructor with different parameters to create the systems.
+# Parameters are taken from the supplementary material Table A.
 @named RRP = SermonNPool(τ_pool=1.0, p_pool=0.3)
 @named RP  = SermonNPool(τ_pool=10.0, p_pool=0.02)
 @named RtP = SermonNPool(τ_pool=50.0, p_pool=0.00035)
 
 # Next, create the DBS stimulator block. This block is a simple pulse with a given width
-# and period. 
-@named DBS = SermonDBS(start_time=50.0)
+# and period. To try and recreate Fig. 2A, start stimulation at 50s.
+@named DBS = SermonDBS(start_time=50.0, pulse_width=0.0001)
 
 # Next, create a collection of 50 Kuramotor oscillators with noise.
-# Set the number of oscillators
+# Set the number of oscillators (default from the paper)
 N = 50
 
 # Define the natural distribution of oscillator frequencies
+# Parameters from supplementary material Table A
 Ω = 249
 σ = 26.317
 
@@ -33,36 +32,27 @@ ks_blocks = [KuramotoOscillator(name=Symbol("KO$i"),
                                 ζ=5.920,
                                 include_noise=true) for i in 1:N]
 
-# Create a graph and add all the oscillators to it
+# Create a graph and add all the blocks as nodes
 g = MetaDiGraph()
 add_blox!.(Ref(g), vcat(ks_blocks, DBS, RRP, RP, RtP))
 
-# p = @parameters M_RRP=1.0 M_RP=0.5 M_RtP = 0.3 k_μ = 11976.0 
-
-# out_RRP = namespace_expr(RRP.output, get_namespaced_sys(RRP))
-# out_RP = namespace_expr(RP.output, get_namespaced_sys(RP))
-# out_RtP = namespace_expr(RtP.output, get_namespaced_sys(RtP))
-
-# # Connect all oscillators to each other
-# for i in 1:N
-#     for j in 1:N
-#         add_edge!(g, i, j, Dict(:weight => k_μ * max(M_RRP * out_RRP, M_RP * out_RP, M_RtP * out_RtP)/N, :sermon_rule => true))
-#     end
-# end
-
-p = @parameters M_RRP=1.0 M_RP=0.5 M_RtP = 0.3 k_μ = 11976.0 
+# Additional connection parameters specified from supplementary material Table A
+p = @parameters M_RRP=1.0 M_RP=0.5 M_RtP = 0.3 k_μ = 11976.0 I =72053.0
 
 # Connect all oscillators to each other
 for i in 1:N
     for j in 1:N
-        add_edge!(g, i, j, Dict(:weight => k_μ/N, :sermon_rule => true, :extra_bloxs => [RRP, RP, RtP], :extra_params => p))
+        # Weight is k_μ to be multiplied by the other terms in equation 6
+        add_edge!(g, i, j, Dict(:weight => k_μ, :sermon_rule => true, :extra_bloxs => [RRP, RP, RtP], :extra_params => p))
     end
 end
 
+# Connect DBS stimulator to oscillators
 for i in 1:N
-    add_edge!(g, N+1, i, Dict(:weight => 72053.0))
+    add_edge!(g, N+1, i, Dict(:weight => I))
 end
 
+# Connect DBS stimulator to neurotransmitter pools
 add_edge!(g, N+1, N+2, Dict(:weight => 1.0))
 add_edge!(g, N+1, N+3, Dict(:weight => 1.0))
 add_edge!(g, N+1, N+4, Dict(:weight => 1.0))
@@ -70,3 +60,81 @@ add_edge!(g, N+1, N+4, Dict(:weight => 1.0))
 @named sys = system_from_graph(g, p)
 sys = structural_simplify(sys)
 
+# Simulate for the length of Fig. 2A
+sim_len = 150.0
+@time prob = SDEProblem(sys, [], (0.0, sim_len))
+# EulerHeun because you need to force timestops at the pulse width of the DBS stimulator
+# so fixed timestep solver is the easiest way to do this. Maxiter problems can happen with 
+# the adaptive ones I've tried.
+@time sol = solve(prob, EulerHeun(), dt=0.0001, saveat=0.001)
+
+# Helper function because θ from the Kuramoto oscillators is unbounded and needs wrapping
+# to get the correct spectrogram
+function wrapTo2Pi(x)
+    posinput = x .> 0
+    wrapped = mod.(x, 2*π)
+    wrapped[wrapped .== 0 .&& posinput] .= 2*π
+    return wrapped
+end
+
+# Helper code to plot the spectrogram averaged acrossed the 50 oscillators
+data = Array(sol)
+spec = spectrogram(wrapTo2Pi(data[1, :]), div(150001, 200); fs=1000)
+
+freq = spec.freq
+time = spec.time
+hmmpower = spec.power
+
+for i = 2:50
+    spec = spectrogram(wrapTo2Pi(data[i, :]), div(150001, 200); fs=1000)
+    hmmpower .+= spec.power
+end
+
+using Plots
+# This should reproduce Fig. 2A. Instead, what I'm seeing is oscillating around 
+# the mean value (~20 Hz) as in the paper for the first 50s, then a jump up to a mean
+# around the stimulation frequency (130Hz) but no emergent coherence at the 300Hz range.
+heatmap(time, freq, pow2db.(hmmpower))
+
+# To check if the coupling matches the values in Fig. 2D, compute equation 5 from the three
+# neurotransmitter pools.
+kₜ = zeros(length(sol))
+for i = 1:length(sol)
+    kₜ[i] = max(data[N+1, i], data[N+2, i]*0.5, data[N+3, i]*0.3)
+end
+
+# This should reproduce the blue line in Fig. 2D. The shape is roughly accurate, but the 
+# scale is entirely off. I think this is because the authors dropped a factor in connecting
+# the DBS stimulator to the neurotransmitter pools. Running just that simulation shows that
+# they're very far off the actual values while preserving the shape of the dynamics.
+plot(sol.t, kₜ .* 11976.0)
+
+# Transmitter pool only simulation
+# This *should* reproduce Figure 1B, but it doesn't. That's why I'm convinced there's a term
+# missing in equation 4. The shape is correct, but the scale is off.
+
+# Setup neurotransmitter pools
+@named RRP = SermonNPool(τ_pool=1.0, p_pool=0.3)
+@named RP  = SermonNPool(τ_pool=10.0, p_pool=0.02)
+@named RtP = SermonNPool(τ_pool=50.0, p_pool=0.00035)
+
+# Next, create the DBS stimulator block. This block is a simple pulse with a given width
+# and period. 
+@named DBS = SermonDBS(start_time=50.0, pulse_width=0.0001)
+
+p = @parameters M_RRP=1.0 M_RP=0.5 M_RtP = 0.3 k_μ = 11976.0 
+
+g = MetaDiGraph()
+add_blox!.(Ref(g), vcat(DBS, RRP, RP, RtP))
+add_edge!(g, 1, 2, Dict(:weight => 1.0))
+add_edge!(g, 1, 3, Dict(:weight => 1.0))
+add_edge!(g, 1, 4, Dict(:weight => 1.0))
+
+@named sys = system_from_graph(g, p)
+sys = structural_simplify(sys)
+
+prob = ODEProblem(sys, [], (0.0, 150.0))
+sol = solve(prob, Euler(), dt=0.0001, saveat=0.001)
+
+# This should reproduce Fig. 1B. It doesn't :()
+plot(sol)