Getting Started with Neuroblox

This tutorial will introduce you to simulating brain dynamics using Neuroblox.

Example 1 : Building an oscillating circuit from two Wilson-Cowan Neural Mass Models

The Wilson–Cowan model describes the dynamics of interactions between populations of excitatory and inhibitory neurons. Each Wilson-Cowan Blox is described by the follwoing equations:

\[\frac{dE}{dt} = \frac{-E}{\tau_E} + \frac{1}{1 + \text{exp}(-a_E*(c_{EE}*E - c_{IE}*I - \theta_E + \eta*(\sum{jcn}))}\\[10pt] \frac{dI}{dt} = \frac{-I}{\tau_I} + \frac{1}{1 + exp(-a_I*(c_{EI}*E - c_{II}*I - \theta_I)}\]

Our first example is to simply combine two Wilson-Cowan Blox to build an oscillatory circuit

using Neuroblox
using DifferentialEquations
using Graphs
using MetaGraphs
using Plots

@named WC1 = WilsonCowan()
@named WC2 = WilsonCowan()

g = MetaDiGraph()
add_blox!.(Ref(g), [WC1, WC2])

adj = [-1 6; 6 -1]
create_adjacency_edges!(g, adj)

First, we create the two Wilson-Cowan Blox: WC1 and WC2. Next, we add the two Blox into a directed graph as nodes and then we are creating weighted edges between the two nodes using an adjacency matrix.

Now we are ready to build the ModelingToolkit System. Structural simplify creates the final set of equations in which all substiutions are made.

@named sys = system_from_graph(g)
sys = structural_simplify(sys)

\[ \begin{align} \frac{\mathrm{d} WC1_{+}E\left( t \right)}{\mathrm{d}t} =& \frac{ - WC1_{+}E\left( t \right)}{WC1_+\tau_{E}} + \frac{1}{1 + e^{ - WC1_{+}a_{E} \left( - WC1_+\theta_{E} + WC1_{+}c_{EE} WC1_{+}E\left( t \right) - WC1_{+}c_{IE} WC1_{+}I\left( t \right) + WC1_+\eta WC1_{+}jcn\left( t \right) \right)}} \\ \frac{\mathrm{d} WC1_{+}I\left( t \right)}{\mathrm{d}t} =& \frac{1}{1 + e^{ - WC1_{+}a_{I} \left( - WC1_+\theta_{I} + WC1_{+}c_{EI} WC1_{+}E\left( t \right) - WC1_{+}c_{II} WC1_{+}I\left( t \right) \right)}} + \frac{ - WC1_{+}I\left( t \right)}{WC1_+\tau_{I}} \\ \frac{\mathrm{d} WC2_{+}E\left( t \right)}{\mathrm{d}t} =& \frac{ - WC2_{+}E\left( t \right)}{WC2_+\tau_{E}} + \frac{1}{1 + e^{ - WC2_{+}a_{E} \left( - WC2_+\theta_{E} + WC2_{+}c_{EE} WC2_{+}E\left( t \right) - WC2_{+}c_{IE} WC2_{+}I\left( t \right) + WC2_+\eta WC2_{+}jcn\left( t \right) \right)}} \\ \frac{\mathrm{d} WC2_{+}I\left( t \right)}{\mathrm{d}t} =& \frac{ - WC2_{+}I\left( t \right)}{WC2_+\tau_{I}} + \frac{1}{1 + e^{ - WC2_{+}a_{I} \left( - WC2_+\theta_{I} + WC2_{+}c_{EI} WC2_{+}E\left( t \right) - WC2_{+}c_{II} WC2_{+}I\left( t \right) \right)}} \end{align} \]

To solve the system, we first create an ODEProblem and then solve it over the tspan of (0,100) using a stiff solver. The solution is saved every 0.1ms. The unit of time in Neuroblox is 1ms.

prob = ODEProblem(sys, [], (0.0, 100), [])
sol = solve(prob, Rodas4(), saveat=0.1)
plot(sol)

Example 2 : Building a Brain Circuit from literature using Neural Mass Models

In this example, we will construct a Parkinsons model from eight Jansen-Rit Neural Mass Models as described in Liu et al. (2020). DOI: 10.1016/j.neunet.2019.12.021. The Jansen-Rit Neural Mass model is defined by the following differential equations:

\[\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]\]

using Neuroblox
using DifferentialEquations
using Graphs
using MetaGraphs
using Plots

The original paper units are in seconds we therefore need to multiply all parameters with a common factor

τ_factor = 1000

@named Str = JansenRit(τ=0.0022*τ_factor, H=20/τ_factor, λ=300, r=0.3)
@named GPE = JansenRit(τ=0.04*τ_factor, cortical=false) # all default subcortical except τ
@named STN = JansenRit(τ=0.01*τ_factor, H=20/τ_factor, λ=500, r=0.1)
@named GPI = JansenRit(cortical=false) # default parameters subcortical Jansen Rit blox
@named Th  = JansenRit(τ=0.002*τ_factor, H=10/τ_factor, λ=20, r=5)
@named EI  = JansenRit(τ=0.01*τ_factor, H=20/τ_factor, λ=5, r=5)
@named PY  = JansenRit(cortical=true) # default parameters cortical Jansen Rit blox
@named II  = JansenRit(τ=2.0*τ_factor, H=60/τ_factor, λ=5, r=5)
blox = [Str, GPE, STN, GPI, Th, EI, PY, II]

8-element Vector{JansenRit}:
 JansenRit(Any[τ, H, λ, r], x(t), jcn(t), ODESystem(0x000000000000000d, Equation[Differential(t)(x(t)) ~ y(t) + (-2x(t)) / τ, Differential(t)(y(t)) ~ (-x(t)) / (τ^2) + (H*((2λ) / (1 + exp(-r*jcn(t))) - λ)) / τ], t, SymbolicUtils.BasicSymbolic{Real}[x(t), y(t), jcn(t)], SymbolicUtils.BasicSymbolic{Real}[τ, H, λ, r], nothing, Dict{Any, Any}(:H => H, :y => y(t), :λ => λ, :τ => τ, :jcn => jcn(t), :r => r, :x => x(t)), Any[], Equation[], Base.RefValue{Vector{Num}}(Num[]), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), :Str, ODESystem[], Dict{Any, Any}(τ => 2.2, H => 0.02, x(t) => 1.0, λ => 300, r => 0.3, jcn(t) => 0.0, y(t) => 1.0), Dict{Any, Any}(), nothing, nothing, Equation[], nothing, nothing, nothing, ModelingToolkit.SymbolicContinuousCallback[], ModelingToolkit.SymbolicDiscreteCallback[], nothing, nothing, nothing, nothing, nothing, false, nothing, nothing, nothing, nothing, nothing), nothing)
 JansenRit(Any[τ, H, λ, r], x(t), jcn(t), ODESystem(0x000000000000000e, Equation[Differential(t)(x(t)) ~ y(t) + (-2x(t)) / τ, Differential(t)(y(t)) ~ (-x(t)) / (τ^2) + (H*((2λ) / (1 + exp(-r*jcn(t))) - λ)) / τ], t, SymbolicUtils.BasicSymbolic{Real}[x(t), y(t), jcn(t)], SymbolicUtils.BasicSymbolic{Real}[τ, H, λ, r], nothing, Dict{Any, Any}(:H => H, :y => y(t), :λ => λ, :τ => τ, :jcn => jcn(t), :r => r, :x => x(t)), Any[], Equation[], Base.RefValue{Vector{Num}}(Num[]), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), :GPE, ODESystem[], Dict{Any, Any}(τ => 40.0, H => 0.02, x(t) => 1.0, λ => 400.0, r => 0.1, jcn(t) => 0.0, y(t) => 1.0), Dict{Any, Any}(), nothing, nothing, Equation[], nothing, nothing, nothing, ModelingToolkit.SymbolicContinuousCallback[], ModelingToolkit.SymbolicDiscreteCallback[], nothing, nothing, nothing, nothing, nothing, false, nothing, nothing, nothing, nothing, nothing), nothing)
 JansenRit(Any[τ, H, λ, r], x(t), jcn(t), ODESystem(0x000000000000000f, Equation[Differential(t)(x(t)) ~ y(t) + (-2x(t)) / τ, Differential(t)(y(t)) ~ (-x(t)) / (τ^2) + (H*((2λ) / (1 + exp(-r*jcn(t))) - λ)) / τ], t, SymbolicUtils.BasicSymbolic{Real}[x(t), y(t), jcn(t)], SymbolicUtils.BasicSymbolic{Real}[τ, H, λ, r], nothing, Dict{Any, Any}(:H => H, :y => y(t), :λ => λ, :τ => τ, :jcn => jcn(t), :r => r, :x => x(t)), Any[], Equation[], Base.RefValue{Vector{Num}}(Num[]), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), :STN, ODESystem[], Dict{Any, Any}(τ => 10.0, H => 0.02, x(t) => 1.0, λ => 500, r => 0.1, jcn(t) => 0.0, y(t) => 1.0), Dict{Any, Any}(), nothing, nothing, Equation[], nothing, nothing, nothing, ModelingToolkit.SymbolicContinuousCallback[], ModelingToolkit.SymbolicDiscreteCallback[], nothing, nothing, nothing, nothing, nothing, false, nothing, nothing, nothing, nothing, nothing), nothing)
 JansenRit(Any[τ, H, λ, r], x(t), jcn(t), ODESystem(0x0000000000000010, Equation[Differential(t)(x(t)) ~ y(t) + (-2x(t)) / τ, Differential(t)(y(t)) ~ (-x(t)) / (τ^2) + (H*((2λ) / (1 + exp(-r*jcn(t))) - λ)) / τ], t, SymbolicUtils.BasicSymbolic{Real}[x(t), y(t), jcn(t)], SymbolicUtils.BasicSymbolic{Real}[τ, H, λ, r], nothing, Dict{Any, Any}(:H => H, :y => y(t), :λ => λ, :τ => τ, :jcn => jcn(t), :r => r, :x => x(t)), Any[], Equation[], Base.RefValue{Vector{Num}}(Num[]), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), :GPI, ODESystem[], Dict{Any, Any}(τ => 14, H => 0.02, x(t) => 1.0, λ => 400.0, r => 0.1, jcn(t) => 0.0, y(t) => 1.0), Dict{Any, Any}(), nothing, nothing, Equation[], nothing, nothing, nothing, ModelingToolkit.SymbolicContinuousCallback[], ModelingToolkit.SymbolicDiscreteCallback[], nothing, nothing, nothing, nothing, nothing, false, nothing, nothing, nothing, nothing, nothing), nothing)
 JansenRit(Any[τ, H, λ, r], x(t), jcn(t), ODESystem(0x0000000000000011, Equation[Differential(t)(x(t)) ~ y(t) + (-2x(t)) / τ, Differential(t)(y(t)) ~ (-x(t)) / (τ^2) + (H*((2λ) / (1 + exp(-r*jcn(t))) - λ)) / τ], t, SymbolicUtils.BasicSymbolic{Real}[x(t), y(t), jcn(t)], SymbolicUtils.BasicSymbolic{Real}[τ, H, λ, r], nothing, Dict{Any, Any}(:H => H, :y => y(t), :λ => λ, :τ => τ, :jcn => jcn(t), :r => r, :x => x(t)), Any[], Equation[], Base.RefValue{Vector{Num}}(Num[]), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), :Th, ODESystem[], Dict{Any, Any}(τ => 2.0, H => 0.01, x(t) => 1.0, λ => 20, r => 5, jcn(t) => 0.0, y(t) => 1.0), Dict{Any, Any}(), nothing, nothing, Equation[], nothing, nothing, nothing, ModelingToolkit.SymbolicContinuousCallback[], ModelingToolkit.SymbolicDiscreteCallback[], nothing, nothing, nothing, nothing, nothing, false, nothing, nothing, nothing, nothing, nothing), nothing)
 JansenRit(Any[τ, H, λ, r], x(t), jcn(t), ODESystem(0x0000000000000012, Equation[Differential(t)(x(t)) ~ y(t) + (-2x(t)) / τ, Differential(t)(y(t)) ~ (-x(t)) / (τ^2) + (H*((2λ) / (1 + exp(-r*jcn(t))) - λ)) / τ], t, SymbolicUtils.BasicSymbolic{Real}[x(t), y(t), jcn(t)], SymbolicUtils.BasicSymbolic{Real}[τ, H, λ, r], nothing, Dict{Any, Any}(:H => H, :y => y(t), :λ => λ, :τ => τ, :jcn => jcn(t), :r => r, :x => x(t)), Any[], Equation[], Base.RefValue{Vector{Num}}(Num[]), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), :EI, ODESystem[], Dict{Any, Any}(τ => 10.0, H => 0.02, x(t) => 1.0, λ => 5, r => 5, jcn(t) => 0.0, y(t) => 1.0), Dict{Any, Any}(), nothing, nothing, Equation[], nothing, nothing, nothing, ModelingToolkit.SymbolicContinuousCallback[], ModelingToolkit.SymbolicDiscreteCallback[], nothing, nothing, nothing, nothing, nothing, false, nothing, nothing, nothing, nothing, nothing), nothing)
 JansenRit(Any[τ, H, λ, r], x(t), jcn(t), ODESystem(0x0000000000000013, Equation[Differential(t)(x(t)) ~ y(t) + (-2x(t)) / τ, Differential(t)(y(t)) ~ (-x(t)) / (τ^2) + (H*((2λ) / (1 + exp(-r*jcn(t))) - λ)) / τ], t, SymbolicUtils.BasicSymbolic{Real}[x(t), y(t), jcn(t)], SymbolicUtils.BasicSymbolic{Real}[τ, H, λ, r], nothing, Dict{Any, Any}(:H => H, :y => y(t), :λ => λ, :τ => τ, :jcn => jcn(t), :r => r, :x => x(t)), Any[], Equation[], Base.RefValue{Vector{Num}}(Num[]), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), :PY, ODESystem[], Dict{Any, Any}(τ => 1, H => 0.02, x(t) => 1.0, λ => 5.0, r => 0.15, jcn(t) => 0.0, y(t) => 1.0), Dict{Any, Any}(), nothing, nothing, Equation[], nothing, nothing, nothing, ModelingToolkit.SymbolicContinuousCallback[], ModelingToolkit.SymbolicDiscreteCallback[], nothing, nothing, nothing, nothing, nothing, false, nothing, nothing, nothing, nothing, nothing), nothing)
 JansenRit(Any[τ, H, λ, r], x(t), jcn(t), ODESystem(0x0000000000000014, Equation[Differential(t)(x(t)) ~ y(t) + (-2x(t)) / τ, Differential(t)(y(t)) ~ (-x(t)) / (τ^2) + (H*((2λ) / (1 + exp(-r*jcn(t))) - λ)) / τ], t, SymbolicUtils.BasicSymbolic{Real}[x(t), y(t), jcn(t)], SymbolicUtils.BasicSymbolic{Real}[τ, H, λ, r], nothing, Dict{Any, Any}(:H => H, :y => y(t), :λ => λ, :τ => τ, :jcn => jcn(t), :r => r, :x => x(t)), Any[], Equation[], Base.RefValue{Vector{Num}}(Num[]), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), :II, ODESystem[], Dict{Any, Any}(τ => 2000.0, H => 0.06, x(t) => 1.0, λ => 5, r => 5, jcn(t) => 0.0, y(t) => 1.0), Dict{Any, Any}(), nothing, nothing, Equation[], nothing, nothing, nothing, ModelingToolkit.SymbolicContinuousCallback[], ModelingToolkit.SymbolicDiscreteCallback[], nothing, nothing, nothing, nothing, nothing, false, nothing, nothing, nothing, nothing, nothing), nothing)

Again, we create a graph and add the Blox as nodes

g = MetaDiGraph()
add_blox!.(Ref(g), blox)

params = @parameters C_Cor=60 C_BG_Th=60 C_Cor_BG_Th=5 C_BG_Th_Cor=5

\[ \begin{equation} \left[ \begin{array}{c} C_{Cor} \\ C_{BG\_Th} \\ C_{Cor\_BG\_Th} \\ C_{BG\_Th\_Cor} \\ \end{array} \right] \end{equation} \]

ModelingToolkit allows us to create parameters that can be passed into the equations symbolically.

We add edges as specified in Table 2 of Liu et al. We only implemented a subset of the nodes and edges to describe a less complex version of the model. Edges can also be created using an adjacency matrix as in the previous example.

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,1,1,:weight, -0.5*C_BG_Th)
add_edge!(g,1,2,:weight, C_BG_Th)
add_edge!(g,2,1,:weight, -0.5*C_BG_Th)
add_edge!(g,2,5,:weight, C_Cor_BG_Th)
add_edge!(g,3,1,:weight, -0.5*C_BG_Th)
add_edge!(g,3,2,:weight, C_BG_Th)
add_edge!(g,4,3,:weight, -0.5*C_BG_Th)
add_edge!(g,4,4,:weight, C_BG_Th_Cor)

true

Now we are ready to build the ModelingToolkit System and apply structural simplification to the equations.

@named final_system = system_from_graph(g)
final_system_sys = structural_simplify(final_system)

\[ \begin{align} \frac{\mathrm{d} Str_{+}x\left( t \right)}{\mathrm{d}t} =& \frac{ - 2 Str_{+}x\left( t \right)}{Str_+\tau} + Str_{+}y\left( t \right) \\ \frac{\mathrm{d} Str_{+}y\left( t \right)}{\mathrm{d}t} =& \frac{ - Str_{+}x\left( t \right)}{Str_+\tau^{2}} + \frac{Str_{+}H \left( - Str_+\lambda + \frac{2 Str_+\lambda}{1 + e^{ - Str_{+}r Str_{+}jcn\left( t \right)}} \right)}{Str_+\tau} \\ \frac{\mathrm{d} GPE_{+}x\left( t \right)}{\mathrm{d}t} =& GPE_{+}y\left( t \right) + \frac{ - 2 GPE_{+}x\left( t \right)}{GPE_+\tau} \\ \frac{\mathrm{d} GPE_{+}y\left( t \right)}{\mathrm{d}t} =& \frac{GPE_{+}H \left( - GPE_+\lambda + \frac{2 GPE_+\lambda}{1 + e^{ - GPE_{+}r GPE_{+}jcn\left( t \right)}} \right)}{GPE_+\tau} + \frac{ - GPE_{+}x\left( t \right)}{GPE_+\tau^{2}} \\ \frac{\mathrm{d} STN_{+}x\left( t \right)}{\mathrm{d}t} =& \frac{ - 2 STN_{+}x\left( t \right)}{STN_+\tau} + STN_{+}y\left( t \right) \\ \frac{\mathrm{d} STN_{+}y\left( t \right)}{\mathrm{d}t} =& \frac{STN_{+}H \left( - STN_+\lambda + \frac{2 STN_+\lambda}{1 + e^{ - STN_{+}r STN_{+}jcn\left( t \right)}} \right)}{STN_+\tau} + \frac{ - STN_{+}x\left( t \right)}{STN_+\tau^{2}} \\ \frac{\mathrm{d} GPI_{+}x\left( t \right)}{\mathrm{d}t} =& \frac{ - 2 GPI_{+}x\left( t \right)}{GPI_+\tau} + GPI_{+}y\left( t \right) \\ \frac{\mathrm{d} GPI_{+}y\left( t \right)}{\mathrm{d}t} =& \frac{ - GPI_{+}x\left( t \right)}{GPI_+\tau^{2}} + \frac{GPI_{+}H \left( - GPI_+\lambda + \frac{2 GPI_+\lambda}{1 + e^{ - GPI_{+}r GPI_{+}jcn\left( t \right)}} \right)}{GPI_+\tau} \\ \frac{\mathrm{d} Th_{+}x\left( t \right)}{\mathrm{d}t} =& \frac{ - 2 Th_{+}x\left( t \right)}{Th_+\tau} + Th_{+}y\left( t \right) \\ \frac{\mathrm{d} Th_{+}y\left( t \right)}{\mathrm{d}t} =& \frac{Th_{+}H \left( - Th_+\lambda + \frac{2 Th_+\lambda}{1 + e^{ - Th_{+}r Th_{+}jcn\left( t \right)}} \right)}{Th_+\tau} + \frac{ - Th_{+}x\left( t \right)}{Th_+\tau^{2}} \\ \frac{\mathrm{d} EI_{+}x\left( t \right)}{\mathrm{d}t} =& \frac{ - 2 EI_{+}x\left( t \right)}{EI_+\tau} + EI_{+}y\left( t \right) \\ \frac{\mathrm{d} EI_{+}y\left( t \right)}{\mathrm{d}t} =& \frac{ - EI_{+}x\left( t \right)}{EI_+\tau^{2}} + \frac{EI_{+}H \left( - EI_+\lambda + \frac{2 EI_+\lambda}{1 + e^{ - EI_{+}r EI_{+}jcn\left( t \right)}} \right)}{EI_+\tau} \\ \frac{\mathrm{d} PY_{+}x\left( t \right)}{\mathrm{d}t} =& \frac{ - 2 PY_{+}x\left( t \right)}{PY_+\tau} + PY_{+}y\left( t \right) \\ \frac{\mathrm{d} PY_{+}y\left( t \right)}{\mathrm{d}t} =& \frac{PY_{+}H \left( - PY_+\lambda + \frac{2 PY_+\lambda}{1 + e^{ - PY_{+}r PY_{+}jcn\left( t \right)}} \right)}{PY_+\tau} + \frac{ - PY_{+}x\left( t \right)}{PY_+\tau^{2}} \\ \frac{\mathrm{d} II_{+}x\left( t \right)}{\mathrm{d}t} =& \frac{ - 2 II_{+}x\left( t \right)}{II_+\tau} + II_{+}y\left( t \right) \\ \frac{\mathrm{d} II_{+}y\left( t \right)}{\mathrm{d}t} =& \frac{ - II_{+}x\left( t \right)}{II_+\tau^{2}} + \frac{II_{+}H \left( - II_+\lambda + \frac{2 II_+\lambda}{1 + e^{ - II_{+}r II_{+}jcn\left( t \right)}} \right)}{II_+\tau} \end{align} \]

Our Jansen-Rit model allows delayed edges, and we therefore need to collect those delays (in our case all delays are zero). Then we build a Delayed Differential Equations Problem (DDEProblem).

sim_dur = 1000.0 # Simulate for 1 second
prob = ODEProblem(final_system_sys,
    [],
    (0.0, sim_dur))

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1000.0)
u0: 16-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0

We select an algorithm and solve the system

alg = Tsit5()
sol_dde_no_delays = solve(prob, alg, saveat=1)
plot(sol_dde_no_delays)

In a later tutorial, we will show how to introduce edge delays.