Neuron, Neural Mass and Composite Blox

HHNeuronExci(; I_bg = 0)

Excitatory neuron using the Hodgkin-Huxley formalism [1].

Equations were based on the supplementary material of [2] and the threshold values α and β were based on [3]:

\[\frac{dV}{dt} = - G_Na m^3 h (V-E_\text{Na}) - G_K n^4 (V-E_\text{K}) - G_L (V-E_L) + I_\text{bg} + I_\text{syn} + \text{jcn} \\ \frac{dn}{dt} = ϕ*(α_n(V)*(1-n)-β_n(V)*n) \\ \frac{dm}{dt} = ϕ*(α_m(V)*(1-m)-β_m(V)*m) \\ \frac{dh}{dt} = ϕ*(α_h(V)*(1-h)-β_h(V)*h) \\ α_n(V) = 0.01 \frac{V+34}{1- e ^{-\frac{V+34}{10}}} \\ β_n(V) = 0.125 e^{-\frac{V+44}{80}} \\ α_m(V) = 0.1 \frac{V+30}{1 - e^{-\frac{V+30}{10}}} β_m(V) = 4 e^{-\frac{V+55}{18}} α_h(V) = 0.07 e^{-\frac{V+44}{20}} β_h(V) = \frac{1}{1 + e^{-\frac{V+14}{10}}}\]

Model parameters :

• Based on [3] :
  • G_Na = 52 [mV, Na channel conductance]

- `G_K = 20`  [mV, K channel conductance] 
- `E_Na = 55` [mV, Na channel reversal potential]
- `E_K = -90` [mV, K channel reversal potential]
- `E_L = -60` [mV, leak reversal potential]
- `ϕ = 5` [temperature effect on timescale]

Arguments :

• I_bg [μA, background current]

References:

1. AL Hodgkin, AF Huxley, A quantitative description of membrane current and its application to conduction and excitation in nerve. The J. Physiol., 117, 500–544, 1952.
2. Pathak, A., Brincat, S.L., Organtzidis, H. et al., Biomimetic model of corticostriatal micro-assemblies discovers a neural code. Nat Commun. 2025.
3. XJ Wang, Ionic basis for intrinsic 40 Hz neuronal oscillations. NeuroReport 5, 221–224, 1993.

HHNeuronInhib(; I_bg = 0)

Inhibitory neuron using the Hodgkin-Huxley formalism [1].

Equations were based on the supplementary material of [2] and the threshold values α and β were based on [3]:

\[\frac{dV}{dt} = - G_Na m^3 h (V-E_\text{Na}) - G_K n^4 (V-E_\text{K}) - G_L (V-E_L) + I_\text{bg} + I_\text{syn} + \text{jcn} \\ \frac{dn}{dt} = ϕ*(α_n(V)*(1-n)-β_n(V)*n) \\ \frac{dm}{dt} = ϕ*(α_m(V)*(1-m)-β_m(V)*m) \\ \frac{dh}{dt} = ϕ*(α_h(V)*(1-h)-β_h(V)*h) \\ α_n(V) = 0.01 \frac{V+34}{1- e ^{-\frac{V+34}{10}}} \\ β_n(V) = 0.125 e^{-\frac{V+44}{80}} \\ α_m(V) = 0.1 \frac{V+30}{1 - e^{-\frac{V+30}{10}}} β_m(V) = 4 e^{-\frac{V+55}{18}} α_h(V) = 0.07 e^{-\frac{V+44}{20}} β_h(V) = \frac{1}{1 + e^{-\frac{V+14}{10}}}\]

Model parameters based on [3] : - G_Na = 52 [mV, Na channel conductance] - G_K = 20 [mV, K channel conductance] - E_Na = 55 [mV, Na channel reversal potential] - E_K = -90 [mV, K channel reversal potential] - E_L = -60 [mV, leak reversal potential] - ϕ = 5 [temperature effect on timescale]

Arguments :

• I_bg [μA, background current]

References:

1. AL Hodgkin, AF Huxley, A quantitative description of membrane current and its application to conduction and excitation in nerve. The J. Physiol. 117, 500–544, 1952.
2. Pathak, A., Brincat, S.L., Organtzidis, H. et al., Biomimetic model of corticostriatal micro-assemblies discovers a neural code. Nat Commun. 2025.
3. XJ Wang, Pacemaker Neurons for the Theta Rhythm and Their Synchronization in the Septohippocampal Reciprocal Loop. J. Neurophysiol. 87, 889–900, 2002.

BaxterSensoryNeuron(; C_m=0.13, E_Na=60.0, E_Ca=60.0, E_K=-70.0, E_KCa=-65.0, E_L=-40.0, ...)

Aplysia sensory neuron conductance-based model with kinase-dependent parameter “corner” selection. This neuron model implements the Baxter et al. (1999) sensory neuron equations, including:

• a membrane potential equation with multiple ionic currents (Na, Ca-N, Ca-L, K-A, K-V, K,S-V, K,S-I, K,Ca-F, K,Ca-S, leak),
• a Ca pool with Ca-dependent attenuation of Ca currents,
• fast/slow Ca-driven enhancers for Ca-activated K currents, and
• voltage-dependent gating variables with the paper’s steady-state and time-constant forms.

The model is designed to be paired with the HTR5 mode selector blox: binary inputs PKA and PKC (ideally 0/1) are used to bilinearly interpolate (or exactly select, when binary) between the four Table-2 parameter corners: Control (00), PKA-only (10), PKC-only (01), and 5-HT (11).

Core membrane equation (paper Appendix A1 form):

\[C_m \frac{dV}{dt} = I_{inj} - \sum I_{ion}\]

Arguments:

• C_m [nF]: Membrane capacitance.
• E_Na, E_Ca, E_K, E_KCa, E_L [mV]: Reversal potentials.
• Table-2 corner parameters: conductances and selected time-constant extrema that depend on (PKA, PKC).
• Table-1 fixed maximal conductances: gNa, gCaN, gKA, gKCaF, gL.
• Ca pool and Ca-dependent regulation parameters: K_Ca, τ_Ca, τ_gbr_F, τ_gbr_S, etc.
• Gating parameters: (hA, sA, tAmax, tAmin, ...) for each channel as listed in the code.

Inputs:

• jcn [nA]: Generic Neuroblox synaptic / network current input.
• I_stim [nA]: External injected stimulus current.
• I_bias_extra [nA]: Optional extra bias current.
• PKA [-]: Kinase flag (0/1 recommended).
• PKC [-]: Kinase flag (0/1 recommended).

Outputs:

• V [mV]: Membrane potential.

References:

1. Baxter, D. A., Byrne, J. H. (1999). Serotonergic modulation of two potassium currents in Aplysia sensory neurons. Journal of Neurophysiology, 82(6), 2914–2926. (doi:10.1152/jn.1999.82.6.2914)

TRNNeuron(; C_M=1.0, g_Na=30.0, g_K=25.0, g_NaL=0.0247, g_KL=0.1855, g_ClL=0.49,
          g_TCa=1.75, g_KCa=10.0, g_CAN=0.25, g_NaCAN=0.124, g_KCAN=0.126, E_CAN=-20.0,
          I_NCXmax=1600.0, K_m_Na=87.5, K_m_Ca=1.38, k_sat=0.1, γ_NCX=0.35, F=96485.332, R=8.314, T=310.0,
          P_max=1.25, U_kcc2=0.3, G_glia=8.0, ε_k=0.33, K_bath=4.0, Na_gi=18.0,
          τ=0.001, γ_con=0.0445, RTF=26.64, RTF_Ca=13.32)

Thalamic reticular nucleus (TRN) neuron model with explicit ion concentration dynamics and coupling to an α7 nAChR + ER Ca store module. This neuron implements:

• a conductance-based membrane equation with Na/K/leak/Cl, T-type Ca, KCa, and CAN currents,
• explicit extracellular/intracellular ion concentration state variables (Na, K, Cl, Ca),
• an NCX exchanger current, Na/K pump, KCC2 cotransporter, glial uptake, and diffusion terms,
• Ca absorption, and
• external coupling inputs from an α7 receptor blox: I_α7 (membrane current) and J_ER (Ca flux).

Membrane equation (sign conventions follow the implementation):

\[\frac{dV_m}{dt} = -\frac{1}{C_M}\left(I_{Na}+I_{K}+I_{ClL}+I_T+I_{NCX}+I_{KCa}+I_{CAN}+I_{\alpha7}-I_{app}\right)\]

Arguments:

• C_M [µF/cm²]: Membrane capacitance density.
• g_* [-]: Maximal conductances for the listed currents (units chosen so I = g*(V-E) yields µA/cm²).
• E_CAN [mV]: CAN reversal potential (fixed by the model table).
• NCX parameters: I_NCXmax, K_m_Na, K_m_Ca, k_sat, γ_NCX, F, R, T.
• Pump/transport parameters: P_max, U_kcc2, G_glia, ε_k, K_bath, Na_gi.
• Conversion factors: τ (ms→s), γ_con (current→concentration).
• Nernst factors: RTF, RTF_Ca.

Inputs:

• I_app [µA/cm²]: Applied current.
• jcn [µA/cm²]: Generic Neuroblox synaptic input (added to I_app).
• I_α7 [µA/cm²]: α7 current from Alpha7ERnAChR.
• J_ER [mM/ms]: ER→cytosol Ca flux contribution from Alpha7ERnAChR.

Outputs:

• V_m [mV]: Membrane potential.
• Ca_i [mM]: Cytosolic calcium concentration (used by α7/ER module and Ca-dependent currents).

References:

1. King, J. R., et al. (2017). (α7 nAChR / Ca signaling framework used for the coupled receptor module). Molecular Pharmacology. (doi:10.1124/mol.117.111401)

MuscarinicNeuron(; C_M=1.0, g_Na=30.0, g_K=25.0, g_NaL=0.0247, g_KL=0.1855, g_ClL=0.49,
                 g_TCa=1.75, g_KCa=10.0, g_CAN=0.25, g_NaCAN=0.124, g_KCAN=0.126, E_CAN=-20.0,
                 I_NCXmax=1600.0, K_m_Na=87.5, K_m_Ca=1.38, k_sat=0.1, γ_NCX=0.35, F=96485.332, R=8.314, T=310.0,
                 P_max=1.25, U_kcc2=0.3, G_glia=8.0, ε_k=0.33, K_bath=4.0, Na_gi=18.0,
                 τ=0.001, γ_con=0.0445, RTF=26.64, RTF_Ca=13.32)

Thalamic reticular nucleus (TRN)-style conductance-based neuron with explicit ion concentration dynamics, extended with a muscarinic-activated NCM/INCM current input that supports Na/K apportionment for ion conservation.

This neuron is a new neuron type (keeps your original TRNNeuron unchanged) and adds three input ports intended to be wired from MuscarinicR:

• I_NCM : total muscarinic NCM current (added to the membrane current balance),
• I_NaNCM : Na portion of that current (added to Na ion bookkeeping),
• I_KNCM : K portion of that current (added to K ion bookkeeping).

Key design choices:

• NCM is treated as Ca-impermeable, so it does not appear in the Ca equation.
• If I_NCM is large and you are conserving ions explicitly, using I_NaNCM/I_KNCM avoids violating Na/K mass balance by accounting for the ionic carriers of the NCM current.

Membrane equation (sign conventions follow the implementation):

\[\frac{dV_m}{dt} = -\frac{1}{C_M}\left(I_{Na}+I_{K}+I_{ClL}+I_T+I_{NCX}+I_{KCa}+I_{CAN}+I_{\alpha7}+I_{NCM}-I_{app}\right)\]

Arguments:

• C_M [µF/cm²]: Membrane capacitance density.
• g_* [-]: Maximal conductances for Na/K/leak/Cl, T-type Ca, KCa, and CAN currents (units chosen so I = g*(V-E) yields µA/cm²).
• E_CAN [mV]: CAN reversal potential (fixed by the model table).
• NCX parameters: I_NCXmax, K_m_Na, K_m_Ca, k_sat, γ_NCX, F, R, T.
• Pump/transport parameters: P_max, U_kcc2, G_glia, ε_k, K_bath, Na_gi.
• Conversion factors: τ (ms→s), γ_con (µA/cm² → mM/s).
• Nernst factors: RTF (monovalent), RTF_Ca (divalent).

Inputs:

• I_app [µA/cm²]: Applied current.
• jcn [µA/cm²]: Generic Neuroblox synaptic input (added to I_app).
• I_α7 [µA/cm²]: α7 current from an α7 receptor module (if used).
• J_ER [mM/ms]: ER→cytosol Ca flux contribution from an ER Ca store module (if used).
• I_NCM [µA/cm²]: Total muscarinic NCM/INCM current from MuscarinicR.
• I_NaNCM [µA/cm²]: Na portion of I_NCM for Na ion bookkeeping.
• I_KNCM [µA/cm²]: K portion of I_NCM for K ion bookkeeping.

Outputs:

• V_m [mV]: Membrane potential.
• Ca_i [mM]: Cytosolic calcium concentration.

References:

1. Fransen, E., Alonso, A. A., & Hasselmo, M. E. (2002). Simulations of the Role of the Muscarinic-Activated Calcium-Sensitive Nonspecific Cation Current INCM in Entorhinal Neuronal Activity during Delayed Matching Tasks. Journal of Neuroscience, 22(3), 1081-1097.

VTADANeuron(; C_m=1.0, g_Na=120.0, g_K=36.0, g_L=0.3, E_Na=50.0, E_K=-77.0, E_L=-54.4, I_bias=0.0)

VTA dopaminergic (DA) neuron model for use with β2-nAChR receptors.

Implements a Hodgkin-Huxley-type neuron with sodium, potassium, and leak currents, plus an input for nAChR-mediated current. Uses unified outward-positive current convention.

\[C_m \frac{dV}{dt} = I_{app} - I_{Na} - I_K - I_L - I_{ACh}\]

Arguments:

• C_m [µF/cm²]: Membrane capacitance.
• g_Na, g_K, g_L [mS/cm²]: Maximal conductances.
• E_Na, E_K, E_L [mV]: Reversal potentials.
• I_bias [µA/cm²]: Constant bias current.

Inputs:

• I_app [µA/cm²]: Applied current.
• I_ACh [µA/cm²]: nAChR current (outward positive).

Outputs:

• V [mV]: Membrane potential.

References:

1. Morozova, E. O., et al. (2020). Distinct Temporal Structure of Nicotinic ACh Receptor Activation Determines Responses of VTA Neurons to Endogenous ACh and Nicotine. eNeuro, 7(4). (doi:10.1523/ENEURO.0418-19.2020)

VTAGABANeuron(; C_m=1.0, g_Na=120.0, g_K=36.0, g_L=0.3, E_Na=50.0, E_K=-77.0, E_L=-54.4, I_bias=0.0)

VTA GABAergic neuron model for use with β2-nAChR receptors.

Identical to VTADANeuron but typed as AbstractInhNeuron for network wiring.

Arguments:

• C_m [µF/cm²]: Membrane capacitance.
• g_Na, g_K, g_L [mS/cm²]: Maximal conductances.
• E_Na, E_K, E_L [mV]: Reversal potentials.
• I_bias [µA/cm²]: Constant bias current.

Inputs:

• I_app [µA/cm²]: Applied current.
• I_ACh [µA/cm²]: nAChR current (outward positive).

Outputs:

• V [mV]: Membrane potential.

References:

1. Morozova, E. O., et al. (2020). Distinct Temporal Structure of Nicotinic ACh Receptor Activation Determines Responses of VTA Neurons to Endogenous ACh and Nicotine. eNeuro, 7(4). (doi:10.1523/ENEURO.0418-19.2020)

PINGNeuronExci(name, namespace, C, g_Na, V_Na, g_K, V_K, g_L, V_L, I_ext, τ_R, τ_D)

Create an excitatory neuron from Borgers et al. (2008). The formal definition of this blox is:

\[\frac{dV}{dt} = \frac{1}{C}(-g_{Na} m_{*}^3*h*(V - V_{Na}) - g_K*n^4*(V - V_K) - g_L*(V - V_L) + I_{ext} + jcn) \\ \m_{*} = \frac{a_m(V)}{a_m(V) + b_m(V)} \\ \frac{dn}{dt} = a_n(V)*(1 - n) - b_n(V)*n \\ \frac{dh}{dt} = a_h(V)*(1 - h) - b_h(V)*h \\ \frac{ds}{dt} = \frac{1}{2}*(1 + \text{tanh}(V/10))*(\frac{1 - s}{\tau_R} - \frac{s}{\tau_D})\]

where $jcn$ is any input to the blox. Note that this is a modified Hodgkin-Huxley formalism with an additional synaptic accumulation term. Synapses are added into the $jcn$ term by connecting the postsynaptic neuron's voltage to the presynaptic neuron's output:

\[jcn = w*s*(V_E - V)\]

where $w$ is the weight of the synapse and $V_E$ is the reversal potential of the excitatory synapse.

Inputs:

• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• C: Membrane capacitance (defaults to 1.0).
• g_Na: Sodium conductance (defaults to 100.0).
• V_Na: Sodium reversal potential (defaults to 50.0).
• g_K: Potassium conductance (defaults to 80.0).
• V_K: Potassium reversal potential (defaults to -100.0).
• g_L: Leak conductance (defaults to 0.1).
• V_L: Leak reversal potential (defaults to -67.0).
• I_ext: External current (defaults to 0.0).
• τ_R: Rise time of synaptic conductance (defaults to 0.2).
• τ_D: Decay time of synaptic conductance (defaults to 2.0).

PINGNeuronInhib(name, namespace, C, g_Na, V_Na, g_K, V_K, g_L, V_L, I_ext, τ_R, τ_D)

Create an inhibitory neuron from Borgers et al. (2008). The formal definition of this blox is:

\[\frac{dV}{dt} = \frac{1}{C}(-g_{Na}*m_{*}^3*h*(V - V_{Na}) - g_K*n^4*(V - V_K) - g_L*(V - V_L) + I_{ext} + jcn) \\ \m_{*} = \frac{a_m(V)}{a_m(V) + b_m(V)} \\ \frac{dn}{dt} = a_n(V)*(1 - n) - b_n(V)*n \\ \frac{dh}{dt} = a_h(V)*(1 - h) - b_h(V)*h \\ \frac{ds}{dt} = \frac{1}{2}*(1 + \tanh(V/10))*(\frac{1 - s}{\tau_R} - \frac{s}{\tau_D})\]

where $jcn$ is any input to the blox. Note that this is a modified Hodgkin-Huxley formalism with an additional synaptic accumulation term. Synapses are added into the $jcn$ term by connecting the postsynaptic neuron's voltage to the presynaptic neuron's output:

\[jcn = w*s*(V_I - V)\]

where $w$ is the weight of the synapse and $V_I$ is the reversal potential of the inhibitory synapse.

Inputs:

• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• C: Membrane capacitance (defaults to 1.0).
• g_Na: Sodium conductance (defaults to 35.0).
• V_Na: Sodium reversal potential (defaults to 55.0).
• g_K: Potassium conductance (defaults to 9.0).
• V_K: Potassium reversal potential (defaults to -90.0).
• g_L: Leak conductance (defaults to 0.1).
• V_L: Leak reversal potential (defaults to -65.0).
• I_ext: External current (defaults to 0.0).
• τ_R: Rise time of synaptic conductance (defaults to 0.5).
• τ_D: Decay time of synaptic conductance (defaults to 10.0).

IFNeuron(name, namespace, C, θ, Eₘ, I_in)

Create a basic integrate-and-fire neuron. This follows Lapicque's equation (see Abbott [1], with parameters chosen to match the LIF/QIF neurons implemented as well):

\[\frac{dV}{dt} = \frac{I_{in} + jcn}{C}\]

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.
• C: Membrane capicitance (Default: 1 μF).
• θ: Threshold voltage (Default: -50 mV).
• Eₘ: Resting membrane potential (Default: -70 mV).
• I_in: External current input (Default: 0 μA).
• dtmax: Maximum timestep allowed for adaptive ODE solvers (Default: 0.05 ms).

References:

1. Abbott, L. Lapicque's introduction of the integrate-and-fire model neuron (1907). Brain Res Bull 50, 303-304 (1999).

LIFNeuron(name, namespace, C, θ, Eₘ, I_in)

Create a leaky integrate-and-fire neuron.

This largely follows the formalism and parameters given in Chapter 8 of Sterratt et al. [1], with the following equations:

\[\frac{dV}{dt} = \frac{\frac{-(V-E_m)}{R_m} + I_{in} + jcn}{C} \frac{dG}{dt} = -\frac{1}{\tau}G\]

where $jcn$ is any synaptic input to the blox (presumably a current G from another neuron).

Keyword Arguments:

• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• C: Membrane capicitance (Default: 1 μF).
• Eₘ: Resting membrane potential (Default -70 mV).
• Rₘ: Membrane resistance (Default: 10 kΩ).
• τ: Synaptic time constant (Default 10 ms).
• θ: Threshold voltage (Default: -50 mV).
• E_syn: Synaptic reversal potential (Default: -70 mV).
• G_syn: Synaptic conductance (Default: 0.002 μA/mV).
• I_in: External current input (Default: 0 μA).
• dtmax: Maximum timestep allowed for adaptive ODE solvers (Default: 0.05 ms).

References:

1. Sterratt, D., Graham, B., Gillies, A., & Willshaw, D. (2011). Principles of Computational Modelling in Neuroscience. Cambridge University Press.

QIFNeuron(; name, namespace, kwargs...)

Create a quadratic integrate-and-fire neuron.

The formal definition of this blox is:

\[ \frac{dV}{dt} = ((V - Eₘ)^2 / (Rₘ^2) + I_in + jcn) / C \frac{dG}{dt} = (-1 / τ₂)G + z \frac{dz}{dt} = (-1 / τ₁)z\]

where $jcn$ is any input to the blox.

Keyword Arguments:

• C: Membrane capacitance (μF).
• Rₘ: membrane resistance (kΩ).
• E_syn: Synaptic reversal potential (mV).
• G_syn: Synaptic conductance (mS).
• τ₁: Timescale of decay of synaptic conductance (ms).
• τ₂: Timescale of decay of synaptic spike variable (ms).
• I_in: External current input (μA).
• Eₘ: Resting membrane potential (mV).
• Vᵣₑₛ: Post action potential (mV).
• θ: Threshold potential (mV).
• dtmax: Maximum timestep allowed for adaptive ODE solvers (Default: 0.05 ms).

IzhikevichNeuron(; name, namespace, kwargs...)

Create an Izhikevich neuron, largely following the implementation in Chen and Campbell (2022), with synaptic decay.

The formal definition of this blox is:

\[\frac{dV}{dt} = V(V - α) - w + η + jcn \\ \frac{dw}{dt} = a(bV - w) \\ \frac{dG}{dt} = -(1 / τ)G + z \\ \frac{dz}{dt} = -(1 / τ)z\]

Arguments:

• α: The firing rate parameter (defaults to 0.6215).
• η: Intrinsic current (defaults to 0.12 mA).
• a: Timescale of the recovery variable (defaults to 0.0077 ms).
• b: Sensitivity of the recovery variable to fluctuations in the voltage (defaults to −0.0062).
• θ: Threshold voltage (defaults to 200.0 mV).
• vᵣ: Reset potential after a spike (defaults to -200.0 mV).
• wⱼ: The jump in the recovery variable after a spike (defaults to 0.0189).
• sⱼ: Reset value for the synaptic spike variable after a spike (defaults to 1.2308).
• gₛ: The synaptic conductance (defaults to 1.2308 mS).
• eᵣ: The synaptic reversal potential (defaults to 1.0 mV).
• τ: The synaptic decay time constant (defaults to 2.6 ms).
• dtmax: Maximum timestep allowed for adaptive ODE solvers (Default: 0.01 ms).

References:

1. Izhikevich, E. (2003). Simple model of spiking neurons. IEEE Transactions on Neural Networks, 14(6), 1569–1572.
2. Chen, L., & Campbell, S. A. (2022). Exact mean-field models for spiking neural networks with adaptation. Journal of Computational Neuroscience, 50(4), 445-469.

LIFExciNeuron(; name, namespace, kwargs...)

Create an excitatory leaky integrate-and-fire neuron. This is a model that uses LIF equations for the voltage dynamics, but explicitly models the gating dynamics of glutamate and GABA receptors and their synaptic decay.

The formal definition of this blox is:

\[\frac{dV}{dt} = (-g_L(V - V_L) - S_\text{AMPA, ext} g_\text{AMPA, ext} (V - V_E) - S_\text{GABA} g_\text{GABA}(V - V_I) - S_\text{AMPA} g_\text{AMPA} (V - V_E) - jcn) / C D(S_\text{AMPA}) = - S_\text{AMPA} / τ_\text{AMPA}\\ D(S_\text{GABA}) = - S_\text{GABA} / τ_\text{GABA}\\ D(S_\text{NMDA}) = - S_\text{NMDA} / τ_\text{NMDA, decay} + αx(1 - S_\text{NMDA}) \\ \frac{dx}{dt} = -x / τ_\text{NMDA, rise} \\ \frac{dS_\text{AMPA, ext}}{dt} = -S_\text{AMPA, ext} / τ_\text{AMPA}\]

Keyword Arguments:

• g_L: Leak conductance (mS).
• V_L: Leak reversal potential (mV).
• V_E: Excitatory reversal potential (mV).
• V_I: Inhibitory reversal potential (mV).
• θ: Firing threshold potential (mV).
• V_reset: Reset potential after a spike (mV).
• C: Membrane capacitance (μF).
• τ_AMPA: time constant for the closing of AMPA receptors (ms).
• τ_GABA: time constant for the closing of GABA receptors (ms).
• τNMDAdecay: time constant for the closing of NMDA receptors (ms).
• τNMDArise: time constant for the decay of NMDA current post spike (ms).
• t_refract: Refractory period after a spike (ms).
• α: Scaling for the rise of NMDA current (ms^-1).
• g_AMPA: Synaptic conductance for AMPA glutamate receptors (mS).
• gAMPAext: Synaptic conductance for external current input through AMPA receptors (mS).
• g_GABA: Synaptic conductance for GABA glutamate receptors (mS).
• g_NMDA: Synaptic conductance for NMDA glutamate receptors (mS).
• Mg: Magnesium ion concentration (mM).
• exciscalingfactor: Excitatory scaling factor.
• inhscalingfactor: Inhibitory scaling factor.
• dtmax: Maximum timestep allowed for adaptive ODE solvers (Default: 0.05 ms).

LIFInhNeuron(; name, namespace, kwargs...)

Create an inhibitory leaky integrate-and-fire neuron. This is a model that uses LIF equations for the voltage dynamics, but explicitly models the gating dynamics of AMPA and GABA receptors and their synaptic decay.

The formal definition of this blox is:

\[\frac{dV}{dt} = (1 - \mathbb{1}_\text{refrac}) (-g_L(V - V_L) - S_\text{AMPA, ext} g_\text{AMPA, ext} (V - V_E) - S_\text{GABA} g_\text{GABA}(V - V_I) - S_\text{AMPA} g_\text{AMPA} (V - V_E) - jcn) / C \\ D(S_\text{AMPA}) = - S_\text{AMPA} / τ_\text{AMPA} \\ D(S_\text{GABA}) = - S_\text{GABA} / τ_\text{GABA} \\ D(S_\text{AMPA, ext}) = - S_\text{AMPA, ext} / τ_\text{AMPA}\]

Arguments:

• g_L: Leak conductance (mS).
• V_L: Leak reversal potential (mV).
• V_E: Excitatory reversal potential (mV).
• V_I: Inhibitory reversal potential (mV).
• θ: Firing threshold potential (mV).
• V_reset: Reset potential after a spike (mV).
• C: Membrane capacitance (μF).
• τ_AMPA: time constant for the closing of AMPA receptors (ms).
• τ_GABA: time constant for the closing of GABA receptors (ms).
• t_refract: Refractory period after a spike (ms).
• α: Scaling for the rise of NMDA current (ms^-1).
• g_AMPA: Synaptic conductance for AMPA glutamate receptors (mS).
• gAMPAext: Synaptic conductance for external current input through AMPA receptors (mS).
• g_GABA: Synaptic conductance for GABA glutamate receptors (mS).
• g_NMDA: Synaptic conductance for NMDA glutamate receptors (mS).
• Mg: Magnesium ion concentration (mM).
• exciscalingfactor: Excitatory scaling factor.
• inhscalingfactor: Inhibitory scaling factor.
• dtmax: Maximum timestep allowed for adaptive ODE solvers (Default: 0.05 ms).

NGNMM_Izh(name, namespace, ...)

This is the basic Izhikevich next-gen neural mass as described in [1]. The corresponding connector is set up to allow for connections between masses, but the user must add their own $ \kappa $ values to the connection weight as there is no good way of accounting for this weight within/between regions.

Currently, the connection weights include the presynaptic $ g_s $, but this could be changed.

Equations: To be added once we have a final form that we like here.

Arguments:

• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• include_noise: (default false) determines if the system supports brownian noise
• ζ: (default 0.0) strength of the Brownian noise term (if include_noise == true)
• Other parameters: See reference for full list. Note that parameters are scaled so that units of time are in milliseconds.

Citation:

1. Chen, L., & Campbell, S. A. (2022). Exact mean-field models for spiking neural networks with adaptation. Journal of Computational Neuroscience, 50(4), 445-469.

NGNMM_QIF(name, namespace, ...)

This is the basic QIF next-gen neural mass as described in [1]. This includes the connections via firing rate as described in [1].

Equations: To be added once we have a final form that we like here.

Arguments:

• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• include_noise: (default false) determines if the system supports brownian noise
• A: (default 0.0) strength of the Brownian noise term (if include_noise == true)
• Other parameters: See reference for full list. Note that parameters are scaled so that units of time are in milliseconds.

Citation: Theta-nested gamma bursts by Torcini group.

LinearNeuralMass(name, namespace)

Create standard linear neural mass blox with a single internal state. There are no parameters in this blox. This is a blox of the sort used for spectral DCM modeling. The formal definition of this blox is:

\[\frac{d}{dx} = \sum{jcn}\]

where $jcn$` is any input to the blox.

Arguments:

• name: Options containing specification about deterministic.
• namespace: Additional namespace above name if needed for inheritance.

HarmonicOscillator(name, namespace, ω, ζ, k, h)

Create a harmonic oscillator blox with the specified parameters.

The formal definition of this blox is:

\[\frac{dx}{dt} = y-(2*\omega*\zeta*x)+ k*(2/\pi)*(\text{atan}((\sum{jcn})/h) \\ \frac{dy}{dt} = -(\omega^2)*x\]

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.
• ω: Base frequency. Note the default value is scaled to give oscillations in milliseconds to match other blocks.
• ζ: Damping ratio.
• k: Gain.
• h: Threshold.

JansenRit(name, namespace, τ, H, λ, r, cortical, delayed)

Create a Jansen Rit blox as described in Liu et al. The formal definition of this blox is:

\[\frac{dx}{dt} = y-\frac{2}{\tau}x \\ \frac{dy}{dt} = -\frac{x}{\tau^2} + \frac{H}{\tau} [\frac{2\lambda}{1 + e^{-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. Defaults to 1 for cortical regions, 14 for subcortical.
• H: See equation for use. Defaults to 0.02 for both cortical and subcortical regions.
• λ: See equation for use. Defaults to 5 for cortical regions, 400 for subcortical.
• r: See equation for use. Defaults to 0.15 for cortical regions, 0.1 for subcortical.
• cortical: Boolean to determine whether to use cortical or subcortical parameters. Specifying any of the parameters above will override this.
• delayed: Boolean to indicate whether states are delayed

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.

WilsonCowan(name, namespace, τ_E, τ_I, a_E, a_I, c_EE, c_IE, c_EI, c_II, θ_E, θ_I, η)

Create a standard Wilson Cowan blox. The formal definition of this blox is:

\[\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}))} \\ \frac{dI}{dt} = \frac{-I}{\tau_I} + \frac{1}{1 + e^{-a_I*(c_{EI}*E - c_{II}*I) - \theta_I}}\]

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.
• Others: See equation for use.

LarterBreakspear(name, namespace, ...)

Create a Larter Breakspear blox described in Endo et al. For a full list of the parameters used see the reference. If you need to modify the parameters, see Chesebro et al. and van Nieuwenhuizen et al. for physiological ranges.

Arguments:

• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• Other parameters: See reference for full list. Note that parameters are scaled so that units of time are in milliseconds.

Citations:

1. Endo H, Hiroe N, Yamashita O. Evaluation of Resting Spatio-Temporal Dynamics of a Neural Mass Model Using Resting fMRI Connectivity and EEG Microstates. Front Comput Neurosci. 2020 Jan 17;13:91. doi: 10.3389/fncom.2019.00091.
2. Chesebro AG, Mujica-Parodi LR, Weistuch C. Ion gradient-driven bifurcations of a multi-scale neuronal model. Chaos Solitons Fractals. 2023 Feb;167:113120. doi: 10.1016/j.chaos.2023.113120.
3. van Nieuwenhuizen, H, Chesebro, AG, Polis, C, Clarke, K, Strey, HH, Weistuch, C, Mujica-Parodi, LR. Ketosis regulates K+ ion channels, strengthening brain-wide signaling disrupted by age. Preprint. bioRxiv 2023.05.10.540257; doi: https://doi.org/10.1101/2023.05.10.540257.

KuramotoOscillator(name, namespace, ω, ζ, include_noise=false)

Simple implementation of the Kuramoto oscillator as described in the original paper [1]. Useful for general models of synchronization and oscillatory behavior. The general form of the Kuramoto oscillator is given by: Equations:

\[ \begin{equation} \dot{\theta_i} = \omega_i + \frac{1}{N}\sum_{j=1}^N{K_{i, j}\text{sin}(\theta_j - \theta_i)} \end{equation}\]

Where this describes the connection between regions $i$ and $j$. An alternative form which includes a noise term for each region is also provided, taking the form:

\[ \begin{equation} \dot{\theta_i} = \omega_i + \zeta dW_i \frac{1}{N}\sum_{j=1}^N{K_{i, j}\text{sin}(\theta_j - \theta_i)} \end{equation}\]

where $W_i$ is a Wiener process and $\zeta_i$ is the noise strength.

Keyword arguments:

• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• include_noise: (default false) determines if brownian noise is included in the dynamics of the blox.
• Other parameters: See reference for full list. Note that parameters are scaled so that units of time are in milliseconds. Default parameter values are taken from [2].

Citations:

1. Kuramoto, Y. (1975). Self-entrainment of a population of coupled non-linear oscillators. In: Araki, H. (eds) International Symposium on Mathematical Problems in Theoretical Physics. Lecture Notes in Physics, vol 39. Springer, Berlin, Heidelberg. https://doi.org/10.1007/BFb0013365

2. Sermon JJ, Wiest C, Tan H, Denison T, Duchet B. Evoked resonant neural activity long-term dynamics can be reproduced by a computational model with vesicle depletion. Neurobiol Dis. 2024 Jun 14;199:106565. doi: 10.1016/j.nbd.2024.106565. Epub ahead of print. PMID: 38880431.

VanDerPol(; name, namespace = nothing, θ=1.0, include_noise=false)

Create a neural mass model whose activity variable follows the dynamics of the (stochastic) van der Pol oscillator.

The formal definition of this blox is:

\[\frac{dx}{dt} = y \frac{dy}{dt} = θ(1 - x^2)y - x + ϕ ξ + jcn\]

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.
• include_noise: (default false) controls whether the system includes stochastic noise
• θ: damping strength
• ϕ: strength of the Brownian motion (if include_noise == true)

Generic2dOscillator(name, namespace, ...)

The Generic2dOscillator model is a generic dynamic system with two state variables. The dynamic equations of this model are composed of two ordinary differential equations comprising two nullclines. The first nullcline is a cubic function as it is found in most neuron and population models; the second nullcline is arbitrarily configurable as a polynomial function up to second order. The manipulation of the latter nullcline's parameters allows to generate a wide range of different behaviours.

Equations:

\[ \begin{align} \dot{V} &= d \, \tau (-f V^3 + e V^2 + g V + \alpha W + \gamma I) \\ \dot{W} &= \dfrac{d}{\tau}\,\,(c V^2 + b V - \beta W + a) \end{align}\]

Arguments:

• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• Other parameters: See reference for full list. Note that parameters are scaled so that units of time are in milliseconds.

Citations: FitzHugh, R., Impulses and physiological states in theoretical models of nerve membrane, Biophysical Journal 1: 445, 1961.

Nagumo et.al, An Active Pulse Transmission Line Simulating Nerve Axon, Proceedings of the IRE 50: 2061, 1962.

Stefanescu, R., Jirsa, V.K. Reduced representations of heterogeneous mixed neural networks with synaptic coupling. Physical Review E, 83, 2011.

Jirsa VK, Stefanescu R. Neural population modes capture biologically realistic large-scale network dynamics. Bulletin of Mathematical Biology, 2010.

Stefanescu, R., Jirsa, V.K. A low dimensional description of globally coupled heterogeneous neural networks of excitatory and inhibitory neurons. PLoS Computational Biology, 4(11), 2008).

OUProcess(; name, namespace, μ, σ, τ)

Create a neural mass model whose activity variable follows an Ornstein-Uhlenbeck process.

The formal definition of this blox is:

\[\frac{dx}{dt} = (-x + μ + jcn)/τ + \sqrt{2 / τ} σw\]

Where w is a Brownian variable and jcn is the input to the blox.

Arguments:

• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• τ: relaxation time
• μ: Mean of the OU process
• σ: Strength of the Brownian motion (variance of OUProcess process is τ*σ^2/2)

NGNMM_theta(Cₑ=30, Cᵢ=30, Δₑ=0.5, Δᵢ=0.5, η_0ₑ=10, η_0ᵢ=0, v_synₑₑ=10,
            v_synₑᵢ=-10, v_synᵢₑ=10, v_synᵢᵢ=-10, alpha_invₑₑ=10,
            alpha_invₑᵢ=0.8, alpha_invᵢₑ=10, alpha_invᵢᵢ=0.8, kₑₑ=0, kₑᵢ=0.5,
            kᵢₑ=0.65, kᵢᵢ=0)

A next-generation neural mass model of coupled theta neuron populations. Parameter description and their values are reported in [1]. Each mass consists of a population of two neurons a and b, coupled using different synaptic terms $g$. The entire expression of these is given by:

\[ \frac{a_e}{dt} = \frac{1}{C_e}(b_e*(a_e-1) - (\Delta_e/2)*((a_e+1)^2-b_e^2) - \eta_{0e}*b_e*(a_e+1) - (v_{syn, ee}*g_{ee}+v_{syn, ei}*g_{ei})*(b_e*(a_e+1)) - (g_{ee}/2+g_{ei}/2)*(a_e^2-b_e^2-1)) \\ \frac{b_e}{dt} = \frac{1}{C_e}*((b_e^2-(a_e-1)^2)/2 - \Delta_e*b_e*(a_e+1) + (\eta_{0e}/2)*((a_e+1)^2-b_e^2) + (v_{syn, ee}*(g_{ee}/2)+v_{syn, ei}*(g_{ei}/2))*((a_e+1)^2-b_e^2) - a_e*b_e*(g_{ee}+g_{ei})) \\ \frac{a_i}{dt} = \frac{1}{C_i}(b_i*(a_i-1) - (\Delta_i/2)*((a_i+1)^2-b_i^2) - \eta_{0i}*b_i*(a_i+1) - (v_{syn, ie}*g_{ie}+v_{syn, ii}*g_{ii})*(b_i*(a_i+1)) - (g_{ie}/2+g_{ii}/2)*(a_i^2-b_i^2-1)) \\ \frac{b_i}{dt} = \frac{1}{C_i}*((b_i^2-(a_i-1)^2)/2 - \Delta_i*b_i*(a_i+1) + (\eta_{0i}/2)*((a_i+1)^2-b_i^2) + (v_{syn, ie}*(g_{ie}/2)+v_{syn, ii}*(g_{ii}/2))*((a_i+1)^2-b_i^2) - a_i*b_i*(g_{ie}+g_{ii})) \\ \frac{g_ee}{dt} = \alpha_{inv, ee} (\frac{k_{ee}}{C_e \pi} \frac{1-a_e^2-b_e^2}{(1+2*a_e+a_e^2+b_e^2)} - g_{ee}) \\ \frac{g_ei}{dt} = \alpha_{inv, ei} (\frac{k_{ei}}{C_i \pi} \frac{1-a_i^2-b_i^2}{(1+2*a_i+a_i^2+b_i^2)} - g_{ei}) \\ \frac{g_ie}{dt} = \alpha_{inv, ie} (\frac{k_{ie}}{C_e \pi} \frac{1-a_e^2-b_e^2}{(1+2*a_e+a_e^2+b_e^2)} - g_{ie}) \\ \frac{g_ii}{dt} = \alpha_{inv, ii} (\frac{k_{ii}}{C_i \pi} \frac{1-a_i^2-b_i^2}{(1+2*a_i+a_i^2+b_i^2)} - g_{ii})\]

Alternatively this blox is called by NextGenerationEI(), but this is deprecated and will be removed in future updates.

Citations:

1. Byrne Á, O'Dea RD, Forrester M, Ross J, Coombes S. Next-generation neural mass and field modeling. J Neurophysiol. 2020 Feb 1;123(2):726-742. doi: 10.1152/jn.00406.2019.

NGNMM_theta(Cₑ=30, Cᵢ=30, Δₑ=0.5, Δᵢ=0.5, η_0ₑ=10, η_0ᵢ=0, v_synₑₑ=10,
            v_synₑᵢ=-10, v_synᵢₑ=10, v_synᵢᵢ=-10, alpha_invₑₑ=10,
            alpha_invₑᵢ=0.8, alpha_invᵢₑ=10, alpha_invᵢᵢ=0.8, kₑₑ=0, kₑᵢ=0.5,
            kᵢₑ=0.65, kᵢᵢ=0)

A next-generation neural mass model of coupled theta neuron populations. Parameter description and their values are reported in [1]. Each mass consists of a population of two neurons a and b, coupled using different synaptic terms $g$. The entire expression of these is given by:

\[ \frac{a_e}{dt} = \frac{1}{C_e}(b_e*(a_e-1) - (\Delta_e/2)*((a_e+1)^2-b_e^2) - \eta_{0e}*b_e*(a_e+1) - (v_{syn, ee}*g_{ee}+v_{syn, ei}*g_{ei})*(b_e*(a_e+1)) - (g_{ee}/2+g_{ei}/2)*(a_e^2-b_e^2-1)) \\ \frac{b_e}{dt} = \frac{1}{C_e}*((b_e^2-(a_e-1)^2)/2 - \Delta_e*b_e*(a_e+1) + (\eta_{0e}/2)*((a_e+1)^2-b_e^2) + (v_{syn, ee}*(g_{ee}/2)+v_{syn, ei}*(g_{ei}/2))*((a_e+1)^2-b_e^2) - a_e*b_e*(g_{ee}+g_{ei})) \\ \frac{a_i}{dt} = \frac{1}{C_i}(b_i*(a_i-1) - (\Delta_i/2)*((a_i+1)^2-b_i^2) - \eta_{0i}*b_i*(a_i+1) - (v_{syn, ie}*g_{ie}+v_{syn, ii}*g_{ii})*(b_i*(a_i+1)) - (g_{ie}/2+g_{ii}/2)*(a_i^2-b_i^2-1)) \\ \frac{b_i}{dt} = \frac{1}{C_i}*((b_i^2-(a_i-1)^2)/2 - \Delta_i*b_i*(a_i+1) + (\eta_{0i}/2)*((a_i+1)^2-b_i^2) + (v_{syn, ie}*(g_{ie}/2)+v_{syn, ii}*(g_{ii}/2))*((a_i+1)^2-b_i^2) - a_i*b_i*(g_{ie}+g_{ii})) \\ \frac{g_ee}{dt} = \alpha_{inv, ee} (\frac{k_{ee}}{C_e \pi} \frac{1-a_e^2-b_e^2}{(1+2*a_e+a_e^2+b_e^2)} - g_{ee}) \\ \frac{g_ei}{dt} = \alpha_{inv, ei} (\frac{k_{ei}}{C_i \pi} \frac{1-a_i^2-b_i^2}{(1+2*a_i+a_i^2+b_i^2)} - g_{ei}) \\ \frac{g_ie}{dt} = \alpha_{inv, ie} (\frac{k_{ie}}{C_e \pi} \frac{1-a_e^2-b_e^2}{(1+2*a_e+a_e^2+b_e^2)} - g_{ie}) \\ \frac{g_ii}{dt} = \alpha_{inv, ii} (\frac{k_{ii}}{C_i \pi} \frac{1-a_i^2-b_i^2}{(1+2*a_i+a_i^2+b_i^2)} - g_{ii})\]

Alternatively this blox is called by NextGenerationEI(), but this is deprecated and will be removed in future updates.

Citations:

1. Byrne Á, O'Dea RD, Forrester M, Ross J, Coombes S. Next-generation neural mass and field modeling. J Neurophysiol. 2020 Feb 1;123(2):726-742. doi: 10.1152/jn.00406.2019.

Matrisome(; t_event=180.0)

A model of matrisome, which is part of the Striatum blox, with discrete dynamics. At t_event time during simulation it will store the Cortical input to the Striatum and the internal activity of the Striatum, which can be later used to apply a plasticity rule in the corticostriatal connection.

References:

1. Pathak, A., Brincat, S.L., Organtzidis, H. et al. Biomimetic model of corticostriatal micro-assemblies discovers a neural code. Nat Commun 2025.

See also Striatum

Striosome()

A model of striosome, which is part of the Striatum blox, with discrete dynamics. It calculates a discretized estimate of the internal activity of a Striatum component, which can be used in learning and decision-making tasks, e.g. by comparing the activities of multiple striatal modules and choosing the highest one as the task action.

References:

1. Pathak, A., Brincat, S.L., Organtzidis, H. et al. Biomimetic model of corticostriatal micro-assemblies discovers a neural code. Nat Commun 2025.

See also Striatum

SNc(; κ_DA=1, DA_reward=10, λ_DA=0.33, t_event=90.0)

A model of Substantia Nigra core (SNc) with discrete dynamics. It represents the average activity of neurons within SNc and releases dopamine in an activity-based manner which in turn can modulate the plasticity in connection with HebbianModulationPlasticity.

Arguments :

• κ_DA : Baseline activity of dopaminergic neurons.
• DA_reward : Amount of dopamine released when a reward is received.
• λ_DA : Gain term on the cumulative activity of the source populations that project to SNc.
• t_event : Timepoint during simulation at which the cumulative input to SNc is stored in order to then be used to modulate plasticity in HebbianModulationPlasticity connections.

References:

1. Pathak, A., Brincat, S.L., Organtzidis, H. et al. Biomimetic model of corticostriatal micro-assemblies discovers a neural code. Nat Commun 2025.

TAN(; κ=100, λ=1, rng=Xoshiro(rand(Int)))

A model of TAN (Tonically Active Neurons), which is part of the Striatum, with discrete dynamics that projects to each striatal module of a model to add a random activity value to striatal Matrisome.

Arguments :

• κ : Baseline activity of TAN neurons.
• λ : Gain term on the cumulative activity of the source populations that project to TAN. Typically these are Striosome components of each Striatum module in a model.

References:

1. Pathak, A., Brincat, S.L., Organtzidis, H. et al. Biomimetic model of corticostriatal micro-assemblies discovers a neural code. Nat Commun 2025.

See also Striatum, Striosome and Matrisome.

WinnerTakeAll(; N_exci=5, E_syn_exci=0, E_syn_inhib=-70, G_syn_exci=3, G_syn_inhib=3, I_bg=zeros(N_exci), τ_exci=5, τ_inhib=70)

A winner-take-all microcircuit found in neocortex, also known as a lateral inhibition microcircuit [1]. Multiple pyramidal (excitatory) neurons send synapses to a single interneuron (inhibitory) and receive feedback inhibition from that interneuron. The excitatory neurons are HHNeuronExci and and feedback interneuron is HHNeuronInhib and their parameters are based on [2].

Arguments :

• N_exci : Number of excitatory neurons.
• E_syn_exci [mV, reversal potential for AMPA receptors]
• E_syn_inhib [mV, reversal potential for GABA A receptors]
• G_syn_exci [mV, AMPA receptor conductance]
• G_syn_inhib [mV, GABA A receptor conductance]
• I_bg : [μA] Background current applied to excitatory neurons. If it is a single value then the same current is applied to all N_exci excitatory neurons. If it is a Vector then it needs to be of length N_exci and each element is applied to each neuron.
• τ_exci : [ms, decay time constant for AMPA receptor conductance]
• τ_inhib : [ms, decay time constant for GABA A receptor conductance]

References:

1. Coultrip, Robert, Richard Granger, and Gary Lynch. “A Cortical Model of Winner-Take-All Competition via Lateral Inhibition.” Neural Networks 5, no. 1 (January 1, 1992): 47-54.
2. Pathak, A., Brincat, S.L., Organtzidis, H. et al. Biomimetic model of corticostriatal micro-assemblies discovers a neural code. Nat Commun 2025.

Cortical(; N_wta=20, N_exci=5, E_syn_exci=0, E_syn_inhib=-70, G_syn_exci=3, G_syn_inhib=3,, G_syn_ff_inhib=3.5, I_bg_ar=0, τ_exci=5, τ_inhib=70)

A component representing a layer of neocortex. It contains multiple lateral inhibition WinnerTakeAll microcircuits and a single feedforward inhibitory neuron. Excitatory neurons across WinnerTakeAll microcircuits are connected together using a hypergeometric rule [1]. The feedforward neuron inhibits each one of the N_wta local interneurons. The excitatory neurons are HHNeuronExci and the feedforward and feedback interneurons are HHNeuronInhib and their parameters are based on [2].

Arguments :

• N_wta : Number of WinnerTakeAll microcircuits in the cortical circuit.
• N_exci : Number of excitatory neurons in each one of the N_wta WinnerTakeAll components.
• E_syn_exci [mV, reversal potential for AMPA receptors]
• E_syn_inhib [mV, reversal potential for GABA A receptors]
• G_syn_exci [mV, AMPA receptor conductance]
• G_syn_inhib [mV, GABA A receptor conductance for lateral inhibition]
• G_syn_ff_inhib [mV, GABA A receptor conductance for feedforward inhibition]
• I_bg_ar : [μA] Background current applied to excitatory neurons. If it is a single value then the same current is applied to all [N_wta x N_exci] excitatory neurons. If it is a Vector then it needs to be of length N_exci and each element is applied to each neuron of every WinnerTakeAll microcircuit.
• τ_exci : [ms, decay time constant for AMPA receptor conductance]
• τ_inhib : [ms, decay time constant for GABA A receptor conductance]

References:

1. Andrew C. Felch and Richard H. Granger. The hypergeometric connectivity hypothesis: Divergent performance of brain circuits with different synaptic connectivity distributions. Brain Research, 1202:3–13, April 2008.
2. Pathak, A., Brincat, S.L., Organtzidis, H. et al. Biomimetic model of corticostriatal micro-assemblies discovers a neural code. Nat Commun 2025.

Striatum(N_inhib = 25,
        E_syn_inhib=-70,
        G_syn_inhib=1.2,
        I_bg=zeros(N_inhib),
        τ_inhib=70.0
)

A model of the Striatum that contains three components; Matrisome, Striosome and an additional group of HHNeuronInhib neurons representing Medium Spiny Neurons (MSNs). There are no connections between these clusters, however they may receive connections from or project to other blox components in a model, depending on what other regions are part of the model and the connection rules between Striatum and these regions.

Arguments :

• N_inhib : Number of inhibitory neurons.
• E_syn_inhib [mV, reversal potential for GABA A receptors]
• G_syn_inhib [mV, GABA A receptor conductance]
• I_bg : [μA] Background current applied to the additional cluster of inhibitory neurons. If it is a single value then the same current is applied to all N_inhib inhibitory neurons. If it is a Vector then it needs to be of length N_inhib and each element is applied to one neuron.
• τ_inhib : [ms, decay time constant for GABA A receptor conductance]

References:

1. Pathak, A., Brincat, S.L., Organtzidis, H. et al. Biomimetic model of corticostriatal micro-assemblies discovers a neural code. Nat Commun 2025.

See also Matrisome, Striosome and HHNeuronInhib.

GPi(; N_inhib = 25,
    E_syn_inhib =-70,
    G_syn_inhib=8,
    I_bg=4*ones(N_inhib),
    τ_inhib=70.0

)

A model of Globus Pallidus Internus (GPi). A group of HHNeuronInhib neurons, which are typically project to the inhibitory neurons in Striatum.

Arguments :

• N_inhib : Number of inhibitory neurons.
• E_syn_inhib [mV, reversal potential for GABA A receptors]
• G_syn_inhib [mV, GABA A receptor conductance]
• I_bg : [μA] Background current applied to all neurons. If it is a single value then the same current is applied to all N_inhib inhibitory neurons. If it is a Vector then it needs to be of length N_inhib and each element is applied to one neuron.
• τ_inhib : [ms, decay time constant for GABA A receptor conductance]

References:

1. Pathak, A., Brincat, S.L., Organtzidis, H. et al. Biomimetic model of corticostriatal micro-assemblies discovers a neural code. Nat Commun 2025.

GPe(; N_inhib = 25,
    E_syn_inhib =-70,
    G_syn_inhib=8,
    I_bg=4*ones(N_inhib),
    τ_inhib=70.0

)

A model of Globus Pallidus Externus (GPe). A group of HHNeuronInhib neurons, which are typically project to the inhibitory neurons in Striatum and to the ones in GPi.

Arguments :

• N_inhib : Number of inhibitory neurons.
• E_syn_inhib [mV, reversal potential for GABA A receptors]
• G_syn_inhib [mV, GABA A receptor conductance]
• I_bg : [μA] Background current applied to all neurons. If it is a single value then the same current is applied to all N_inhib inhibitory neurons. If it is a Vector then it needs to be of length N_inhib and each element is applied to one neuron.
• τ_inhib : [ms, decay time constant for GABA A receptor conductance]

References:

1. Pathak, A., Brincat, S.L., Organtzidis, H. et al. Biomimetic model of corticostriatal micro-assemblies discovers a neural code. Nat Commun 2025.

Thalamus(N_exci = 25,
        E_syn_exci=0,
        G_syn_exci=3,
        I_bg=3*ones(N_exci),
        τ_exci=5.0,
        density=0.0,
        weight=1.0,
        rng=default_rng(),
        connection_matrix=subcortical_connection_matrix(density, N_exci, weight; rng)
)

A Thalamus blox used to model either the Thalamus core or matrix subregions. It contains Nexci HHNeuronExci neurons connected via a connectionmatrix.

Arguments :

• N_exci : Number of excitatory neurons.
• E_syn_exci [mV, reversal potential for AMPA receptors]
• G_syn_exci [mV, AMPA receptor conductance]
• I_bg : [μA] Background current applied to excitatory neurons. If it is a single value then the same current is applied to all N_exci excitatory neurons. If it is a Vector then it needs to be of length N_exci and each element is applied to one neuron.
• τ_exci : [ms, decay time constant for AMPA receptor conductance]
• density : Value in range [0,1]. Connection density that determines the sparseness of the Thalamus circuit.
• weight : Connection weight to be used in all connections made between the excitatory neurons of Thalamus.
• rng : Random number generator. This may be used in the connection_matrix to sample connections based on the density value.
• connection_matrix : The connection rule that is applied between the excitatory neurons. By default subcorticalconnectionmatrix determines whether a connection between every possible pair of [N_exci x N_exci] neurons is made by flipping a coin with a probability of success equal to the density value.

References:

1. Pathak, A., Brincat, S.L., Organtzidis, H. et al. Biomimetic model of corticostriatal micro-assemblies discovers a neural code. Nat Commun 2025.

STN(; N_exci = 25,
    E_syn_exci=0,
    G_syn_exci=3,
    I_bg=3*ones(N_exci),
    τ_exci=5.0
)

A model of the Subthalamic Nucleus (STN). A group of HHNeuronExci neurons, which are typically project to the inhibitory neurons in GPe and to the ones in GPi.

Arguments :

• N_exci : Number of excitatory neurons.
• E_syn_exci [mV, reversal potential for AMPA receptors]
• G_syn_exci [mV, AMPA receptor conductance]
• I_bg : [μA] Background current applied to excitatory neurons. If it is a single value then the same current is applied to all N_exci excitatory neurons. If it is a Vector then it needs to be of length N_exci and each element is applied to one neuron.
• τ_exci : [ms, decay time constant for AMPA receptor conductance]

References:

1. Pathak, A., Brincat, S.L., Organtzidis, H. et al. Biomimetic model of corticostriatal micro-assemblies discovers a neural code. Nat Commun 2025.

ConstantInput(; namespace, namespace, I)

Create a constant input current.

Inputs:

• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• I: External current input (mA).

PoissonSpikeTrain(rate, tspan; name, namespace, N_trains, prob_dt, rng)

Create an input that generates spikes according to a Poisson process.

Inputs:

• rate: A single rate or vector of rates for the Poisson processes.
• tspan: A single timespan tuple or vector of timespan tuples.
• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• N_trains: Number of generated spike trains.
• prob_dt: Timestep of the problem, used to determine the time step for the Poisson process.
• rng: Random number generator.

VoltageClampSource(clamp_schedule::Vector{<:NamedTuple{(:t, :V)}})

A blox to emulate a voltage clamp experiment.

Arguments :

• clampschedule : A vector where each element is a Tuple of the form (t = timevalueA, V = voltagevalueA), representing a schedule for the clamping experiment. Each value corresponds to setting the voltage V of the target neuron to a value voltagevalueA at the timepoint timevalue_A.

PulsesInput(; baseline=0, 
            pulse_amp=1,
            t_start=[0], 
            pulse_width=100, 
)

A source that applies square wave pulses as input currents to neuron blox.

Arguments :

• baseline : Baseline current value that applies when the pulse is off.
• pulse_amp : [μA] Current pulse amplitude.
• tstart : [ms] Vector of values that denote the begin of each pulse in a sequence. there will be as many pulses as elements in tstart.
• pulse_width : [ms] Duration of each pulse in a sequence.

DBS(; name, namespace=nothing, frequency=130.0, amplitude=2.5, pulse_width=0.066, 
    offset=0.0, start_time=0.0, smooth=1e-4)

Create a continuous deep brain stimulation (DBS) stimulus with regular pulses.

Arguments:

• name: Name given to ODESystem object within the blox
• namespace: Additional namespace above name if needed for inheritance
• frequency: Pulse frequency in Hz
• amplitude: Pulse amplitude in arbitrary units
• pulse_width: Duration of each pulse in ms
• offset: Baseline value of the signal between pulses
• start_time: Time delay before stimulation begins in ms
• smooth: Smoothing parameter for pulse transitions, set to 0 for sharp transitions

Returns a DBS stimulus blox that outputs square pulses with specified parameters.

ProtocolDBS(; name, namespace=nothing, frequency=130.0, amplitude=2.5,
              pulse_width=0.066, offset=0.0, start_time=0.0, smooth=1e-4,
              pulses_per_burst=10, bursts_per_block=12, 
              pre_block_time=200.0, inter_burst_time=200.0)

Create a deep brain stimulation (DBS) stimulus consisting of a block of pulse bursts.

Arguments:

• name: Name given to ODESystem object within the blox
• namespace: Additional namespace above name if needed for inheritance
• frequency: Pulse frequency in Hz
• amplitude: Pulse amplitude in arbitrary units
• pulse_width: Duration of each pulse in ms
• offset: Baseline value of the signal between pulses
• start_time: Time delay before stimulation begins in ms
• smooth: Smoothing parameter for pulse transitions, set to 0 for sharp transitions
• pulsesperburst: Number of pulses in each burst
• burstsperblock: Number of bursts in the block
• preblocktime: Time before the block starts in ms
• interbursttime: Time between bursts in ms

Returns a DBS stimulus blox that outputs a block of pulse bursts.

ImageStimulus(data::Union{String, StructArray}; t_stimulus, t_pause)

A blox to emulate the presentation of images during a behavioral task.

Arguments :

• data: Either a path to a CSV file, or a CSV file read in as a StructArray. Each row of the CSV represents a separate image and each column is a pixel value of the corresponding image. The images are flattened from a Matrix to a Vector to fit this format.
• t_stimulus : [ms] Timepoint at which an image is presented. This is a single value and applied to all images in data.
• t_pause : [ms] Timepoint at which an image disappears. This emulates pauses that are often added in behavioral tasks after an image is presented and before the participant makes a choice and/or the next trial begins.

BalloonModel(;name, namespace=nothing, lnκ=0.0, lnτ=0.0)

Create a balloon model blox which computes the hemodynamic responses to some underlying neuronal activity The formal definition of this blox is:

\[ \frac{ds}{dt} = \text{jcn} - \kappa s - \gamma (u-1) \frac{du}{dt} = s \frac{d\nu}{dt} = u - v^{1/\alpha} \frac{dq}{dt} = u E(u, E_0)/E_0 - v^{1/\alpha} q/v\]

where $jcn$ is any input to the blox (represents the neuronal activity).

Includes also the bold measurement function:

\[ \lambda(\nu, q) = V_0 \left[ k_1 (1-q) + k_2 \left( 1 - \frac{q}{v} \right) + k_3 (1-v)\right]\]

Arguments:

• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• lnκ: logarithmic prefactor to signal decay H[1], set to 0 for standard parameter value.
• lnτ: logarithmic prefactor to transit time H[3], set to 0 for standard parameter value.
• lnϵ: logarithm of ratio of intra- to extra-vascular signal

NB: the prefix ln of the variables u, ν, q as well as the parameters κ, τ denotes their transformation into logarithmic space to enforce their positivity. This transformation is considered in the derivates of the model equations below.

Citations:

1. Stephan K E, Weiskopf N, Drysdale P M, Robinson P A, and Friston K J. Comparing Hemodynamic Models with DCM. NeuroImage 38, no. 3 (2007): 387–401. doi: 10.1016/j.neuroimage.2007.07.040
2. Hofmann D, Chesebro A G, Rackauckas C, Mujica-Parodi L R, Friston K J, Edelman A, and Strey H H. Leveraging Julia's Automated Differentiation and Symbolic Computation to Increase Spectral DCM Flexibility and Speed, 2023. doi: 10.1101/2023.10.27.564407