Anandpathak31/issue5

docs/Project.toml

@@ -1,5 +1,6 @@
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"

docs/make.jl

@@ -3,6 +3,8 @@ using Documenter
 
 cp("./docs/Manifest.toml", "./docs/src/assets/Manifest.toml", force = true)
 cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true)
+#cp("./Manifest.toml", "./src/assets/Manifest.toml", force = true)
+#cp("./Project.toml", "./src/assets/Project.toml", force = true)
 
 DocMeta.setdocmeta!(Neuroblox, :DocTestSetup, :(using Neuroblox); recursive = true)
 

docs/pages.jl

@@ -3,7 +3,15 @@
 pages = ["index.md",
     "Getting Started" => "getting_started.md",
     "Tutorials" => Any[
-        "tutorials/resting_state_wb.md"
+        "tutorials/beginner_tutorial.md",
+        "tutorials/resting_state_wb.md",
+        #"tutorials/resting_state_wb_delays.md"
+        #"tutorials/DCM.md"
+        #"tutorials/Adam_basal_ganglia.md"
+        #"tutorial/winner_takes_all.md"
+        "tutorials/neural_assembly.md"
+        #"tutorial/reinforcement_learning.md"
+        #"tutorial/cortical.md"
     ],
     #"Manual" => Any[],
     "API" => "api.md",

docs/src/getting_started.md

@@ -2,16 +2,34 @@
 
 This tutorial will introduce you to simulating brain dynamics using Neuroblox.
 
-## Example 1 : Building an oscillating circuit from two Wilson-Cowan Neural Mass Models
+## Example 1: Creating a Simple Neural Circuit
 
-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:
+In this example, we'll 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:
+
+```math
+\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:
 
 ```math
-\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)}
+S_k(x) = \frac{1}{1 + exp(-a_kx - \theta_k)}
 ```
 
-Our first example is to simply combine two Wilson-Cowan Blox to build an oscillatory circuit
+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:
 
 ```@example Wilson-Cowan
 using Neuroblox
@@ -20,79 +38,128 @@ using Graphs
 using MetaGraphs
 using Plots
 
+# 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])
 
-adj = [-1 6; 6 -1]
+# Define the connectivity between the neural masses
+adj = [-1 7; 4 -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.
+Here, we've created two Wilson-Cowan Blox and connected them as nodes in a directed graph. The `adj` matrix defines the weighted edges between these nodes. Each entry `adj[i,j]` represents how the output of blox `j` influences the input of blox `i`:
+
+- Diagonal elements (`adj[1,1]` and `adj[2,2]`): Self-connections, adding feedback to each blox.
+- Off-diagonal elements (`adj[1,2]` and `adj[2,1]`): Inter-blox connections, determining how each blox influences the other.
 
-Now we are ready to build the ModelingToolkit System.  Structural simplify creates the final set of equations in which all substiutions are made.
+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:
 
 ```@example Wilson-Cowan
 @named sys = system_from_graph(g)
 sys = structural_simplify(sys)
 ```
 
-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.
+This creates a differential equations system from our graph representation using ModelingToolkit and symbolically simplifies it for efficient computation.
+
+### Simulating the Model
+
+Finally, let's 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.
 
 ```@example Wilson-Cowan
 prob = ODEProblem(sys, [], (0.0, 100), [])
 sol = solve(prob, Rodas4(), saveat=0.1)
 plot(sol)
 ```
 
+[[1] Wilson, H. R., & Cowan, J. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophysical journal, 12(1), 1-24.](https://www.cell.com/biophysj/fulltext/S0006-3495(72)86068-5)
+
+[[2] Destexhe, A., & Sejnowski, T. J. (2009). The Wilson–Cowan model, 36 years later. Biological cybernetics, 101(1), 1-2.](https://link.springer.com/article/10.1007/s00422-009-0328-3)
+
+
+
 ## 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:
+In this example, we'll construct a model of Parkinson's disease using eight Jansen-Rit Neural Mass Models, based on the work of Liu et al. (2020) [1].
+
+### The Jansen-Rit Neural Mass Model
+
+The Jansen-Rit model [2] is another popular neural mass model that, like the Wilson-Cowan model from [Example 1](#example-1-creating-a-simple-neural-circuit), describes the average activity of neural populations. Each Jansen-Rit unit is defined by the following differential equations:
+
+```math
+\begin{align}
+\frac{dx}{dt} &= y-\frac{2}{\tau}x \\[10pt]
+\frac{dy}{dt} &= -\frac{x}{\tau^2} + \frac{H}{\tau} \left[2\lambda S(\textstyle\sum{jcn}) - \lambda\right]
+\end{align}
+```
+
+where $x$ represents the average postsynaptic membrane potential of the neural population, $y$ is an auxiliary variable, $\tau$ is the membrane time constant, $H$ is the maximum postsynaptic potential amplitude, $\lambda$ determines the maximum firing rate, and $\sum{jcn}$ represents the sum of all synaptic inputs to the population. The sigmoid function $S(x)$ models the population's firing rate response to input and is defined as:
 
 ```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]
+S(x) = \frac{1}{1 + \text{exp}(-rx)}
 ```
 
+where $r$ controls the steepness of the sigmoid, affecting the population's sensitivity to input.
+
+### Setting Up the Model
+
+Let's start by importing the necessary libraries and defining our neural masses:
+
 ```@example Jansen-Rit
 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
 
-```@example Jansen-Rit
+# Convert time units from seconds to milliseconds
 τ_factor = 1000
 
+# Define Jansen-Rit neural masses for different brain regions
 @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 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 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 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]
 ```
-Again, we create a graph and add the Blox as nodes
+
+Here, we've created eight Jansen-Rit neural masses representing different brain regions involved in Parkinson's disease. The `τ_factor` is used to convert time units from seconds (as in the original paper) to milliseconds (Neuroblox's default time unit).
+
+### Building the Circuit
+
+Now, let's create a graph representing our brain circuit:
 
 ```@example Jansen-Rit
 g = MetaDiGraph()
 add_blox!.(Ref(g), blox)
+```
+
+We've created a MetaDiGraph and added our neural masses as nodes. Next, we'll define the connections between these nodes based on the known anatomy of the basal ganglia-thalamocortical circuit.
+
+ModelingToolkit allows us to create parameters that can be passed into the equations symbolically:
 
+```@example Jansen-Rit
+# Define connection strength parameters
 params = @parameters C_Cor=60 C_BG_Th=60 C_Cor_BG_Th=5 C_BG_Th_Cor=5
 ```
-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.
+As an alternative to creating edges with an adjacency matrix (as shown in the previous example), here we demonstrate a different approach by adding edges one by one. In this case, we set the connections specified in Table 2 of Liu et al. [1], although we only implement a subset of the nodes and edges to describe a simplified version of the model:
 
 ```@example Jansen-Rit
+# Add connections between brain regions
 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))
@@ -116,23 +183,27 @@ 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)
 ```
-Now we are ready to build the ModelingToolkit System and apply structural simplification to the equations.
+
+### Creating the Model
+
+Let's build the complete model:
 
 ```@example Jansen-Rit
 @named final_system = system_from_graph(g)
 final_system_sys = structural_simplify(final_system)
 ```
-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).
+
+This creates a differential equations system from our graph representation using ModelingToolkit and symbolically simplifies it for efficient computation.
+
+### Simulating the Model
+
+Lastly, we create the `ODEProblem` for our system, select an algorithm, in this case `Tsit5()`, and simulate 1 second of brain activity.
+
 ```@example Jansen-Rit
 sim_dur = 1000.0 # Simulate for 1 second
-prob = ODEProblem(final_system_sys,
-    [],
-    (0.0, sim_dur))
-```
-We select an algorithm and solve the system
-```@example Jansen-Rit
-alg = Tsit5()
-sol_dde_no_delays = solve(prob, alg, saveat=1)
-plot(sol_dde_no_delays)
+prob = ODEProblem(final_system_sys, [], (0.0, sim_dur))
+sol = solve(prob, Tsit5(), saveat=1)
+plot(sol)
 ```
-In a later tutorial, we will show how to introduce edge delays.
+
+[[1] Liu, C., Zhou, C., Wang, J., Fietkiewicz, C., & Loparo, K. A. (2020). The role of coupling connections in a model of the cortico-basal ganglia-thalamocortical neural loop for the generation of beta oscillations. Neural Networks, 123, 381-392.](https://doi.org/10.1016/j.neunet.2019.12.021)

docs/src/index.md

@@ -1,14 +1,27 @@
 # Neuroblox
 
 ## About
-Neuroblox.jl is designed for computational neuroscience and psychiatry applications. Our tools range from control circuit system identification to brain circuit simulations bridging scales from spiking neurons to fMRI-derived circuits, parameter-fitting models to neuroimaging data, interactions between the brain and other physiological systems, experimental optimization, and scientific machine learning.
+### Overview
+Neuroblox is designed for computational neuroscience and psychiatry applications. Our tools range from brain circuit simulations to control circuit system identification, bridging scales from spiking neurons to fMRI-derived circuits, interactions between the brain and other physiological systems, experimental optimization, and scientific machine learning.
 
-## Description
-Neuroblox.jl is based on a library of modular computational building blocks (“blox”) in the form of systems of symbolic dynamic differential equations that can be combined to describe large-scale brain dynamics.  Once a model is built, it can be simulated efficiently and fit electrophysiological and neuroimaging data.  Moreover, the circuit behavior of multiple model variants can be investigated to aid in distinguishing between competing hypotheses.
-We employ ModelingToolkit.jl to describe the dynamical behavior of blox as symbolic (stochastic/delay) differential equations.  Our libraries of modular blox consist of individual neurons (Hodgkin-Huxley, IF, QIF, LIF, etc.), neural mass models (Jansen-Rit, Wilson-Cowan, Lauter-Breakspear, Next Generation, microcanonical circuits etc.) and biomimetically-constrained control circuit elements.  A GUI designed to be intuitive to neuroscientists allows researchers to build models that automatically generate high-performance systems of numerical ordinary/stochastic differential equations from which one can run stimulations with parameters fit to experimental data.  Our benchmarks show that the increase in speed for simulation often exceeds a factor of 100 as compared to neural mass model implementation by the Virtual Brain (python) and similar packages in MATLAB.  For parameter fitting of brain circuit dynamical models, we use Turing.jl to perform probabilistic modeling, including Hamilton-Monte-Carlo sampling and Automated Differentiation Variational Inference.
+### Features
+Neuroblox is based on a library of modular computational building blocks (“blox”) in the form of systems of symbolic dynamic differential equations, which can be flexibly combined to describe large-scale brain dynamics. Our libraries of modular blox consist of individual neurons (Hodgkin-Huxley, IF, QIF, LIF, etc.), neural mass models (Jansen-Rit, Wilson-Cowan, Lauter-Breakspear, Next Generation, microcanonical circuits, etc.), and biomimetically-constrained control circuit elements.
+
+Once a model is built, it can be simulated efficiently and used to fit electrophysiological and neuroimaging data. Moreover, the circuit behavior of multiple model variants can be investigated to aid in distinguishing between competing hypotheses.
+
+### Interface
+Users can interface with Neuroblox either via Julia code or using a simple drag-and-drop GUI designed to be intuitive to neuroscientists. Both interfaces allow researchers to automatically generate high-performance models from which one can run stimulations with parameters fit to experimental data.
+
+### Performance
+The Neuroblox back-end (Neuroblox.jl) is built using [Julia](https://julialang.org/), an open-source, high-level scripting language designed for high-performance in computation-intensive applications. Our benchmarks show greater than 100x increase in speed over neural mass model implementations using the Virtual Brain (Python) and similar packages in MATLAB.
+
+### Implementation
+Under the hood, we employ ModelingToolkit.jl to describe the dynamical behavior of blox as symbolic (stochastic/delay) differential equations. For parameter fitting of brain circuit dynamical models, we use Turing.jl to perform probabilistic modeling, including Hamilton-Monte-Carlo sampling and Automated Differentiation Variational Inference.
 
 ## Installation
 
+Neuroblox requires a valid [Julia](https://julialang.org/downloads/) installation and [JuliaHub](https://juliahub.com/ui/Home) user account. Using Neuroblox via the GUI does *not* require any programming knowledge, but interested users can learn more about Julia [here](https://julialang.org/learning/).
+
 To install Neuroblox.jl, first add the JuliaHubRegistry and then use the Julia package manager:
 
 ```julia
@@ -22,6 +35,6 @@ Pkg.add("Neuroblox")
 
 ## Licensing
 
-Neuroblox is free for non-commerical and academic use. For full details of the license, please see 
+Neuroblox is free for non-commerical and academic use. For full details of the license, please see
 [the Neuroblox EULA](https://github.com/Neuroblox/NeurobloxEULA). For commercial use, get in contact
 with [email protected].

docs/src/tutorials/QIF_movie.jl

@@ -0,0 +1,28 @@
+using Neuroblox
+using MetaGraphs
+using DifferentialEquations
+using Random
+using Plots
+using Statistics
+
+qif1 = QIFNeuron(name=Symbol("qif1"), I_in=0.5)
+g = MetaDiGraph()
+add_blox!.(Ref(g), [qif1])
+
+@named sys = system_from_graph(g)
+sys = structural_simplify(sys)
+length(unknowns(sys))
+
+prob = ODEProblem(sys, [], (0.0, 1000), [])
+sol = solve(prob, Vern7(), saveat=0.1)
+plot(sol.t,sol[1,:],xlims=(0,1000),xlabel="time (ms)",ylabel="mV")
+pI_in = parameters(sys)[findall(x->contains(string(x),"I_in"),parameters(sys))][1]
+
+I = 0:0.01:0.1
+anim = @animate for i in I
+    @show i
+    prob_new = remake(prob,p=[pI_in => i])
+    sol_new = solve(prob_new, Vern7(), saveat=0.1)
+    plot_qif = plot(sol_new.t,sol_new[1,:],xlims=(0,1000),xlabel="time (ms)",ylabel="mV",label="I_in=$(i)",legend=:top)
+end
+gif(anim, "anim_fps5.gif", fps = 5)