diff --git a/examples/jansen_rit_liu_et_al.jl b/examples/jansen_rit_liu_et_al.jl index 555a7199..dcfe5be0 100644 --- a/examples/jansen_rit_liu_et_al.jl +++ b/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) diff --git a/src/Neuroblox.jl b/src/Neuroblox.jl index 5cddbcc5..05623a89 100644 --- a/src/Neuroblox.jl +++ b/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! diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index ca83bd9f..c88756db 100644 --- a/src/blox/neural_mass.jl +++ b/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, η) diff --git a/test/components.jl b/test/components.jl index ec4087ca..9a91972c 100644 --- a/test/components.jl +++ b/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 diff --git a/test/jansen_rit_component_tests.jl b/test/jansen_rit_component_tests.jl index a809737a..f06309ac 100644 --- a/test/jansen_rit_component_tests.jl +++ b/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)