Switching default Jansen-Rit to ODE

Current Jansen-Rit is renamed JansenRitDelayed, and JansenRit is now back to being an ODE

examples/jansen_rit_liu_et_al.jl

@@ -50,19 +50,27 @@ final_system_sys = structural_simplify(final_system)
 # Collect the graph delays and create a DDEProblem. This will all be zero in this case.
 final_delays = graph_delays(g)
 sim_dur = 1000.0 # Simulate for 1 second
-prob = DDEProblem(final_system_sys,
+prob = ODEProblem(final_system_sys,
     [],
-    (0.0, sim_dur),
-    constant_lags = final_delays)
+    (0.0, sim_dur))
 
 # Select the algorihm. MethodOfSteps will return the same as Vern7() in this case because there are no non-zero delays, but is required since this is a DDEProblem.
-alg = MethodOfSteps(Vern7())
+alg = Vern7()
 sol_dde_no_delays = solve(prob, alg, saveat=1)
 
 
 # Example of delayed connections
 
 # First, recreate the graph to remove previous connections
+@named Str = JansenRitDelayed(τ=0.0022*τ_factor, H=20, λ=300, r=0.3)
+@named GPE = JansenRitDelayed(τ=0.04*τ_factor, cortical=false) # all default subcortical except τ
+@named STN = JansenRitDelayed(τ=0.01*τ_factor, H=20, λ=500, r=0.1)
+@named GPI = JansenRitDelayed(cortical=false) # default parameters subcortical Jansen Rit blox
+@named Th  = JansenRitDelayed(τ=0.002*τ_factor, H=10, λ=20, r=5)
+@named EI  = JansenRitDelayed(τ=0.01*τ_factor, H=20, λ=5, r=5)
+@named PY  = JansenRitDelayed(cortical=true) # default parameters cortical Jansen Rit blox
+@named II  = JansenRitDelayed(τ=2.0*τ_factor, H=60, λ=5, r=5)
+blox = [Str, GPE, STN, GPI, Th, EI, PY, II]
 g = MetaDiGraph()
 add_blox!.(Ref(g), blox)
 

src/Neuroblox.jl

@@ -185,7 +185,7 @@ export JansenRitSPM12, next_generation, qif_neuron, if_neuron, hh_neuron_excitat
     hh_neuron_inhibitory, synaptic_network, van_der_pol
 export IFNeuronBlox, LIFNeuronBlox, QIFNeuronBlox, HHNeuronExciBlox, HHNeuronInhibBlox,
     CanonicalMicroCircuitBlox, WinnerTakeAllBlox, CorticalBlox, SuperCortical
-export LinearNeuralMass, HarmonicOscillator, JansenRit, WilsonCowan, LarterBreakspear, NextGenerationBlox, NextGenerationResolvedBlox, NextGenerationEIBlox
+export LinearNeuralMass, HarmonicOscillator, JansenRit, JansenRitDelayed, WilsonCowan, LarterBreakspear, NextGenerationBlox, NextGenerationResolvedBlox, NextGenerationEIBlox
 export Matrisome, Striosome, Striatum, GPi, GPe, Thalamus, STN, TAN, SNc
 export HebbianPlasticity, HebbianModulationPlasticity
 export Agent, ClassificationEnvironment, GreedyPolicy, reset!

src/blox/neural_mass.jl

@@ -156,7 +156,7 @@ end
 """
     JansenRit(name, namespace, τ, H, λ, r, cortical)
 
-    Create a Jansen Rit blox as described in Liu et al.
+    Create a Jansen Rit blox with optional delays as described in Liu et al.
     The formal definition of this blox is:
 
 ```math
@@ -179,13 +179,13 @@ Citations:
 1. Liu C, Zhou C, Wang J, Fietkiewicz C, Loparo KA. The role of coupling connections in a model of the cortico-basal ganglia-thalamocortical neural loop for the generation of beta oscillations. Neural Netw. 2020 Mar;123:381-392. doi: 10.1016/j.neunet.2019.12.021.
 
 """
-struct JansenRit <: NeuralMassBlox
+struct JansenRitDelayed <: NeuralMassBlox
     params
     output
     jcn
     odesystem
     namespace
-    function JansenRit(;name,
+    function JansenRitDelayed(;name,
                         namespace=nothing,
                         τ=nothing, 
                         H=nothing, 
@@ -211,6 +211,64 @@ struct JansenRit <: NeuralMassBlox
     end
 end
 
+"""
+    JansenRit(name, namespace, τ, H, λ, r, cortical)
+
+    Create a Jansen Rit blox as described in Liu et al.
+    The formal definition of this blox is:
+
+```math
+\\frac{dx}{dt} = y-\\frac{2}{\\tau}x
+\\frac{dy}{dt} = -\\frac{x}{\\tau^2} + \\frac{H}{\\tau} [\\frac{2\\lambda}{1+\\text{exp}(-r*\\sum{jcn})} - \\lambda]
+```
+
+where ``jcn`` is any input to the blox.
+
+Arguments:
+- name: Name given to ODESystem object within the blox.
+- namespace: Additional namespace above name if needed for inheritance.
+- τ: Time constant. This is changed from the original source as the time constant was in seconds, while all our blocks are in milliseconds.
+- H: See equation for use.
+- λ: See equation for use.
+- r: See equation for use.
+- cortical: Boolean to determine whether to use cortical or subcortical parameters. Specifying any of the parameters above will override this.
+
+Citations:
+1. Liu C, Zhou C, Wang J, Fietkiewicz C, Loparo KA. The role of coupling connections in a model of the cortico-basal ganglia-thalamocortical neural loop for the generation of beta oscillations. Neural Netw. 2020 Mar;123:381-392. doi: 10.1016/j.neunet.2019.12.021.
+
+"""
+struct JansenRit <: NeuralMassBlox
+    params
+    output
+    jcn
+    odesystem
+    namespace
+    function JansenRit(;name,
+                        namespace=nothing,
+                        τ=nothing, 
+                        H=nothing, 
+                        λ=nothing, 
+                        r=nothing, 
+                        cortical=true)
+
+        τ = isnothing(τ) ? (cortical ? 1 : 14) : τ
+        H = isnothing(H) ? 20.0 : H # H doesn't have different parameters for cortical and subcortical
+        λ = isnothing(λ) ? (cortical ? 5.0 : 400.0) : λ
+        r = isnothing(r) ? (cortical ? 0.15 : 0.1) : r
+
+        # p = progress_scope(τ, H, λ, r)
+        p = paramscoping(τ=τ, H=H, λ=λ, r=r)
+        τ, H, λ, r = p
+        sts = @variables x(t)=1.0 [output=true] y(t)=1.0 jcn(t)=0.0 [input=true] 
+        eqs = [D(x) ~ y - ((2/τ)*x),
+               D(y) ~ -x/(τ*τ) + (H/τ)*((2*λ)/(1 + exp(-r*(jcn))) - λ)]
+        sys = System(eqs, name=name)
+        #can't use outputs because x(t) is Num by then
+        #wrote inputs similarly to keep consistent
+        new(p, sts[1], sts[3], sys, namespace)
+    end
+end
+
 """
     WilsonCown(name, namespace, τ_E, τ_I, a_E, a_I, c_EE, c_IE, c_EI, c_II, θ_E, θ_I, η)
 

test/components.jl

@@ -90,11 +90,10 @@ New Jansen-Rit tests
     final_delays = graph_delays(g)
     sim_dur = 2000.0 # Simulate for 2 Seconds
     final_system_sys = structural_simplify(final_system)
-    prob = DDEProblem(final_system_sys,
+    prob = ODEProblem(final_system_sys,
         [],
-        (0.0, sim_dur),
-        constant_lags = final_delays)
-    alg = MethodOfSteps(Vern7())
+        (0.0, sim_dur))
+    alg = Vern7()
     sol_dde_no_delays = solve(prob, alg, saveat=1)
     @test sol_dde_no_delays.retcode == ReturnCode.Success
 end

test/jansen_rit_component_tests.jl

@@ -97,11 +97,10 @@ add_edge!(g, 8, 8, Dict(:weight => 3.3*C_Cor))
 final_delays = graph_delays(g)
 sim_dur = 2000.0 # Simulate for 2 Seconds
 final_system_sys = structural_simplify(final_system)
-prob = DDEProblem(final_system_sys,
+prob = ODEProblem(final_system_sys,
     [],
-    (0.0, sim_dur),
-    constant_lags = final_delays)
-alg = MethodOfSteps(Vern7())
+    (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)