GitHub

Getting Started

Getting Started with Julia

Here we would like to summarize some resources for people that are interested in learning more about the Julia language before or while exploring Neuroblox. Please follow the links below for introductory material on the language that is inclusive to all users; people familiar with programming or not, people with a mathematics, engineering, or science background :

Getting Started with Neuroblox

This example will introduce you to simulating brain dynamics using Neuroblox. We will create a simple oscillating circuit using two Wilson-Cowan neural mass models [1]. The Wilson-Cowan model is one of the most influential models in computational neuroscience [2], describing the dynamics of interactions between populations of excitatory and inhibitory neurons.

The Wilson-Cowan Model

Each Wilson-Cowan neural mass is described by the following equations:

\[\begin{align} \nonumber \frac{dE}{dt} &= \frac{-E}{\tau_E} + S_E(c_{EE}E - c_{IE}I + \eta\textstyle\sum{jcn})\\[10pt] \nonumber \frac{dI}{dt} &= \frac{-I}{\tau_I} + S_I\left(c_{EI}E - c_{II}I\right) \end{align}\]

where $E$ and $I$ denote the activity levels of the excitatory and inhibitory populations, respectively. The terms $\frac{dE}{dt}$ and $\frac{dI}{dt}$ describe the rate of change of these activity levels over time. The parameters $\tau_E$ and $\tau_I$ are time constants analogous to membrane time constants in single neuron models, determining how quickly the excitatory and inhibitory populations respond to changes. The coefficients $c_{EE}$ and $c_{II}$ represent self-interaction (or feedback) within excitatory and inhibitory populations, while $c_{IE}$ and $c_{EI}$ represent the cross-interactions between the two populations. The term $\eta\sum{jcn}$ represents external input to the excitatory population from other brain regions or external stimuli, with $\eta$ acting as a scaling factor. While $S_E$ and $S_I$ are sigmoid functions that represent the responses of neuronal populations to input stimuli, defined as:

\[S_k(x) = \frac{1}{1 + exp(-a_kx - \theta_k)}\]

where $a_k$ and $\theta_k$ determine the steepness and threshold of the response, respectively.

Building the Circuit

Let's create an oscillating circuit by connecting two Wilson-Cowan neural masses:

using Neuroblox
using OrdinaryDiffEq
using CairoMakie

# Create two Wilson-Cowan blox
@named WC1 = WilsonCowan()
@named WC2 = WilsonCowan()

# Create a graph to represent our circuit
g = MetaDiGraph()
add_blox!.(Ref(g), [WC1, WC2])

# Define the connectivity between the neural masses
add_edge!(g, WC1 => WC1; weight = -1) ## recurrent connection from WC1 to itself
add_edge!(g, WC1 => WC2; weight = 7) ## connection from WC1 to WC2
add_edge!(g, WC2 => WC1; weight = 4) ## connection from WC2 to WC1
add_edge!(g, WC2 => WC2; weight = -1) ## recurrent connection from WC2 to itself
true

Here, we've created two Wilson-Cowan blox and connected them as nodes in a directed graph. Each blox connects to itself and to the other blox.

By default, the output of each Wilson-Cowan blox is its excitatory activity (E). The negative self-connections (-1) provide inhibitory feedback, while the positive inter-blox connections (6) provide strong excitatory coupling. This setup creates an oscillatory dynamic between the two Wilson-Cowan units.

Creating the Model

Now, let's build the complete model:

@named sys = system_from_graph(g)

\[ \begin{align} \frac{\mathrm{d} \mathtt{WC1.E}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{WC1.E}\left( t \right)}{\mathtt{WC1.\tau\_E}} + \frac{1}{1 + e^{ - \mathtt{WC1.a\_E} \left( - \mathtt{WC1.\theta\_E} + \mathtt{WC1.c\_EE} \mathtt{WC1.E}\left( t \right) - \mathtt{WC1.c\_IE} \mathtt{WC1.I}\left( t \right) + \mathtt{WC1.\eta} \mathtt{WC1.jcn}\left( t \right) \right)}} \\ \frac{\mathrm{d} \mathtt{WC1.I}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{WC1.I}\left( t \right)}{\mathtt{WC1.\tau\_I}} + \frac{1}{1 + e^{\mathtt{WC1.a\_I} \left( \mathtt{WC1.\theta\_I} - \mathtt{WC1.c\_EI} \mathtt{WC1.E}\left( t \right) + \mathtt{WC1.c\_II} \mathtt{WC1.I}\left( t \right) \right)}} \\ \frac{\mathrm{d} \mathtt{WC2.E}\left( t \right)}{\mathrm{d}t} &= \frac{1}{1 + e^{\mathtt{WC2.a\_E} \left( \mathtt{WC2.\theta\_E} - \mathtt{WC2.c\_EE} \mathtt{WC2.E}\left( t \right) + \mathtt{WC2.c\_IE} \mathtt{WC2.I}\left( t \right) - \mathtt{WC2.\eta} \mathtt{WC2.jcn}\left( t \right) \right)}} + \frac{ - \mathtt{WC2.E}\left( t \right)}{\mathtt{WC2.\tau\_E}} \\ \frac{\mathrm{d} \mathtt{WC2.I}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{WC2.I}\left( t \right)}{\mathtt{WC2.\tau\_I}} + \frac{1}{1 + e^{\mathtt{WC2.a\_I} \left( \mathtt{WC2.\theta\_I} - \mathtt{WC2.c\_EI} \mathtt{WC2.E}\left( t \right) + \mathtt{WC2.c\_II} \mathtt{WC2.I}\left( t \right) \right)}} \end{align} \]

This creates a differential equations system from our graph representation using ModelingToolkit and symbolically simplifies it for efficient computation.

Simulating the Model

We are now ready to simulate our model. The following code creates and solves an ODEProblem for our system, simulating 100 time units of activity. In Neuroblox, the default time unit is milliseconds. We use Rodas4, a solver efficient for stiff problems. The solution is saved every 0.1 ms, allowing us to observe the detailed evolution of the system's behavior.

prob = ODEProblem(sys, [], (0.0, 100), [])
sol = solve(prob, Rodas4(), saveat=0.1)
retcode: Success
Interpolation: 1st order linear
t: 1001-element Vector{Float64}:
   0.0
   0.1
   0.2
   0.3
   0.4
   0.5
   0.6
   0.7
   0.8
   0.9
   ⋮
  99.2
  99.3
  99.4
  99.5
  99.6
  99.7
  99.8
  99.9
 100.0
u: 1001-element Vector{Vector{Float64}}:
 [1.0, 1.0, 1.0, 1.0]
 [0.9493257635566694, 0.999997227549365, 0.9968352574060094, 0.9999983614072379]
 [0.8972980729622233, 0.9999897131895668, 0.9922642579108254, 0.9999967482029593]
 [0.8440391811260224, 0.9999681723259077, 0.9854388758920339, 0.999995076307697]
 [0.7898349351204718, 0.9999082303969463, 0.9750726533220467, 0.9999931987811771]
 [0.7351468545581609, 0.9997317040527284, 0.9591819783728336, 0.9999908111737857]
 [0.6805567621145077, 0.999147043910833, 0.9353051118230763, 0.9999870244686659]
 [0.6267161119059087, 0.9976995169055252, 0.9005839696994439, 0.9999806570867977]
 [0.5745445542719343, 0.9931773367150695, 0.8545965249945568, 0.9999637048536514]
 [0.5248014014381632, 0.9822251893286915, 0.7986897069764031, 0.9999118449047143]
 ⋮
 [0.28279518476405013, 0.2324195285296703, 0.31706361763697516, 0.3471838455135596]
 [0.278177945925614, 0.22328371239125838, 0.31279792715965477, 0.3332261043931271]
 [0.27448088654236413, 0.21430539986619204, 0.309809926269029, 0.31991542196228917]
 [0.27181878342800286, 0.20567351013756544, 0.30820728524811536, 0.30756838705022393]
 [0.27027557823071185, 0.19759133790123104, 0.30805675975403085, 0.29650516666840665]
 [0.2698099045024096, 0.19023716779344882, 0.30928326784170596, 0.28692620357578835]
 [0.27050297529084016, 0.1837315093688662, 0.31192346211382116, 0.2790718675463416]
 [0.27238765332993337, 0.17825467421576802, 0.3159222625765795, 0.27326609552613756]
 [0.27549680135361837, 0.17398697392244247, 0.3212245892361824, 0.26983282446125106]

Plotting simulation results

Finally, let us plot the E and I states of the first component, WC1. To do this we will use the state_timeseries function that extracts the timeseries of a blox state from the solution object.

E1 = state_timeseries(WC1, sol, "E")
I1 = state_timeseries(WC1, sol, "I")

fig = Figure()
ax = Axis(fig[1,1]; xlabel = "time (ms)")
lines!(ax, sol.t, E1, label = "E")
lines!(ax, sol.t, I1, label = "I")
Legend(fig[1,2], ax)
fig
Example block output

[1] Wilson, H. R., & Cowan, J. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophysical journal, 12(1), 1-24.

[2] Destexhe, A., & Sejnowski, T. J. (2009). The Wilson–Cowan model, 36 years later. Biological cybernetics, 101(1), 1-2.

This page was generated using Literate.jl.