diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index b347e12f..e1a42359 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -14,26 +14,22 @@ jobs:
         group:
           - Basics
           - Advanced
-          - GraphDynamics
+          - GraphDynamics1
+          - GraphDynamics2
+          - GraphDynamics3
         version:
           - '1'
     steps:
-      - uses: actions/checkout@v3
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
-      - uses: actions/cache@v3
-        env:
-          cache-name: cache-artifacts
-        with:
-          path: ~/.julia/artifacts
-          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
-          restore-keys: |
-            ${{ runner.os }}-test-${{ env.cache-name }}-
-            ${{ runner.os }}-test-
-            ${{ runner.os }}-
+      - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
+        env:
+          JULIA_NUM_THREADS: 4
+          GROUP: ${{ matrix.group }}
       - uses: julia-actions/julia-processcoverage@v1
       - uses: codecov/codecov-action@v3
         with:
diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml
index cde17110..3465aef6 100644
--- a/.github/workflows/Documentation.yml
+++ b/.github/workflows/Documentation.yml
@@ -16,6 +16,7 @@ jobs:
       - uses: julia-actions/setup-julia@latest
         with:
           version: '1'
+      - uses: julia-actions/cache@v2
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy
diff --git a/test/GraphDynamicsTests/runtests.jl b/test/GraphDynamicsTests/runtests.jl
index 140ca507..0d957460 100644
--- a/test/GraphDynamicsTests/runtests.jl
+++ b/test/GraphDynamicsTests/runtests.jl
@@ -1,16 +1,27 @@
 include(joinpath(@__DIR__(), "test_suite.jl"))
 
-basic_smoketest()
-neuron_and_neural_mass_comparison_tests()
-basic_hh_network_tests()
-stochastic_hh_network_tests()
-ngei_test()
-kuramato_test()
-wta_tests()
-dbs_circuit_components()
-dbs_circuit()
-discrete()
-striatum_tests()
-lif_exci_inh_tests()
-decision_making_test()
-ping_tests()
+const GROUP = get(ENV, "GROUP", "All")
+@show GROUP
+
+if GROUP == "All" || GROUP == "GraphDynamics1"
+    basic_smoketest()
+    neuron_and_neural_mass_comparison_tests()
+    basic_hh_network_tests()
+    stochastic_hh_network_tests()
+    ngei_test()
+end
+
+if GROUP == "All" || GROUP == "GraphDynamics2"
+    kuramato_test()
+    wta_tests()
+    dbs_circuit_components()
+    dbs_circuit()
+    discrete()
+    striatum_tests()
+end
+
+if GROUP == "All" || GROUP == "GraphDynamics3"
+    lif_exci_inh_tests()
+    decision_making_test()
+    ping_tests()
+end
diff --git a/test/bayesian_fitting.jl b/test/bayesian_fitting.jl
deleted file mode 100644
index 2274d154..00000000
--- a/test/bayesian_fitting.jl
+++ /dev/null
@@ -1,69 +0,0 @@
-using Neuroblox, StochasticDiffEq ,ForwardDiff, LinearAlgebra, Test
-using Turing
-using Distributions
-using Turing: Variational
-using Random
-using Tracker
-
-@named VdP = van_der_pol()
-
-prob_vdp = SDEProblem(VdP,[0.1,0.1],[0.0, 10.0],[])
-sol = solve(prob_vdp,EM(),dt=0.1)
-time = sol.t
-dt = time[2]-time[1]
-x = sol[1,:]
-y = sol[2,:]
-
-# turing model for fitting SDE without external noise
-@model function fitvpEM(datax, datay)
-  θ ~ Gamma(0.1,10.0)
-  ϕ ~ Uniform(0,0.5)
-  dxh = datay[1:end-1]
-  dyh = θ .* (1 .- datax[1:end-1] .^ 2) .* datay[1:end-1] .- datax[1:end-1]
-  datax[2:end] ~ MvNormal(datax[1:end-1] .+ dxh * dt,ϕ^2*dt * I)
-  datay[2:end] ~ MvNormal(datay[1:end-1] .+ dyh * dt,ϕ^2*dt * I)
-end
-
-Turing.setadbackend(:forwarddiff)
-modelEM = fitvpEM(x,y)
-chain = Turing.sample(modelEM,NUTS(0.65),200)
-@test 0.5<mean(chain[:θ])<1.5
-@test 0.05<mean(chain[:ϕ])<0.15
-
-#turing model for fitting with external noise
-noise = Normal(0,0.3)
-xn = x .+ rand(noise,length(x))
-yn = y .+ rand(noise,length(y))
-
-@model function fitvpEMn(datax, datay)
-  σ = 0.3 # here we explicitly set the noise, but it could also be a parameter
-  θ ~ Gamma(0.25,4.0)
-  ϕ ~ Uniform(0,0.2)
-  # create arrays for hidden variables
-  xh = eltype(θ).(zero(datax))
-  yh = eltype(θ).(zero(datay))
-  xh[1] = 0.1
-  yh[1] = 0.1
-  for i in 2:length(datax)
-      dxh = yh[i-1]
-      dyh = θ * (1.0 - xh[i-1]^2) * yh[i-1] - xh[i-1]
-      xh[i] ~ Normal(xh[i-1] + dxh * dt,ϕ*sqrt(dt))
-      yh[i] ~ Normal(yh[i-1] + dyh * dt,ϕ*sqrt(dt))
-  end
-  datax ~ MvNormal(xh,σ ^ 2 * I)
-  datay ~ MvNormal(yh,σ ^ 2 * I)
-end
-
-# ADVI
-modelEMn = fitvpEMn(xn,yn)
-Turing.setadbackend(:tracker)
-advi = ADVI(10, 1000)
-setchunksize(8)
-q = vi(modelEMn, advi)
-
-# sampling
-z = rand(q, 500)
-avg = vec(mean(z; dims = 2))
-
-@test 0.01<avg[1]<3.0 # testing θ
-@test 0.05<avg[2]<0.25 # tesing ϕ
diff --git a/test/fitting.jl b/test/fitting.jl
deleted file mode 100644
index 92033a6f..00000000
--- a/test/fitting.jl
+++ /dev/null
@@ -1,66 +0,0 @@
-using Neuroblox, LinearAlgebra, OrdinaryDiffEq, Optimization, OptimizationOptimJL, OptimizationOptimisers, Flux, ForwardDiff, Test
-
-@named Str = jansen_ritSC(τ=0.0022, H=20, λ=300, r=0.3)
-@named GPe = jansen_ritSC(τ=0.04, H=20, λ=400, r=0.1)
-@named STN = jansen_ritSC(τ=0.01, H=20, λ=500, r=0.1)
-@named GPi = jansen_ritSC(τ=0.014, H=20, λ=400, r=0.1)
-
-# Connect Regions through Adjacency Matrix
-blox = [Str, GPe, STN, GPi]
-sys = [s.odesystem for s in blox]
-connect = [s.connector for s in blox]
-
-@parameters C_Cor=3 C_BG_Th=3 C_Cor_BG_Th=9.75 C_BG_Th_Cor=9.75
-adj_matrix = [0 0 0 0;
-             -0.5*C_BG_Th -0.5*C_BG_Th C_BG_Th 0;
-              0 -0.5*C_BG_Th 0 0;
-              0 -0.5*C_BG_Th C_BG_Th 0]
-
-@named BG_Circuit = LinearConnections(sys=sys, adj_matrix=adj_matrix, connector=connect)
-
-sim_dur = 2.0 # Simulation time (seconds)
-prob = ODEProblem(structural_simplify(BG_Circuit), [], (0.0, sim_dur), [])
-
-# Fitting Procedure
-# Grab system parameters
-ptrue = prob.p
-# Alter system parameters randomly
-pinit = ptrue .+ ptrue .* randn(length(prob.p)) ./ 100
-# Solve problem as normal
-data = solve(prob, Rodas4(), abstol=1e-9, reltol=1e-9, saveat=0.05)
-
-# Define loss function
-lb = 0.1 .* ptrue
-ub = 2 .* ptrue
-internalnorm(u,t) = sqrt(ForwardDiff.value(sum(abs2,u)))/length(u)
-function loss(p,_)
-    sol = solve(remake(prob,p=p), Rodas4(), abstol=1e-9, reltol=1e-9, saveat=0.05)
-    if sol.retcode != :Success
-        Inf * abs(p[1])
-    else
-        sqrt(sum(abs2, sol .- data)/length(sol))
-    end
-end
-
-iter = 0
-function cb(p,l)
-    global iter
-    iter += 1
-    @show iter,l
-    l < 0.01
-end
-
-optfun = OptimizationFunction(loss,Optimization.AutoForwardDiff())
-optprob = OptimizationProblem(optfun,pinit)
-
-# loss(pinit,nothing)
-# loss(ptrue,nothing)
-
-# Solve optimization problem
-optsol = Optimization.solve(optprob, OptimizationOptimisers.Adam(0.0001),maxiters = 200, callback = cb)
-optprob2 = remake(optprob,u0 = optsol.u)
-optsol2 = Optimization.solve(optprob2, BFGS(initial_stepnorm = 0.0001),maxiters = 200, callback = cb)
-
-# Calculate squared error of fit
-@test norm(optsol.u - ptrue) < 20 # balance threshold value and maxiters for run time
-#@test_broken !iszero(norm(optsol.u - pinit))
diff --git a/test/graph_to_dataframe.jl b/test/graph_to_dataframe.jl
deleted file mode 100644
index 4dd55c8d..00000000
--- a/test/graph_to_dataframe.jl
+++ /dev/null
@@ -1,70 +0,0 @@
-using Neuroblox, Graphs, MetaGraphs, DataFrames, Test
-
-"""
-harmonic oscillator neural mass blox
-"""
-@named GPe = harmonic_oscillator(ω=4*2*π, ζ=1, k=(4*2*π)^2, h=0.1)
-@named STN = harmonic_oscillator(ω=4*2*π, ζ=1, k=(4*2*π)^2, h=1)
-
-# Connect Regions through Adjacency Matrix
-@parameters g_GPe_STN=1.0 g_STN_GPe=1.0
-
-#create graph
-g = MetaDiGraph()
-add_blox!(g,GPe)
-add_blox!(g,STN)
-add_edge!(g,1,1,:weight,1.0)
-add_edge!(g,1,2,:weight,g_STN_GPe)
-add_edge!(g,2,1,:weight,g_STN_GPe*g_GPe_STN)
-add_edge!(g,2,2,:weight,1.0)
-
-@named two_regions_gr = ODEfromGraph(g)
-
-sim_dur = 5.0 # Simulate for 10 Seconds
-
-# returns dataframe with time series for 4 outputs
-sol = simulate(structural_simplify(two_regions_gr), [], (0.0, sim_dur), [])
-@test typeof(sol) == DataFrame
-
-"""
-jansen rit neural mass blox
-"""
-
-@named Str = jansen_ritSC(τ=0.0022, H=20, λ=300, r=0.3)
-@named GPe = jansen_ritSC(τ=0.04, H=20, λ=400, r=0.1)
-@named STN = jansen_ritSC(τ=0.01, H=20, λ=500, r=0.1)
-@named GPi = jansen_ritSC(τ=0.014, H=20, λ=400, r=0.1)
-
-# Connect Regions through Adjacency Matrix
-@parameters C_Cor=3 C_BG_Th=3 C_Cor_BG_Th=9.75 C_BG_Th_Cor=9.75
-
-# Create Graph
-g = MetaDiGraph()
-add_blox!(g,Str)
-add_blox!(g,GPe)
-add_blox!(g,STN)
-add_blox!(g,GPi)
-add_edge!(g,1,1,:weight,0.0)
-add_edge!(g,1,2,:weight,0.0)
-add_edge!(g,1,3,:weight,0.0)
-add_edge!(g,1,4,:weight,0.0)
-add_edge!(g,2,1,:weight,-0.5*C_BG_Th)
-add_edge!(g,2,2,:weight,-0.5*C_BG_Th)
-add_edge!(g,2,3,:weight,C_BG_Th)
-add_edge!(g,2,4,:weight,0.0)
-add_edge!(g,3,1,:weight,0.0)
-add_edge!(g,3,2,:weight,-0.5*C_BG_Th)
-add_edge!(g,3,3,:weight,0.0)
-add_edge!(g,3,4,:weight,0.0)
-add_edge!(g,4,1,:weight,0.0)
-add_edge!(g,4,2,:weight,-0.5*C_BG_Th)
-add_edge!(g,4,3,:weight, C_BG_Th)
-add_edge!(g,4,4,:weight,0.0)
-
-@named four_regions_gr = ODEfromGraph(g)
-
-sim_dur = 5.0 # Simulate for 5 Seconds
-
-# returns dataframe with time series for 2*4 outputs
-sol = simulate(structural_simplify(four_regions_gr), [], (0.0, sim_dur), [])
-@test typeof(sol) == DataFrame
diff --git a/test/jansen_rit_component_tests.jl b/test/jansen_rit_component_tests.jl
deleted file mode 100644
index f06309ac..00000000
--- a/test/jansen_rit_component_tests.jl
+++ /dev/null
@@ -1,125 +0,0 @@
-using Neuroblox, DifferentialEquations, DataFrames, Test, Distributions, Statistics, LinearAlgebra, Graphs, MetaGraphs, Random
-
-τ_factor = 1000 #needed because the old units were in seconds, and we need ms to be consistent
-
-# Create Regions
-@named Str = jansen_ritC(τ=0.0022*τ_factor, H=20, λ=300, r=0.3)
-@named gpe = jansen_ritC(τ=0.04*τ_factor, H=20, λ=400, r=0.1)
-@named stn = jansen_ritC(τ=0.01*τ_factor, H=20, λ=500, r=0.1)
-@named gpi = jansen_ritSC(τ=0.014*τ_factor, H=20, λ=400, r=0.1)
-@named Th  = jansen_ritSC(τ=0.002*τ_factor, H=10, λ=20, r=5)
-@named EI  = jansen_ritSC(τ=0.01*τ_factor, H=20, λ=5, r=5)
-@named PY  = jansen_ritSC(τ=0.001*τ_factor, H=20, λ=5, r=0.15)
-@named II  = jansen_ritSC(τ=2.0*τ_factor, H=60, λ=5, r=5)
-
-# Connect Regions through Adjacency Matrix
-blox = [Str, gpe, stn, gpi, Th, EI, PY, II]
-sys = [s.odesystem for s in blox]
-connect = [s.connector for s in blox]
-
-@parameters C_Cor=60 C_BG_Th=60 C_Cor_BG_Th=5 C_BG_Th_Cor=5
-
-adj_matrix_lin = [0 0 0 0 0 0 0 0;
-            -0.5*C_BG_Th -0.5*C_BG_Th C_BG_Th 0 0 0 0 0;
-            0 -0.5*C_BG_Th 0 0 0 0 C_Cor_BG_Th 0;
-            0 -0.5*C_BG_Th C_BG_Th 0 0 0 0 0;
-            0 0 0 -0.5*C_BG_Th 0 0 0 0;
-            0 0 0 0 C_BG_Th_Cor 0 6*C_Cor 0;
-            0 0 0 0 0 4.8*C_Cor 0 -1.5*C_Cor;
-            0 0 0 0 0 0 1.5*C_Cor 3.3*C_Cor]
-
-@named CBGTC_Circuit_lin = LinearConnections(sys=sys, adj_matrix=adj_matrix_lin, connector=connect)
-
-sim_dur = 2000.0 # Simulate for 2 seconds
-mysys = structural_simplify(CBGTC_Circuit_lin)
-sol = simulate(mysys, [], (0.0, sim_dur), [], Vern7(); saveat=1)
-@test sol[!, "gpi₊x(t)"][4] ≈ -2219.2560209502685 #updated to new value in ms
-
-"""
-Testing new Jansen-Rit blox
-
-This sets up the exact same system as above, but using the JansenRit with different flags cortical/subcortical flags.
-
-The purpose of this test is to make sure that setting up everything using System rather than ODESystem works as expected.
-It also shows test code for the new system_from_graph calls and handling delays when everything is 0 (MethodOfSteps returns the 
-same thing as the old simulate call with AutoVern7(Rodas4() since there are no delays.)
-"""
-
-# test new Jansen-Rit blox
-@named Str = JansenRit(τ=0.0022*τ_factor, H=20, λ=300, r=0.3)
-@named gpe = JansenRit(τ=0.04*τ_factor, cortical=false) # all default subcortical except τ
-@named stn = JansenRit(τ=0.01*τ_factor, H=20, λ=500, r=0.1)
-@named gpi = JansenRit(cortical=false) # default parameters subcortical Jansen Rit blox
-@named Th  = JansenRit(τ=0.002*τ_factor, H=10, λ=20, r=5)
-@named EI  = JansenRit(τ=0.01*τ_factor, H=20, λ=5, r=5)
-@named PY  = JansenRit(cortical=true) # default parameters cortical Jansen Rit blox
-@named II  = JansenRit(τ=2.0*τ_factor, H=60, λ=5, r=5)
-blox = [Str, gpe, stn, gpi, Th, EI, PY, II]
-
-# test graphs
-g = MetaDiGraph()
-add_blox!.(Ref(g), blox)
-
-# Store parameters to be passed later on
-params = @parameters C_Cor=60 C_BG_Th=60 C_Cor_BG_Th=5 C_BG_Th_Cor=5
-
-add_edge!(g, 2, 1, Dict(:weight => -0.5*C_BG_Th))
-add_edge!(g, 2, 2, Dict(:weight => -0.5*C_BG_Th))
-add_edge!(g, 2, 3, Dict(:weight => C_BG_Th))
-add_edge!(g, 3, 2, Dict(:weight => -0.5*C_BG_Th))
-add_edge!(g, 3, 7, Dict(:weight => C_Cor_BG_Th))
-add_edge!(g, 4, 2, Dict(:weight => -0.5*C_BG_Th))
-add_edge!(g, 4, 3, Dict(:weight => C_BG_Th))
-add_edge!(g, 5, 4, Dict(:weight => -0.5*C_BG_Th))
-add_edge!(g, 6, 5, Dict(:weight => C_BG_Th_Cor))
-add_edge!(g, 6, 7, Dict(:weight => 6*C_Cor))
-add_edge!(g, 7, 6, Dict(:weight => 4.8*C_Cor))
-add_edge!(g, 7, 8, Dict(:weight => -1.5*C_Cor))
-add_edge!(g, 8, 7, Dict(:weight => 1.5*C_Cor))
-add_edge!(g, 8, 8, Dict(:weight => 3.3*C_Cor))
-
-# add_edge!(g, 2, 1, Dict(:weight => -0.5*60, :delay => 0.01))
-# add_edge!(g, 2, 2, Dict(:weight => -0.5*60, :delay => 0.01))
-# add_edge!(g, 2, 3, Dict(:weight => 60, :delay => 0.01))
-# add_edge!(g, 3, 2, Dict(:weight => -0.5*60, :delay => 0.01))
-# add_edge!(g, 3, 7, Dict(:weight => 5, :delay => 0.01))
-# add_edge!(g, 4, 2, Dict(:weight => -0.5*60, :delay => 0.01))
-# add_edge!(g, 4, 3, Dict(:weight => 60, :delay => 0.01))
-# add_edge!(g, 5, 4, Dict(:weight => -0.5*60, :delay => 0.01))
-# add_edge!(g, 6, 5, Dict(:weight => 5, :delay => 0.01))
-# add_edge!(g, 6, 7, Dict(:weight => 6*60, :delay => 0.01))
-# add_edge!(g, 7, 6, Dict(:weight => 4.8*60, :delay => 0.01))
-# add_edge!(g, 7, 8, Dict(:weight => -1.5*60, :delay => 0.01))
-# add_edge!(g, 8, 7, Dict(:weight => 1.5*60, :delay => 0.01))
-# add_edge!(g, 8, 8, Dict(:weight => 3.3*60, :delay => 0.01))
-
-@named final_system = system_from_graph(g, params)
-final_delays = graph_delays(g)
-sim_dur = 2000.0 # Simulate for 2 Seconds
-final_system_sys = structural_simplify(final_system)
-prob = ODEProblem(final_system_sys,
-    [],
-    (0.0, sim_dur))
-alg = Vern7()
-sol_dde_no_delays = solve(prob, alg, saveat=1)
-sol2 = DataFrame(sol_dde_no_delays)
-@test isapprox(sol2[!, "gpi₊x(t)"][500:1000], sol[!, "gpi₊x(t)"][500:1000], rtol=1e-8)
-
-
-# Alternative version using adjacency matrix
-g2 = MetaDiGraph()
-add_blox!.(Ref(g2), blox)
-create_adjacency_edges!(g2, adj_matrix_lin)
-
-@named final_system = system_from_graph(g2, params)
-final_delays = graph_delays(g2)
-sim_dur = 2000.0 # Simulate for 2 Seconds
-final_system_sys = structural_simplify(final_system)
-prob = DDEProblem(final_system_sys,
-    [],
-    (0.0, sim_dur),
-    constant_lags = final_delays)
-alg = MethodOfSteps(Vern7())
-sol_dde_no_delays = solve(prob, alg, saveat=1)
-sol3 = DataFrame(sol_dde_no_delays)
-@test isapprox(sol3[!, "gpi₊x(t)"][500:1000], sol[!, "gpi₊x(t)"][500:1000], rtol=1e-8)
\ No newline at end of file
diff --git a/test/new_LB_blox.jl b/test/new_LB_blox.jl
deleted file mode 100644
index 717b24d5..00000000
--- a/test/new_LB_blox.jl
+++ /dev/null
@@ -1,59 +0,0 @@
-using Plots, Neuroblox, OrdinaryDiffEq, MetaGraphs, Graphs, LinearAlgebra, Random, DataFrames, Test
-
-
-function LinearConnectionsLB(;name, sys=sys, adj_matrix=adj_matrix, connector=connector)
-    adj = adj_matrix .* connector
-    eqs = []
-    for region_num in 1:length(sys)
-		norm_factor = sum(adj_matrix[:, region_num])
-    	push!(eqs, sys[region_num].jcn ~ sum(adj[:, region_num])/norm_factor)
-    end
-    return @named Circuit = ODESystem(eqs, systems = sys)
-end
-
-adj = [0 1 0 0 0;
-       0 0 1 0 0;
-       0 0 0 1 0;
-       0 0 0 0 1;
-       1 0 0 0 0]
-
-blox = []
-
-for i = 1:5
-    lb = LarterBreakspearBlox(name=Symbol("LB$i"))
-    push!(blox,lb)
-end
-
-sys = [b.odesystem for b in blox]
-con = [s.connector for s in blox]
-@named LB_circuit_lin = LinearConnectionsLB(sys=sys, adj_matrix=adj,connector=con)
-mysys = structural_simplify(LB_circuit_lin)
-
-uw = rand(Int(length(mysys.states)/3))
-uv = 0.1*rand(Int(length(mysys.states)/3)) .- 0.05
-uz = 0.1*rand(Int(length(mysys.states)/3)) .- 0.05
-u0 = collect(Iterators.flatten(zip(uv,uz,uw)))
-
-prob = ODEProblem(mysys,u0,(0.0,5e2),[])
-sol = solve(prob,AutoVern7(Rodas4()),saveat=0.1)
-dsol = DataFrame(sol)
-
-blox2 = []
-
-for i = 1:5
-    lb = LarterBreakspear(name=Symbol("LB$i"))
-    push!(blox2,lb)
-end
-
-adj_temp = adj
-
-g = MetaDiGraph()
-add_blox!.(Ref(g), blox2)
-create_adjacency_edges!(g, adj_temp)
-@named mysys2 = system_from_graph(g)
-mysys2 = structural_simplify(mysys2)
-prob2 = ODEProblem(mysys2,u0,(0.0,5e2),[])
-sol2 = solve(prob2,AutoVern7(Rodas4()),saveat=0.1)
-dsol2 = DataFrame(sol2)
-
-@test isapprox(dsol2[!, "LB2₊V(t)"][100:end], dsol[!, "LB2₊V(t)"][100:end], rtol=1e-8)
\ No newline at end of file
diff --git a/test/ode_from_graph.jl b/test/ode_from_graph.jl
deleted file mode 100644
index e2f39473..00000000
--- a/test/ode_from_graph.jl
+++ /dev/null
@@ -1,63 +0,0 @@
-using Neuroblox, Test, SparseArrays, Graphs, MetaGraphs
-
-@named GPe = HarmonicOscillatorBlox()
-@named STN = HarmonicOscillatorBlox()
-@named BP = BandPassFilterBlox()
-
-# Connect Regions through Adjacency Matrix
-blox = [GPe, STN]
-sys = [s.odesystem for s in blox]
-connect = [s.connector for s in blox]
-
-@parameters g_GPe_STN=1.0 g_STN_GPe=1.0
-adj_matrix = [1.0 g_STN_GPe;
-            g_STN_GPe*g_GPe_STN 1.0]
-
-#create equivalent graph
-g = MetaDiGraph()
-add_blox!(g,GPe)
-add_blox!(g,STN)
-add_edge!(g,1,1,:weight,1.0)
-add_edge!(g,1,2,:weight,g_STN_GPe)
-add_edge!(g,2,1,:weight,g_STN_GPe*g_GPe_STN)
-add_edge!(g,2,2,:weight,1.0)
-a = adjmatrixfromdigraph(g)
-
-#create equivalent graph with added Utilities Blox
-ggb = MetaDiGraph()
-add_blox!(ggb,GPe)
-add_blox!(ggb,STN)
-add_blox!(ggb,BP)
-add_edge!(ggb,1,1,:weight,1.0)
-add_edge!(ggb,1,2,:weight,g_STN_GPe)
-add_edge!(ggb,2,1,:weight,g_STN_GPe*g_GPe_STN)
-add_edge!(ggb,2,2,:weight,1.0)
-add_edge!(ggb,2,3,:weight,1.0)
-
-# I am not sure how to test that the two adjacency matrices are equal
-@test isequal(a,adj_matrix)
-
-@named two_regions = LinearConnections(sys=sys,adj_matrix=adj_matrix, connector=connect)
-@named two_regions_gr = ODEfromGraph(g)
-@named two_regions_grdb = ODEfromGraph(ggb)
-
-@test typeof(two_regions) == ODESystem
-@test typeof(two_regions_gr) == ODESystem
-@test typeof(two_regions_grdb) == ODESystem
-@test equations(two_regions) == equations(two_regions_gr)
-@test equations(two_regions_gr) == equations(two_regions_grdb)
-
-# parameter sharing test
-@parameters λ=400.0
-@named jr1 = JansenRitSCBlox(τ=0.014, H=20.0, λ=λ, r=0.1)
-@named jr2 = JansenRitSCBlox(τ=0.014, H=20.0, λ=λ, r=0.1)
-jg = MetaDiGraph()
-add_blox!(jg,jr1)
-add_blox!(jg,jr2)
-add_edge!(jg,1,2,:weight,1.0)
-add_edge!(jg,2,1,:weight,2.0)
-add_edge!(jg,1,1,:weight,-1.0)
-add_edge!(jg,2,2,:weight,-2.0)
-@named jg_graph = ODEfromGraph(jg)
-jg_graph_s = structural_simplify(jg_graph)
-@test length(parameters(jg_graph_s))==7
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 2bcdc82f..59c1ec1f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -8,15 +8,8 @@ if GROUP == "All" || GROUP == "Basics"
     @time @safetestset "Neurograph Tests" begin include("graphs.jl") end
 end
 
-# Commenting out because we no longer use ODEfromGraph for neural masses
-# @time @safetestset "ODE from Graph Tests" begin include("ode_from_graph.jl") end
-
-# @time @safetestset "Neural Signal Measurement Models Tests" begin include("measurementmodels.jl") end
-# @time @safetestset "Spectral Utilities Tests" begin include("spectraltools.jl") end
-
 if GROUP == "All" || GROUP == "Advanced"
     @time @safetestset "Data Fitting Tests" begin include("datafitting.jl") end 
-    #@time @safetestset "ODE from Graph and simulate" begin include("graph_to_dataframe.jl") end #This no longer works because it's built on the old blocks
     @time @safetestset "Learning Tests" begin include("learning.jl") end
     @time @safetestset "Control Tests" begin include("controllers.jl") end
     @time @safetestset "Source Tests" begin include("source_components.jl") end
@@ -24,16 +17,6 @@ if GROUP == "All" || GROUP == "Advanced"
     @time @safetestset "Cort-Cort plasticity Tests" begin include("plasticity.jl") end
 end
 
-if GROUP == "All" || GROUP == "GraphDynamics"
-    @time @safetestset "GraphDynamics vs MTK tests" begin include("GraphDynamicsTests/runtests.jl") end
-end
-
-# fitting tests should be at the end since they take the longest
-# removing them for now until we have real fitting tests
-# @time @safetestset "Fitting Tests" begin include("fitting.jl") end
-# @time @safetestset "Bayesian Fitting Tests" begin include("bayesian_fitting.jl") end
 
-# Extra tests illustrating use - could be removed once examples are written
-#@time @safetestset "New Jansen-Rit Tests" begin include("jansen_rit_component_tests.jl") end
-#@time @safetestset "New Larter-Breakspear Tests" begin include("new_LB_blox.jl") end
+@time @safetestset "GraphDynamics vs MTK tests" begin include("GraphDynamicsTests/runtests.jl") end
 
diff --git a/test/spectraltools.jl b/test/spectraltools.jl
deleted file mode 100644
index a1c6a59e..00000000
--- a/test/spectraltools.jl
+++ /dev/null
@@ -1,105 +0,0 @@
-using Neuroblox, Test
-
-import LinearAlgebra as la
-using Distributions
-using OrdinaryDiffEq, Plots
-using Statistics
-using MAT
-using DSP
-
-
-nd = 2  # number of dimensions
-
-p = 2  # number of time lags of MAR model
-f = 2.0.^(range(0,stop=5)) # frequencies at which to evaluate CSD
-dt = 1/(2*f[end]) # time step, inverse of sampling frequency
-dist = InverseWishart(nd*2, la.Matrix(1.0la.I, nd, nd))
-Σ = rand(dist)   # noise covariance matrix of MAR model
-a = [randn(nd, nd) for i = 1:p]   # MAR model parameters
-mar = Dict([("A", a), ("Σ", Σ), ("p", p)])
-
-csd = mar2csd(mar, f)
-a_est, Σ_est = csd2mar(csd, f, dt, p)
-
-@test_broken a ≈ a_est
-@test_broken Σ ≈ Σ_est
-
-
-"""
-Design Test Case (PowerSpectrum). 
-A circuit model is used, with an explicit paramter setting the oscillation frequency range of the system in radians/second:
-    ω = 4*2*π
-The test will count the number of peaks in a given time window and match that to the frequency containing the most power.
-"""
-freq_of_interest = 4
-# Create Circuit 
-@named STN = harmonic_oscillator(ω=freq_of_interest*2*π, ζ=1, k=(freq_of_interest*2*π)^2, h=5.0)
-# Run Simulation
-
-sys = [STN.odesystem]
-adj_mat = [1.0]
-connect = [STN.connector]
-@named STN_Circuit = LinearConnections(sys=sys, adj_matrix=adj_mat, connector=connect)
-
-sim_dur = 5.0                                                                         # Simulation time (seconds)
-prob = ODEProblem(structural_simplify(STN_Circuit), [], (0.0, sim_dur), [])
-sol = solve(prob, Tsit5(), dt = 0.001)
-
-"""
-powerspectrum test.
-The power spectrum of the solution should have a center frequency around that set by ω in Hz, with some tolerance.
-To Do: Employ a peak detection algorithm, useful for this test and signal processing. Number of peaks should match center frequency.
-"""
-f, pxx = Neuroblox.powerspectrum(sol[1,:], sim_dur, 1000, "auto", hanning)
-# 1) Make sure you can plot the data
-p1 = Plots.plot(sol.t, sol[1,:])
-p2 = Plots.plot(f, pxx, xlims=(0,10))
-# 2) Find the local maxima
-T = sim_dur
-df = 1/T
-find_max = findmax(pxx[1:length(f)-Int(ceil(freq_of_interest/df))])
-index_of_maximum = find_max[2]
-tol = 0.5
-@test index_of_maximum*df>4-tol
-@test index_of_maximum*df<4+tol
-
-"""
-bandpassfilter test
-Compute power spectrum before and after filtering the data.
-"""
-data = matread(joinpath(@__DIR__, "lfp_test_data.mat"))
-data = data["lfp"]
-f, pxx = Neuroblox.powerspectrum(data, length(data), 1000, "periodogram", hanning)
-
-lb = 12
-ub = 30
-signal = Neuroblox.bandpassfilter(data=data, lb=lb, ub=ub, fs=1000, order=4)
-f_signal, pxx_signal = Neuroblox.powerspectrum(signal, length(data), 1000, "periodogram", hanning)
-
-@test pxx_signal[1:lb] < pxx[1:lb]
-@test pxx_signal[ub:100] < pxx[ub:100]
-
-
-"""
-phaseangle test
-the imaginary part of the transformed data must sum to approximately zero
-"""
-data = matread(joinpath(@__DIR__, "lfp_test_data_phase.mat"))
-data = data["lfp_data"]
-phase = Neuroblox.phaseangle(data=data)
-
-@test minimum(phase) > -pi
-@test maximum(phase) < pi
-
-"""
-complexwavelet test
-Wavelets must have values near zero at both ends, as well as a mean value of zero
-"""
-data = matread(joinpath(@__DIR__, "lfp_test_data.mat"))
-wavelets = Neuroblox.complexwavelet(data["lfp"], 0.001, 2, 60)
-tol = 0.2
-@test real(wavelets[1][1]) < tol
-all_wavelets = Statistics.mean(real(wavelets))
-average_over_all = sum(all_wavelets)/length(all_wavelets)
-tol = 0.001
-@test sum(average_over_all) < 0 + tol
\ No newline at end of file