My journey in exploring and understanding Brian2 Codebase ...

Table of contents
- Step 1: Cracking Open the Codebase – Where Do We Even Begin?
- Still Feeling Lost ( Yup felt the same ) ? Let’s Start with an Example
- Wait... How Does Brian2 Even Understand Units?
- But... How Does Brian2 Pull This Off?
- Step 1: Import Everything from NumPy
- Step 2: Override with Safe Unit-Aware Functions
- Back to Our Main Quest: Understanding NeuronGroup
- An Intuitive Analogy: Brian2 as a Restaurant Kitchen
- How Brian2 Manages Variables Like a Pro
- 1. Variable: The Foundation Class
- 2. Constant: Scalar Values That Never Change
- 3. ArrayVariable: The Workhorse for State Variables
- 4. DynamicArrayVariable: For Growing/Shrinking Data
- 5. Subexpression: Computed Variables
- 6. AuxiliaryVariable: Temporary Helper Variables
- 7. VariableView: The User Interface
- 8. Variables: The Container and Manager
- How It All Works Together: A Complete Example
- How Groups Work
- NeuronGroup: Putting It All Together
- How Everything in Brian2 Ties Together
- Stage 1: From Equations to Abstract Code
- Stage 2: Analysis & Optimization
- Stage 3: Code Generation
- Stage 4: Templates — Execution Blueprints
- Stage 5: CodeObject — The Executable Package
- Putting It All Together — One Timestep in Brian2
- The Big Picture: What Are Devices Really?
- Why This Matters
- The Device Architecture: The Master Pattern
- RuntimeDevice: The "Live Theater" Approach
- CPPStandaloneDevice: The "Studio Production" Approach
- 2. The Network Class: The Simulation Controller
- 3. The Scheduling System: Coordinating Updates
- 4. How a Simulation Runs: Step by Step
- 5. Multiple Clocks: Different Timescales
- 6. State Storage and Restoration
- 7. How Network Works Internally
- Putting It All Together: A Complete Example
- MagicNetwork vs Network
- The Magic: Automatic Object Discovery
- The Scope System: Controlling the Magic
- Phew… That Was a Ride!
- 🧩 What’s Left? (A Few Extras)
- Putting It All Together: End-to-End Example
- “Miscellaneous But Mighty” tools
- The Monitor Family Tree
- 1. EventMonitor: The Foundation
- 2. SpikeMonitor: The Specialist
- 3. StateMonitor: Watching Continuous Evolution
- 4. PopulationRateMonitor: The Big Picture View
- The Technical Infrastructure: How Monitors Really Work
- Input Generation
- 1. SpikeGeneratorGroup: The Precision Stimulator
- 2. PoissonGroup: The Random Noise Generator
- 3. PoissonInput: The Efficient Background Noise
- 4. BinomialFunction: The Mathematical Engine
- 5. TimedArray: The Signal Generator
- Namespaces
- Understanding Namespaces: The Address Book Metaphor
- The Namespace Hierarchy: Who Gets Priority?
- Functions in Brian2
- The Function Class: Universal Function Abstraction
- Why Preference Systems Matter
- BrianPreference: The Building Block of Configuration
- The Global prefs Object: Mission Control
Alright… so if you’re reading this, chances are you’ve also strapped yourself in for the wild ride that is Brian2—that brilliant Python library that makes neuron simulations look deceptively simple. But spoiler alert: that tiny run()
function you innocently type? Oh boy, it hides a lot under the hood. Like, “secret underground lair” levels of complexity.
This blog is not just a technical walkthrough—it's a story. A story of how I started with wide-eyed curiosity, got hopelessly lost in a maze of abstractions, and slowly (but triumphantly?) began connecting the dots.
Since I’m a visual learner, expect plenty of diagrams, flowcharts, and annotated code snippets. We’ll uncover how your neat little simulation script gets transformed into fast, optimized code that somehow just works.
I'll also be brutally honest about the parts I found hard to follow, the abstractions I lovingly (or not-so-lovingly) broke down, and what finally clicked for me. If you’re also trying to make sense of this brilliant beast of a codebase—grab a coffee, buckle up, and let’s decode Brian2 together - And yes — this guide is long (about an hour’s worth of reading, phew 😅). But that’s intentional. I didn’t want to write just another high-level blog that skims over the surface. My goal was to give you a deep, intuitive, and practical understanding of how Brian2 really works — beneath the syntax and behind the scenes. :)
Step 1: Cracking Open the Codebase – Where Do We Even Begin?
Alright, let’s do what every explorer must—open the map. Or in our case, the codebase. If you take a peek inside the brian2/
directory, you’ll see it's packed with all the internals that power the magic behind Brian2. But here comes the million-dollar question: Where should we even start?
Do we dive head-first into the mysterious core/
folder that sounds important? Or do we play the classic newbie card and open the sacred __init__.py
hoping for divine guidance?
Been there. Done that
To save you from wandering aimlessly through the maze of modules and magic, I’ve done the legwork—and labeled the starting points that actually make sense. Think of it like me leaving trail markers in the jungle 🐾 so you don’t get eaten by wild abstractions …
Still Feeling Lost ( Yup felt the same ) ? Let’s Start with an Example
Honestly, poking around random files in the codebase isn’t the best way to start—it's like reading a novel by jumping to random chapters. Instead, let’s do something more grounded.
from brian2 import *
# Create a SpikeGeneratorGroup with predetermined spike times
spike_indices = [0, 0, 1, 1, 2] # Which neurons spike
spike_times = [10, 30, 15, 35, 25] * ms # When they spike
# Create the generator
stimulus = SpikeGeneratorGroup(3, spike_indices, spike_times)
# Create target neurons to receive the input
target_neurons = NeuronGroup(3, '''
dv/dt = -v/(10*ms) : volt
''', threshold='v > 0.8*volt', reset='v = 0*volt')
# Connect stimulus to targets
connections = Synapses(stimulus, target_neurons, 'w : volt', on_pre='v += w')
connections.connect(j='i') # One-to-one connection
connections.w = 0.5*volt
run(10*ms) # Run the simulation
Now let’s zoom in on this chunk:
target_neurons = NeuronGroup(3, '''
dv/dt = -v/(10*ms) : volt
''', threshold='v > 0.8*volt', reset='v = 0*volt')
Wait... How Does Brian2 Even Understand Units?
Before diving deeper, there was one question I really needed an answer to:
How on earth does Brian2 understand units like
volt
,ms
, oramp
—and even throw errors when I mess up dimensional consistency?
Let me show you what I mean.
For example, if I mess up and do something like:
dv/dt = -v + I_input
...where v
is in volts and I_input
is in amps, Brian doesn’t just shrug and run it. It hits you with this beauty:
X brian2.units.fundamentalunits.DimensionMismatchError: Inconsistent units in differential equation defining variable 'v': Expression '- v + I_input' uses inconsistent units ('- v' has unit V; 'I_input' has unit A).
That’s right—Brian2 actually performs dimensional analysis at runtime.
If you violate physics, it will let you know.
But... How Does Brian2 Pull This Off?
Normally, Python (and even NumPy) has no idea what a volt or second means. Units are just… words. But Brian2 pulls a genius trick:
It wraps and extends NumPy to become unit-aware. Here's how:
Step 1: Import Everything from NumPy
from numpy import *
This brings in the full NumPy universe: sin
, cos
, exp
, array
, etc.—all your mathematical friends.
Step 2: Override with Safe Unit-Aware Functions
from brian2.units.unitsafefunctions import *
Here comes the magic: Brian2 overrides NumPy’s standard functions with its own unit-checked versions. These replacements have the same names as NumPy’s originals, but come with built-in physics sanity checks.
For example:
sin = wrap_function_dimensionless(np.sin)
This creates a new sin
function that:
Checks that its input is dimensionless (like an angle in radians)
If so, calls the regular
np.sin
If not, raises an error
So , Instead of writing:
import numpy as np
A Brian2 user can write:
import brian2.numpy_ as np
Now when they call np.sin(angle)
, they get the unit-aware version that will catch physics mistakes …. how cool is that :)
Back to Our Main Quest: Understanding NeuronGroup
Okay, time to zoom out and get back on track. You’ll find that Brian2 is full of these clever, subtle systems just like the unit checker—but before we chase all the cool details, let’s first understand the core building blocks.
If we want to truly understand how something like:
NeuronGroup(...)
works, we first need to break down two essential concepts:
Variables in Brian2
Groups, and how they own/manage those variables
An Intuitive Analogy: Brian2 as a Restaurant Kitchen
Let’s make this fun with a metaphor.
Imagine Brian2 is a restaurant kitchen. Here's how the parts map out:
Brian2 Concept | Kitchen Analogy |
Variables | Ingredients and storage containers (flour, salt, bowls) |
VariableOwner | Cabinets or fridges that store ingredients by name |
Group | A cooking station that owns certain ingredients and follows a recipe |
CodeRunner | A chef assigned a specific task (sauce chef, grill master) |
NeuronGroup | A fully equipped station with recipes, ingredients, and staff — ready to cook up simulations! 🍳 |
So, when you define a NeuronGroup
, you're basically saying:
“Hey Brian, set up a fully functional cooking station. Give it some ingredients (variables), a recipe (differential equations, thresholds, resets), and the staff (code runners) to execute it every timestep.”
And behind the scenes, Brian2 orchestrates this beautifully.
How Brian2 Manages Variables Like a Pro
Alright, back to real code now — because let’s face it, Brian2 isn’t flipping burgers at McDonald's 🍔 … it's doing high-performance neural simulation, and doing it beautifully.
Let’s start with one of its most fundamental building blocks: Variables.
Why Do We Even Need a "Variable System"?
Imagine you’re simulating 10,000 neurons. Each has:
A membrane voltage
v
A bunch of ion channels (Na⁺, K⁺, etc.)
Several dynamic synapses connecting to other neurons
Now think:
That’s millions of values to manage efficiently
Some variables grow and shrink (like synapses being added or pruned)
Variables carry physical units (volts, amperes, seconds…)
Some variables are shared across simulation components
All this needs to compile to fast code across CPUs, GPUs, C++, and more
This is where Brian2’s Variable system comes in — smart, modular, and super extensible.
Let’s break it down step-by-step:
1. Variable: The Foundation Class
class Variable(CacheKey):
def __init__(self, name, dimensions=DIMENSIONLESS, owner=None, dtype=None,
scalar=False, constant=False, read_only=False, dynamic=False, array=False):
Think of Variable as a "blueprint" or "metadata container." It doesn't store actual data—it stores information about data.
What it tracks:
Name: What this variable is called (
"v"
,"g_syn"
,"t"
)Dimensions: Physical units (volts, seconds, amperes) remember we discussed this earlier
Owner: Which group "owns" this variable (a neuron group, synapse group, etc.)
Type information: Integer, float, boolean
Properties: Is it constant? Read-only? Scalar vs. array?
Example :
# This describes a membrane voltage variable
v_variable = Variable(
name="v",
dimensions=volt.dim, # This has units of voltage
owner=my_neuron_group, # Belongs to a group of neurons
dtype=float64, # Stored as 64-bit floats
scalar=False, # One value per neuron (array)
constant=False, # Can change during simulation
read_only=False # User can modify it
)
Why this design? Separating metadata from actual data allows Brian2 to:
Generate different code for the same variable on different devices ( we’ll go in depth over this)
Validate operations before they happen
Keep track of units without storing them with every number ( super witty right)
2. Constant: Scalar Values That Never Change
class Constant(Variable):
def __init__(self, name, value, dimensions=DIMENSIONLESS, owner=None):
self.value = value # The actual constant value
Constants store single, unchanging values. Think of mathematical constants like π, or simulation parameters like the number of neurons.
Example:
# Number of neurons in a group
N_const = Constant("N", value=1000, dimensions=DIMENSIONLESS)
# A physical constant
R_membrane = Constant("R_m", value=100e6, dimensions=ohm.dim) # 100 MΩ
Why separate from Variable? Constants are so common and simple that they deserve their own optimized implementation. They can be directly embedded in generated code rather than requiring memory lookups.
3. ArrayVariable: The Workhorse for State Variables
class ArrayVariable(Variable):
def __init__(self, name, owner, size, device, dimensions=DIMENSIONLESS,
dtype=None, constant=False, scalar=False, read_only=False,
dynamic=False, unique=False):
ArrayVariable represents data stored in arrays—the actual "memory" where simulation state lives.
Key concept: Device Abstraction
Notice the device
parameter. ArrayVariable doesn't directly store data—it asks a "device" to manage the actual memory. This allows the same variable to live on:
CPU memory (NumPy arrays)
GPU memory (CUDA arrays)
Generated C++ code
Example:
# Membrane voltages for 1000 neurons
voltages = ArrayVariable(
name="v",
owner=neuron_group,
size=1000, # 1000 neurons
device=current_device, # Handles actual memory
dimensions=volt.dim, # Physical units
dtype=float64,
scalar=False # Array, not single value
)
# Device creates actual array: [-70e-3, -70e-3, -70e-3, ...] (1000 values)
4. DynamicArrayVariable: For Growing/Shrinking Data
class DynamicArrayVariable(ArrayVariable):
def __init__(self, ..., needs_reference_update=False, resize_along_first=False):
self.needs_reference_update = needs_reference_update
DynamicArrayVariable handles data that changes size during simulation.
When you need this:
Synapses: You might add/remove connections during simulation
Spikes: Recording spike times (don't know how many will occur)
Dynamic networks: Neurons that die or are born during simulation
Example:
python# Synaptic weights that can grow as connections are added
weights = DynamicArrayVariable(
name="w",
size=0, # Start with no connections
dimensions=siemens.dim, # Conductance units
needs_reference_update=True # Code objects need fresh pointers
)
# During simulation:
# Initially: weights = []
# After connecting 100 synapses: weights = [0.1nS, 0.2nS, ..., 0.15nS]
# After adding 50 more: weights = [0.1nS, 0.2nS, ..., 0.15nS, 0.3nS, ...]
The reference update problem: When arrays resize, their memory location might change. If generated code holds old pointers, it crashes. The needs_reference_update
flag tells Brian2 to refresh all pointers before each simulation step.
Why not just use lists? Dynamic arrays maintain the performance and device abstraction benefits while adding resize capability.
5. Subexpression: Computed Variables
pythonclass Subexpression(Variable):
def __init__(self, name, owner, expr, device, dimensions=DIMENSIONLESS,
dtype=None, scalar=False):
self.expr = expr.strip() # The mathematical expression
self.identifiers = get_identifiers(expr) # Variables used in expression
Subexpressions define variables that are calculated from other variables rather than stored directly.
Why use subexpressions?
Memory efficiency: Don't store computed values
Always current: Automatically reflect changes in underlying variables
Code clarity: Give names to complex expressions
Example:
# In a neuron model:
base_variables = {
'g_na': ArrayVariable(...), # Sodium conductance
'g_k': ArrayVariable(...), # Potassium conductance
'g_l': ArrayVariable(...) # Leak conductance
}
# Total conductance (computed, not stored)
g_total = Subexpression(
name="g_total",
expr="g_na + g_k + g_l", # Mathematical expression
dimensions=siemens.dim, # Result has conductance units
owner=neuron_group
)
# Usage in equations:
# "dv/dt = -(v - E_L) * g_total / C_m"
How it works: When code generation happens, Brian2 substitutes the expression directly:
# Instead of generating:
# g_total = g_na + g_k + g_l
# dv_dt = -(v - E_L) * g_total / C_m
# It generates:
# dv_dt = -(v - E_L) * (g_na + g_k + g_l) / C_m
So , Subexpressions act like "macros" that get expanded during code generation, providing the benefits of named intermediate results without runtime overhead.
6. AuxiliaryVariable: Temporary Helper Variables
class AuxiliaryVariable(Variable):
def get_value(self):
raise TypeError(f"Cannot get the value for an auxiliary variable ({self.name}).")
AuxiliaryVariables are temporary variables created during code generation that don't correspond to anything users directly create.
When they're created:
Condition checking:
_cond
variables for threshold conditionsIntermediate calculations: Breaking complex expressions into steps
Index variables:
_idx
for array indexing
Example:
# User writes: neurons.v[neurons.v > threshold] = reset_value
# Brian2 internally creates:
auxiliary_vars = {
'_cond': AuxiliaryVariable('_cond', dtype=bool), # Which neurons to reset
'_idx': AuxiliaryVariable('_idx', dtype=int32), # Array indices
'_reset_val': AuxiliaryVariable('_reset_val', ...) # Reset value
}
# Generated code (conceptually):
# _cond = (v > threshold)
# _idx = where(_cond)
# v[_idx] = reset_value
Why separate class? AuxiliaryVariables have special properties:
Never directly accessible to users
Don't need persistent storage
Only exist during code generation
7. VariableView: The User Interface
pythonclass VariableView:
def __init__(self, name, variable, group, dimensions=None):
self.variable = variable # The actual Variable object
self.group = group # Context for indexing
self.name = name # Name in this context
VariableView is what users actually interact with—it's the "interface" to variables that provides indexing, assignment, and unit handling.
The key insight: Same variable, different views
# A NeuronGroup voltage variable
neurons = NeuronGroup(100, 'dv/dt = -v/tau : volt')
# A Synapse connecting to those neurons
synapses = Synapses(source_neurons, neurons, ...)
# Same underlying variable, different views:
neurons.v # VariableView: direct access to neuron voltages
synapses.v_post # VariableView: access via synaptic indexing
How indexing works:
python# Different ways to access the same underlying array:
neurons.v[0:10] # First 10 neurons (simple slice)
neurons.v['v > -50*mV'] # Neurons above threshold (string condition)
synapses.v_post[0:100] # Postsynaptic voltages for first 100 synapses
VariableView provides a NumPy-like interface while handling units, code generation, and different indexing schemes behind the scenes.
8. Variables: The Container and Manager
class Variables(Mapping):
def __init__(self, owner, default_index="_idx"):
self._variables = {} # name -> Variable mapping
self.indices = defaultdict(...) # name -> index_name mapping
self.device = get_device() # For creating arrays
Variables is a dictionary-like container that manages all variables for a group and handles the complex relationships between them.
Key responsibilities:
Variable Creation and Management:
vars = Variables(neuron_group)
# Create different types of variables:
vars.add_array('v', size=1000, dimensions=volt.dim) # ArrayVariable
vars.add_constant('N', value=1000) # Constant
vars.add_subexpression('I_total', 'I_Na + I_K + I_L', ...) # Subexpression
So components lik Groups
don’t use the individual classes we discusses but access everything via class Variables
…
How It All Works Together: A Complete Example
Let's trace through creating a simple neuron group to see how all these pieces interact:
# 1. User creates a neuron group
neurons = NeuronGroup(1000, '''
dv/dt = -(v - E_L)/tau + I_syn/C_m : volt
I_syn = g_syn * (v - E_syn) : amp
g_syn : siemens
E_L : volt
tau : second
C_m : farad
''')
# 2. Brian2 internally creates Variables container
neurons.variables = Variables(neurons)
# 3. Variables creates different Variable types:
neurons.variables.add_array('v', size=1000, dimensions=volt.dim) # ArrayVariable
neurons.variables.add_array('g_syn', size=1000, dimensions=siemens.dim) # ArrayVariable
neurons.variables.add_constant('N', value=1000) # Constant
neurons.variables.add_subexpression('I_syn', 'g_syn * (v - E_syn)', amp.dim) # Subexpression
# 4. Device allocates actual memory:
# device.arrays['neurons_v'] = np.zeros(1000) # Voltage array
# device.arrays['neurons_g_syn'] = np.zeros(1000) # Conductance array
# 5. User access creates VariableViews:
voltage_view = neurons.v # VariableView(name='v', variable=v_array_var, group=neurons)
# 6. Operations go through VariableView:
neurons.v[:] = -70*mV # Calls voltage_view.__setitem__(slice(None), -70*mV)
# Which calls set_with_index_array()
# Which calls device.get_value(v_array_var)[slice(None)] = -0.07
# 7. String expressions trigger code generation:
active = neurons.v['v > -50*mV'] # Creates auxiliary variables, generates code
Here's a quick demo I coded 🚀
Phew… that was a lot to unpack!
And believe it or not—that was just the initial layer of how Brian2 handles variables. Now let’s level up and explore how Brian2 organizes Groups, which will lead us directly into the powerful concept of the NeuronGroup.
How Groups Work
1. VariableOwner: The Storage Manager
class VariableOwner(Nameable):
"""Mix-in class for accessing arrays by attribute."""
def __getattr__(self, name):
# This is what makes neurons.v work!
if name[-1] == "_":
name = name[:-1]
use_units = False
else:
use_units = True
return self.state(name, use_units)
def state(self, name, use_units=True):
"""Return the variable in a way that supports indexing"""
var = self.variables[name] # Get the Variable object
if use_units:
return var.get_addressable_value_with_unit(name, self) # VariableView with units
else:
return var.get_addressable_value(name, self) # VariableView without units
VariableOwner's Job: Make variable access natural and user-friendly. When you type neurons.v
, it automatically creates a VariableView that lets you treat the variable like a numpy array ( remember we discusses this earlier)
2. Group: The Organizer and Controller
class Group(VariableOwner, BrianObject ): # We'll discuss what a BrainObject is ...
"""Owns variables AND can execute code on them"""
def __init__(self):
# Inherits variable storage from VariableOwner
# Inherits timing/scheduling from BrianObject
self.contained_objects = [] # List of CodeRunners
def _resolve(self, identifier, run_namespace):
"""Figure out what names in equations refer to"""
# Look in my variables first
if identifier in self.variables:
return self.variables[identifier]
# Then look in standard functions/units
# Then look in user's namespace
def run_regularly(self, code):
"""Create a CodeRunner to execute custom code"""
runner = CodeRunner(self, "stateupdate", code=code)
self.contained_objects.append(runner)
return runner
Group's Job:
Organize related variables together
Understand equations and resolve variable names
Contain and manage CodeRunners that operate on the variables
Provide the context for code execution
3. CodeRunner: The Worker
class CodeRunner(BrianObject):
"""Runs code that operates on a Group's variables"""
def __init__(self, group, template, code):
self.group = group # Which group's variables to work with
self.template = template # What kind of operation (stateupdate, threshold, etc.)
self.abstract_code = code # What code to execute
self.codeobj = None # Will hold the compiled code
def before_run(self, run_namespace):
"""Compile the code into a CodeObject"""
# Resolve all variables and functions mentioned in the code
# Generate actual executable code for the target device
self.codeobj = create_runner_codeobj(
group=self.group,
code=self.abstract_code,
template_name=self.template,
# ... more parameters
)
def run(self):
"""Execute the compiled code"""
if self.codeobj:
self.codeobj() # Run the generated code
CodeRunner's Job: Take abstract mathematical code and actually execute it on the group's variables.
NeuronGroup: Putting It All Together
Now let's see how NeuronGroup
orchestrates all these components:
# When you create a NeuronGroup:
neurons = NeuronGroup(1000, '''
dv/dt = -(v - E_L)/tau : volt
I_syn : amp
''', threshold='v > -50*mV', reset='v = -70*mV')
Here's what happens internally:
Step 1: Variable Creation
class NeuronGroup(Group):
def __init__(self, N, model, threshold=None, reset=None):
# Parse the equations
equations = Equations(model)
# Create Variables container (inherits from VariableOwner)
self.variables = Variables(self)
# Create variables for each equation
for eq in equations.values():
if eq.type == DIFFERENTIAL_EQUATION:
self.variables.add_array(eq.varname, size=N, dimensions=eq.dim)
elif eq.type == PARAMETER:
self.variables.add_array(eq.varname, size=N, dimensions=eq.dim)
# etc.
# Add special variables
self.variables.add_constant('N', N)
self.variables.add_arange('i', size=N) # Neuron indices
self.variables.add_array('_spikespace', size=N+1, dtype=np.int32)
Step 2: CodeRunner Creation
# Create specialized CodeRunners for different jobs
# 1. StateUpdater: Integrates differential equations
self.state_updater = StateUpdater(self, method='euler')
self.contained_objects.append(self.state_updater)
# 2. Thresholder: Checks spike condition
if threshold:
self.thresholder = Thresholder(self, threshold)
self.contained_objects.append(self.thresholder)
# 3. Resetter: Resets variables after spikes
if reset:
self.resetter = Resetter(self, reset)
self.contained_objects.append(self.resetter)
Step 3: Each CodeRunner Has Specialized Code
class StateUpdater(CodeRunner):
"""Updates state variables using numerical integration"""
def update_abstract_code(self, run_namespace):
# Convert "dv/dt = ..." into "v += dt * (...)"
# This generates something like:
self.abstract_code = """
v += dt * (-(v - E_L)/tau)
I_syn = g_syn * (v - E_syn)
"""
class Thresholder(CodeRunner):
"""Checks threshold condition and records spikes"""
def update_abstract_code(self, run_namespace):
# Convert "v > -50*mV" into spike detection code
self.abstract_code = f"""
_cond = {self.group.events['spike']}
# Code to store spike indices in _spikespace
"""
class Resetter(CodeRunner):
"""Resets variables for neurons that spiked"""
def update_abstract_code(self, run_namespace):
# Reset code only runs for spiking neurons
self.abstract_code = self.group.event_codes['spike'] # "v = -70*mV"
How Everything in Brian2 Ties Together
Let’s walk through the full lifecycle of what happens during one simulation timestep in Brian2 – from your high-level equation to optimized, executable code running thousands of neurons.
User Creates: Brian2 Creates Internally:
NeuronGroup(1000, equations) → Variables Container
├── v: ArrayVariable(1000 elements)
├── I_syn: ArrayVariable(1000 elements)
├── N: Constant(1000)
└── i: ArrayVariable([0,1,2,...,999])
CodeRunners:
├── StateUpdater (integrates dv/dt)
├── Thresholder (checks v > threshold)
└── Resetter (executes reset code)
User Access: What Actually Happens:
neurons.v → VariableView(variable=v_array, group=neurons)
neurons.v[0:10] → device.get_value(v_array)[0:10]
neurons.v = -70*mV → device.fill_with_array(v_array, -0.07)
Here's a quick demo again … The more you can see it , the better you’ll understand :)
Now I left a quite few things like how does Code Generation happens … We’ll we not explore the full code-gen system here as that’s for the next blog … but to give you a decent understanding of it -Let’s now unpack the magic of Code Generation – the core engine that converts your abstract equations into blazing fast code.
Stage 1: From Equations to Abstract Code
When you write:
eqs = '''
dv/dt = (E_L - v + R*I)/tau : volt
'''
Brian2 converts this into a language-agnostic mathematical expression (abstract code):
v = v + dt * ((E_L - v + R*I)/tau)
This happens via:
translation.py
→ make_statements()
Each line turns into a
Statement
object:
Statement(
var='v',
op='=',
expr='v + dt * ((E_L - v + R*I)/tau)',
dtype=float64,
...
)
Stage 2: Analysis & Optimization
Before generating code, Brian2 analyzes:
What variables are being read or written
Which are arrays vs scalars
Which expressions are constant (can be moved out of loops)
All this happens in:
translation.py → analyse_identifiers()
Then optimisation.py
kicks in to do the magic:
Optimization Examples
# Original code:
v = v + dt * (sin(2*pi*freq*t) - v)/tau
# Optimized code:
_lio_1 = sin(2*pi*freq*t) # loop invariant
v = v + dt * (_lio_1 - v)/tau
Stage 3: Code Generation
Brian2 supports multiple target languages:
Generator | Language |
NumpyCodeGenerator | Python + NumPy |
CythonCodeGenerator | Python + C++ |
CPPCodeGenerator | C++ |
They all inherit from:
base.py → CodeGenerator
What Code Generators Do
They convert abstract Statement
objects into real code. For example:
// Input Statement:
v = v + dt * ((E_L - v + R*I)/tau)
// Output C++:
v = v + dt * ((E_L - v + R*I)/tau);
Each generator uses methods like translate_statement()
and translate_expression()
to handle language-specific quirks.
Stage 4: Templates — Execution Blueprints
Templates wrap generated code with the actual execution logic.
Think of them as reusable scaffolding that tells Brian2:
How to loop over neurons
What variables are needed
When/how to run the code
Example Template Logic:
{# USES_VARIABLES { v, tau, dt, N } #}
{# ITERATE_ALL { _idx } #}
for _idx in range(N):
v[_idx] = v[_idx] + dt * ((-v[_idx])/tau)
Different templates are used for:
stateupdate.pyx
→ Integratingdv/dt
threshold.pyx
→ Detecting spikesreset.pyx
→ Resetting post-spike neuronssynapses.pyx
→ Synaptic updates
Stage 5: CodeObject — The Executable Package
Finally, everything gets wrapped in a CodeObject
.
This object handles:
Compilation (for Cython or C++)
Variable/namespace binding
Execution of the final code
Used by:
CodeRunner.run() → self.codeobj()
How It’s Created
self.codeobj = create_runner_codeobj(
group=self.group,
code=self.abstract_code,
template_name=self.template,
...
)
Putting It All Together — One Timestep in Brian2
each timestep:
┌──────────────────────┐
│ neurons.state_updater│ ➤ Integrate dv/dt using Euler
└──────────────────────┘
↓
┌──────────────────────┐
│ neurons.thresholder │ ➤ Detect spikes (v > threshold)
└──────────────────────┘
↓
┌──────────────────────┐
│ neurons.resetter │ ➤ Reset v for spiking neurons
└──────────────────────┘
Each of these is a specialized CodeRunner, powered by generated CodeObjects, executing highly optimized code translated from your original equations.
So far, we’ve tamed one of the key beasts in Brian2—the NeuronGroup
But we’re not done yet. In fact, we’ve only scratched the surface.
There are two more core players in Brian2’s simulation architecture that we must meet before we can truly say we understand what happens when you call that innocent-looking line:
run(100*ms)
These two players are:
Devices – the where and how of execution
Networks – the orchestrator that pulls it all together
The Big Picture: What Are Devices Really?
Imagine you're a director of a movie. You have a script (your Brian simulation code), but you need to decide how to film it. You could:
Film it live on stage - Everything happens in real-time, you see results immediately, but it's limited by the theater's resources
Film it in a studio with full production - Take time to set up elaborate sets, special effects, optimize everything, then film all at once for maximum quality
Brian's devices work exactly like this choice. They're execution strategies that determine how your simulation gets transformed from Python code into running computation.
# Your Brian code (the "script")
from brian2 import *
# The "filming choice" - which device to use
set_device('runtime') # Live theater approach
# OR
set_device('cpp_standalone') # Studio production approach
neurons = NeuronGroup(1000, 'dv/dt = -v/tau : volt')
# ... more simulation code
Why This Matters
Everything you’ve seen so far—equation parsing, variable containers, statement translation, optimizations, templates—funnels down into code that runs on a device.
The device handles:
How variables are stored (Python arrays, C++ buffers, etc.)
How we compile or interpret code
How we run operations in simulation steps
How final results are returned to the user
In essence, the device is the execution engine for the entire Brian2 simulation.
The Device Architecture: The Master Pattern
Brian uses a Strategy Pattern at its core. Let's examine how this works by looking at the base Device
class:
Every device must implement a specific contract defined by the base Device
class. Think of this like a standardized interface that all movie production companies must follow, regardless of whether they're doing live theater or big-budget films.
Here are the critical methods every device must implement:
class Device:
def add_array(self, var):
"""Create storage for a variable"""
def get_array_name(self, var):
"""Get a unique identifier for a variable"""
def code_object(self, owner, name, abstract_code, variables, ...):
"""Convert Brian's abstract code into executable form"""
def seed(self, seed=None):
"""Set random number generator seed"""
def spike_queue(self, source_start, source_end):
"""Create a mechanism for handling spike propagation"""
The genius is that your simulation code never needs to know which device is running it. Brian acts as a translator, converting your high-level neuron equations into device-specific operations.
RuntimeDevice: The "Live Theater" Approach
Let's start with the simpler RuntimeDevice
because it's easier to understand and will build your intuition for how devices work.
Core Philosophy
The RuntimeDevice
executes everything immediately in Python/NumPy. When you create a NeuronGroup
, arrays are allocated right away in memory. When you call run()
, the simulation loop executes immediately.
The Array Management System
Here's where it gets fascinating. Look at how RuntimeDevice
manages memory:
class RuntimeDevice(Device):
def __init__(self):
super().__init__()
# This is the heart of RuntimeDevice - a WeakKeyDictionary
# that maps Variable objects to actual numpy arrays
self.arrays = WeakKeyDictionary()
Why a WeakKeyDictionary
? This is brilliant engineering. When your NeuronGroup
gets garbage collected, its variables automatically disappear from the device's memory too. No memory leaks!
Let's trace through an example:
# What happens when you write this:
neurons = NeuronGroup(100, 'dv/dt = -v/tau : volt')
# Behind the scenes:
# 1. Brian creates Variable objects for 'v', 'tau', etc.
# 2. Device.add_array() gets called for each variable
# 3. RuntimeDevice creates actual numpy arrays:
def add_array(self, var):
if isinstance(var, DynamicArrayVariable):
# For arrays that can grow (like spike monitors)
if var.ndim == 1:
arr = DynamicArray1D(var.size, dtype=var.dtype)
else:
arr = DynamicArray(var.size, dtype=var.dtype)
else:
# For fixed-size arrays (like state variables)
arr = np.empty(var.size, dtype=var.dtype)
# Store the mapping: Variable -> numpy array
self.arrays[var] = arr
When Brian encounters your differential equation dv/dt = -v/tau
, it needs to convert this into executable code. This is where code objects come in.
Think of a code object ( which we discussed in depth earlier) as a compiled "function" that knows how to update your variables for one timestep. For RuntimeDevice
, these are typically Cython functions that operate on numpy arrays.
# Your equation: dv/dt = -v/tau
# Gets converted to a code object that essentially does:
# v[i] += dt * (-v[i]/tau[i]) # for each neuron i
The RuntimeDevice
generates these code objects immediately and keeps them in memory, ready to execute.
The Simulation Loop
When you call run()
, here's what happens in RuntimeDevice
:
Immediate execution: The simulation loop starts right away
Direct array access: Code objects read/write directly to numpy arrays
Python/Cython speed: Each timestep executes compiled code, but coordination happens in Python
Real-time results: You can examine variables at any point during simulation
CPPStandaloneDevice: The "Studio Production" Approach
Now
for the complex but powerful CPPStandaloneDevice
. This is where Brian becomes a code generator rather than a direct executor.
The Fundamental Paradigm Shift
Instead of executing your simulation, CPPStandaloneDevice
generates a complete C++ program that implements your simulation. It's like having Brian write a custom simulator specifically for your problem.
The Deferred Execution Model
Here's the mind-bending part: when you use CPPStandaloneDevice
, nothing actually executes when you think it does. Instead, everything gets recorded for later.
# When you write this with cpp_standalone:
neurons = NeuronGroup(100, 'dv/dt = -v/tau : volt')
neurons.v = -70*mV
monitor = SpikeMonitor(neurons)
run(10*ms)
# Nothing executes yet! Instead:
# 1. Device records: "create array for neurons.v"
# 2. Device records: "set neurons.v to -70mV"
# 3. Device records: "create spike monitoring code"
# 4. Device records: "run simulation for 10ms"
# Only when build() is called does code generation begin
This is stored in various queues and data structures within the device:
class CPPStandaloneDevice(Device):
def __init__(self):
# Store all the operations to be performed
self.main_queue = [] # Sequence of operations
self.arrays = {} # Variable -> array name mapping
self.static_arrays = {} # Constant data to include
self.code_objects = {} # Generated code fragments
# ... many more tracking structures
The Array Naming System
Since we're generating C++ code, we need unique names for all arrays. The device maintains a sophisticated naming system:
def get_array_name(self, var, access_data=True):
# Examples of generated names:
# neurons.v -> "_array_neurons_v"
# synapses._synaptic_pre -> "_array_synapses__synaptic_pre"
# monitor.t -> "_dynamic_array_monitor_t"
if isinstance(var, DynamicArrayVariable):
if access_data:
return self.arrays[var] # "_array_monitor_t"
else:
return self.dynamic_arrays[var] # "_dynamic_array_monitor_t"
The distinction between access_data=True/False
for dynamic arrays is crucial. Dynamic arrays in C++ are std::vector
objects. Sometimes you want the data pointer (&vector[0]
), sometimes you want the vector object itself (for resizing).
The Template System: Code Generation Engine
This is where things get really sophisticated. Brian uses a template system to generate C++ code. Let's look at how this works by examining one of the templates.
From the stateupdate.cpp
template:
{# USES_VARIABLES { N } #}
{# ALLOWS_SCALAR_WRITE #}
{% extends 'common_group.cpp' %}
{% block maincode %}
//// MAIN CODE ////////////
// scalar code
const size_t _vectorisation_idx = -1;
{{scalar_code|autoindent}}
const int _N = {{constant_or_scalar('N', variables['N'])}};
{{openmp_pragma('parallel-static')}}
for(int _idx=0; _idx<_N; _idx++)
{
// vector code
const size_t _vectorisation_idx = _idx;
{{vector_code|autoindent}}
}
{% endblock %}
This template gets filled in with your specific equations. The {{}}
parts are template variables that get replaced with generated code. Let's trace through what happens:
Your equation:
dv/dt = -v/tau
Brian's analysis: This becomes vector code that updates each neuron
Code generation: The template gets filled with C++ code
Final C++: A complete function that updates all neurons
The generated C++ might look like:
void _run_neurongroup_stateupdater() {
// scalar code (stuff that happens once per timestep)
const size_t _vectorisation_idx = -1;
const int _N = 100; // Number of neurons
#pragma omp parallel for schedule(static)
for(int _idx=0; _idx<_N; _idx++) {
// vector code (stuff that happens for each neuron)
const size_t _vectorisation_idx = _idx;
_array_neurons_v[_idx] += dt * (-_array_neurons_v[_idx] / _array_neurons_tau[_idx]);
}
}
The Build Process: From Python to Executable
When you call device.build
()
or run with build_on_run=True
, here's the complete pipeline:
When you build a standalone simulation, Brian creates a complete C++ project. Here's what the directory looks like:
output/ # Project directory
├── main.cpp # Main program entry point
├── objects.cpp/h # Array declarations and management
├── network.cpp/h # Simulation loop coordination
├── run.cpp/h # Startup and cleanup functions
├── makefile # Build instructions
├── code_objects/ # Generated update functions
│ ├── neurongroup_stateupdater.cpp
│ ├── synapses_pre.cpp
│ └── spikemonitor.cpp
├── brianlib/ # Brian's C++ support library
│ ├── clocks.h
│ ├── dynamic_array.h
│ └── spikequeue.h
├── static_arrays/ # Constant data files
│ └── initial_values.dat
└── results/ # Output files
├── neuron_v_data.dat
└── spike_times.dat
Deep Dive: The Objects File
Let's examine how the objects.cpp
file gets generated, because this is where you can really see the device's sophistication. This file declares all your variables as C++ arrays and sets up their initialization.
From the template objects.cpp
:
//////////////// arrays ///////////////////
{% for var, varname in array_specs | dictsort(by='value') %}
{% if not var in dynamic_array_specs %}
{{c_data_type(var.dtype)}} * {{varname}};
const int _num_{{varname}} = {{var.size}};
{% endif %}
{% endfor %}
//////////////// dynamic arrays 1d /////////
{% for var, varname in dynamic_array_specs | dictsort(by='value') %}
std::vector<{{c_data_type(var.dtype)}}> {{varname}};
{% endfor %}
When you have a NeuronGroup(100, 'v:volt')
, this template generates:
// Generated C++ code
double* _array_neurons_v; // Pointer to voltage array
const int _num_array_neurons_v = 100; // Size of the array
For monitors with dynamic arrays (that grow during simulation):
std::vector<double> _dynamic_array_spikemonitor_t; // Times of spikes
std::vector<int32_t> _dynamic_array_spikemonitor_i; // Neuron indices
The device tracks all this information during the "collection phase" and then generates the appropriate C++ declarations.
The Main Queue: Recording Operations
One of the most fascinating aspects of CPPStandaloneDevice
is how it records operations for later execution. Everything gets stored in the main_queue
:
# When you write this Python code:
neurons.v = -70*mV
# The device records this operation:
self.main_queue.append(('set_by_constant', (array_name, -0.07, False)))
# When you write:
neurons.v[0:10] = [-65, -66, -67, ...]*mV
# The device records:
self.main_queue.append(('set_by_array', (array_name, static_array_name, False)))
During code generation, these operations get converted to C++ initialization code:
// From 'set_by_constant' operation:
#pragma omp parallel for schedule(static)
for(int i=0; i<_num_array_neurons_v; i++) {
_array_neurons_v[i] = -0.07;
}
// From 'set_by_array' operation:
#pragma omp parallel for schedule(static)
for(int i=0; i<_num_static_array_init_v; i++) {
_array_neurons_v[i] = _static_array_init_v[i];
}
Notice how the device automatically adds OpenMP parallelization pragmas for performance!
The Network Coordination System
The most complex part of the standalone device is coordinating the simulation loop. Brian needs to figure out the correct order to execute different operations based on the scheduling system.
Looking at the network.cpp
template, you can see how this works: ( Don’t worry we’ll see what network does in depth next )
void Network::run(const double duration, ...) {
// Compute which clocks need to update
compute_clocks();
// Set time intervals for all clocks
for(std::set<Clock*>::iterator i=clocks.begin(); i!=clocks.end(); i++)
(*i)->set_interval(t, t_end);
// Main simulation loop
while(clock && clock->running() && !Network::_globally_stopped) {
t = clock->t[0]; // Update network time
// Execute all objects that should run at this time
for(size_t i=0; i<objects.size(); i++) {
Clock *obj_clock = objects[i].first;
if (curclocks.find(obj_clock) != curclocks.end()) {
codeobj_func func = objects[i].second;
if (func) // Execute the code object
func();
}
}
// Advance all active clocks
for(std::set<Clock*>::iterator i=curclocks.begin(); i!=curclocks.end(); i++)
(*i)->tick();
clock = next_clocks(); // Find next clocks to update
}
}
This C++ code implements the same scheduling logic that Brian uses in Python, but with much better performance since it's compiled.
2. The Network Class: The Simulation Controller
The Network
class is the central component of Brian2 that manages the simulation. It acts as a container for all objects involved in a simulation and controls their execution.
What the Network Class Does:
Contains simulation objects: Neuron groups, synapse groups, monitors, etc.
Controls simulation flow: Determines when different components should update
Manages simulation time: Keeps track of the current simulation time and advances it
Handles scheduling: Ensures objects are updated in the correct order
Manages state: Can store and restore the state of the simulation
3. The Scheduling System: Coordinating Updates
One of the most important aspects of the Network
class is its scheduling system. This determines the order in which objects are updated during each time step of the simulation.
Default Schedule
default_schedule = [
"start",
"groups",
"thresholds",
"synapses",
"resets",
"end"
]
Each object in Brian2 has a when
attribute that determines at which stage in this schedule it will be updated. For example, neuron groups are typically updated in the "groups" slot, while spike thresholds are checked in the "thresholds" slot.
Extended Schedule with Automatic Slots
Brian2 automatically adds more slots to the schedule:
before_start → start → after_start →
before_groups → groups → after_groups →
before_thresholds → thresholds → after_thresholds →
before_synapses → synapses → after_synapses →
before_resets → resets → after_resets →
before_end → end → after_end
This allows for fine-grained control over execution order.
4. How a Simulation Runs: Step by Step
Let's break down how a simulation run works in Brian2:
Step 1: Preparation
When network.run
(duration)
is called:
Objects are sorted according to the schedule
Each object's
before_run
method is calledClocks are set up to cover the simulation duration
Step 2: Main Simulation Loop
The main simulation loop does the following at each step:
Find the clock(s) with the smallest current time
Set the network time to this clock's time
Update all objects associated with these clocks in their scheduled order
Advance these clocks by their time step (dt)
Repeat until the simulation duration is reached or the simulation is stopped
Step 3: Cleanup
After the simulation is complete:
Each object's
after_run
method is calledProfiling information is collected (if enabled)
Visual Timeline of a Simulation Run
Time →
│
├── t=0ms: Initial state
│ ├── Update "start" objects
│ ├── Update "groups" objects
│ ├── Update "thresholds" objects
│ ├── Update "synapses" objects
│ ├── Update "resets" objects
│ └── Update "end" objects
│
├── t=0.1ms: Next time step
│ ├── Update "start" objects
│ ├── Update "groups" objects
│ ...
...
5. Multiple Clocks: Different Timescales
Brian2 allows different components to operate on different timescales using separate clocks with different time steps (dt). The Network
class handles this by:
Identifying which clocks have the earliest time
Updating only objects associated with those clocks
Advancing those clocks, leaving other clocks unchanged
Repeating the process, which means faster clocks will update more frequently than slower ones
Example with Multiple Clocks
# 0.1ms time step for neurons
neuron_clock = Clock(dt=0.1*ms)
neurons = NeuronGroup(100, model='dv/dt = -v/tau : 1', clock=neuron_clock)
# 1ms time step for recording
monitor_clock = Clock(dt=1*ms)
monitor = StateMonitor(neurons, 'v', record=True, clock=monitor_clock)
# During a 10ms simulation:
# - neurons will update 100 times (every 0.1ms)
# - monitor will update 10 times (every 1ms)
6. State Storage and Restoration
The Network
class allows you to save and restore the state of a simulation, which is useful for:
Creating checkpoints in long simulations
Running different scenarios from the same initial conditions
Analyzing how small changes affect outcomes
Methods for State Management:
store(name): Saves the current state under a given name
restore(name): Restores a previously saved state
get_states(): Returns a dictionary with the current values of all state variables
set_states(values): Sets state variables according to the provided values
7. How Network Works Internally
Let me walk you through the Network's anatomy:
pythonclass Network(Nameable):
def __init__(self, *objs, **kwds):
self.objects = set() # All the objects in our simulation
self.t_ = 0.0 # Current simulation time
self._schedule = None # Order of execution
The Network maintains several crucial pieces of state:
1. Object Collection: The objects
set contains all BrianObjects (neuron groups, synapses, monitors, etc.). Think of this as the conductor's roster of all musicians.
2. Time Management: self.t_
tracks the current simulation time. Unlike real time, this is "biological time" - the time that matters to your neurons.
3. Scheduling: The _schedule
determines the order operations happen within each time step.
The Simulation Loop: Heart of the Network
The most important method is run()
. Here's how it orchestrates everything:
pythondef run(self, duration, report=None, report_period=10*second, ...):
# Step 1: Preparation
self.before_run(namespace) # Set up all objects
# Step 2: Main simulation loop
while running and not stopped:
# Update network time to current clock time
self.t_ = current_time
# Step 3: Execute all active objects in scheduled order
for obj in active_objects:
if obj._clock in current_clocks:
obj.run() # Execute this object's update
# Step 4: Advance clocks
for clock in current_clocks:
clock.tick_forward()
# Step 5: Cleanup
self.after_run()
This loop is the heartbeat of your simulation. Let me explain each step:
Step 1 - Preparation: Like a conductor checking that all musicians have their sheet music, before_run()
ensures every object knows what it needs to do.
Step 2 - Time Synchronization: The Network figures out which "clock" (timekeeper) should advance next. If you have different components running at different time steps, it coordinates them.
Step 3 - Execution: This is where the magic happens. Each object's run()
method executes its core function - neurons integrate their equations, synapses propagate spikes, monitors record data.
Step 4 - Time Advancement: Move forward in time by incrementing the appropriate clocks.
Step 5 - Cleanup: Like putting instruments away after a concert.
Putting It All Together: A Complete Example
Let me show you how network works …
pythonimport brian2 as b2
# Create a simple network
neurons = b2.NeuronGroup(1000, '''
dv/dt = (-v + I)/tau : volt
I : amp
tau : second
''', threshold='v > -50*b2.mV', reset='v = -70*b2.mV')
# Set parameters
neurons.v = -70*b2.mV
neurons.I = '0.5*nA + 0.1*nA*randn()' # Random input current
neurons.tau = 10*b2.ms
# Add synapses
synapses = b2.Synapses(neurons, neurons, 'w : 1', on_pre='v += w*mV')
synapses.connect(p=0.1) # 10% connection probability
synapses.w = 0.5
# Add monitoring
spike_mon = b2.SpikeMonitor(neurons)
voltage_mon = b2.StateMonitor(neurons, 'v', record=[0, 1, 2])
# Create network explicitly
net = b2.Network(neurons, synapses, spike_mon, voltage_mon)
# Look at scheduling before running
print("=== SCHEDULING SUMMARY ===")
print(net.scheduling_summary())
# Run with profiling
print("\n=== RUNNING SIMULATION ===")
net.run(100*b2.ms, profile=True)
# Look at profiling results
print("\n=== PROFILING SUMMARY ===")
print(b2.profiling_summary(net))
But still we did not get run(100*b2.ms, profile=True)
but instead net.run(100*b2.ms, profile=True)
pretty close hmm.. so what’s the catch …
Enter the MagicNetwork
This is where Brian2 pulls a little trick out of its hat: the MagicNetwork.
Let’s break down why it exists.
# Without MagicNetwork, you'd have to do this:
neurons = NeuronGroup(100, 'dv/dt = -v/tau : volt')
synapses = Synapses(neurons, neurons, on_pre='v += 0.1*mV')
monitor = SpikeMonitor(neurons)
# Manually create network and add everything
net = Network(neurons, synapses, monitor)
net.run(100*ms)
This becomes tedious quickly. You have to remember every object you create and manually add it to the network. Miss one object? Your simulation behaves unexpectedly. Add an object twice? Brian complains. For simple scripts and interactive exploration, this manual bookkeeping is a burden.
MagicNetwork solves this by automatically discovering what should be simulated:
# With MagicNetwork, Brian figures it out automatically:
neurons = NeuronGroup(100, 'dv/dt = -v/tau : volt')
synapses = Synapses(neurons, neurons, on_pre='v += 0.1*mV')
monitor = SpikeMonitor(neurons)
# Just run - Brian automatically finds and includes everything!
run(100*ms)
The "magic" is that Brian automatically discovers the neurons, synapses, and monitor you just created and includes them in the simulation.
MagicNetwork vs Network
Let's understand the relationship between MagicNetwork
and regular Network
by thinking of them as two different approaches to the same fundamental task.
Regular Network: The Explicit Approach
A regular Network
is like a concert conductor who works with a predetermined orchestra roster:
class Network(Nameable):
def __init__(self, *objs, **kwds):
self.objects = set() # Explicitly managed collection
for obj in objs:
self.add(obj) # Manual addition required
The regular Network gives you complete control but requires explicit management:
# You control exactly what gets included
neurons1 = NeuronGroup(50, eqs)
neurons2 = NeuronGroup(50, eqs)
synapses = Synapses(neurons1, neurons2)
# Only include neurons1 and synapses, not neurons2
net = Network(neurons1, synapses)
net.run(100*ms) # neurons2 won't be simulated
MagicNetwork: The Automatic Approach
MagicNetwork
is like a smart assistant who looks around the room and automatically figures out who should be in the orchestra:
class MagicNetwork(Network):
def __init__(self):
if MagicNetwork._already_created:
raise ValueError("There can be only one MagicNetwork.")
super().__init__(name="magicnetwork*")
self._previous_refs = set()
Notice the singleton pattern - there can only be one MagicNetwork instance. This is crucial because it needs to maintain a global view of what objects exist.
The Magic: Automatic Object Discovery
The heart of MagicNetwork's "magic" lies in its ability to automatically discover Brian objects. This happens through a sophisticated process of stack frame inspection - essentially, Python introspection that looks at what variables exist in your code.
Stack Frame Inspection: Looking Into Your Code
When you call run()
, MagicNetwork needs to figure out what objects you want to simulate. It does this by examining the "stack frame" - the current state of your Python program's execution:
def get_objects_in_namespace(level):
"""Find all BrianObjects in the current scope"""
objects = set()
frame = inspect.stack()[level + 1][0] # Get the calling frame
# Look through both global and local variables
for _, v in itertools.chain(frame.f_globals.items(), frame.f_locals.items()):
if isinstance(v, BrianObject):
objects.add(weakref.ref(v)) # Store weak references
del frame # Clean up to avoid memory leaks
return objects
This function performs a kind of "mind reading" - it looks at all the variables in your current Python session and finds the ones that are Brian objects.
The Collect Function: Gathering the Orchestra
The collect()
function builds on stack frame inspection to create the final list of objects:
def collect(level=0):
"""Return the list of BrianObjects that will be simulated"""
all_objects = set()
for obj in get_objects_in_namespace(level=level + 1):
obj = obj() # Dereference the weak reference
if obj.add_to_magic_network: # Only include objects that want to be included
gk = BrianObject._scope_current_key
k = obj._scope_key
if gk != k:
continue # Skip objects from different scopes
all_objects.add(obj)
return all_objects
This function implements several important filtering mechanisms:
Object Consent: Only objects with add_to_magic_network = True
are included. This allows Brian objects to opt out of automatic inclusion.
Scope Filtering: The scope system (which we'll explore shortly) ensures that only objects from the current "scope" are included.
The Scope System: Controlling the Magic
One of MagicNetwork's most sophisticated features is its scope system, which allows you to control which objects get automatically included. Think of scopes as different "rooms" - objects in different rooms won't interfere with each other.
How Scopes Work
Every BrianObject gets a scope key when it's created:
class BrianObject(Nameable):
#: Global key value for scope management
_scope_current_key = 0
def __init__(self, ...):
# Objects remember which scope they were created in
self._scope_key = self._scope_current_key
The scope key is like a timestamp - all objects created at the same "time" (in the same scope) get the same key.
Creating New Scopes
You can create a new scope explicitly:
def start_scope():
"""Start a new scope for magic functions"""
BrianObject._scope_current_key += 1
This is incredibly useful for organizing complex simulations:
import brian2 as b2
print("=== SCOPE 1: First Experiment ===")
neurons1 = b2.NeuronGroup(10, 'dv/dt = -v/(10*ms) : volt')
monitor1 = b2.SpikeMonitor(neurons1)
print(f"Objects in scope {neurons1._scope_key}")
b2.run(10*b2.ms) # Simulates neurons1 and monitor1
print(f"Spikes in first experiment: {monitor1.count[:]}")
print("\n=== STARTING NEW SCOPE ===")
b2.start_scope() # Everything created after this is in a new scope
print("\n=== SCOPE 2: Second Experiment ===")
neurons2 = b2.NeuronGroup(5, 'dv/dt = -v/(5*ms) : volt') # Different parameters
monitor2 = b2.SpikeMonitor(neurons2)
print(f"Objects in scope {neurons2._scope_key}")
b2.run(10*b2.ms) # Only simulates neurons2 and monitor2
print(f"Spikes in second experiment: {monitor2.count[:]}")
# neurons1 and monitor1 are ignored in the second run!
The scope system prevents accidental mixing of different experiments in the same Python session.
The Update Magic Objects Process
The heart of MagicNetwork is the _update_magic_objects()
method, which is called every time you run a simulation:
def _update_magic_objects(self, level):
# Step 1: Discover all objects in the current scope
objects = collect(level + 1)
# Step 2: Include contained objects (e.g., StateUpdater inside NeuronGroup)
contained_objects = set()
for obj in objects:
for contained in _get_contained_objects(obj):
contained_objects.add(contained)
objects |= contained_objects
# Step 3: Analyze what we found - new objects vs. previously run objects
some_known = False # Have we seen any of these objects before?
some_new = False # Are any objects completely new?
for obj in objects:
if obj._network == self.id:
some_known = True # This object was in our previous run
elif obj._network is None and obj.invalidates_magic_network:
some_new = True # This is a new object that matters
# Step 4: Decide what to do based on what we found
if some_known and some_new:
# Ambiguous situation - raise an error
raise MagicError("Mix of old and new objects - unclear intent")
elif some_new:
# All objects are new - start fresh
self.t_ = 0.0
self.assign_id() # Get a new network ID
# Step 5: Claim ownership of all objects
for obj in objects:
if obj._network is None:
obj._network = self.id
self.objects = objects
This method implements a sophisticated decision tree that handles different scenarios automatically.
Ambiguity Detection: When Magic Fails
MagicNetwork's most complex logic deals with ambiguous situations. Consider this scenario:
# First simulation
neurons1 = NeuronGroup(10, eqs)
run(100*ms) # neurons1.t is now 100ms
# Later, add more objects
neurons2 = NeuronGroup(5, eqs) # Fresh object, t=0ms
# Now what should run() do?
run(50*ms) # Should this continue from 100ms or start fresh at 0ms?
Phew… That Was a Ride!
And with that… we’ve unpacked a big chunk of how Brian2 works under the hood — from NeuronGroup
creation to code generation, execution strategies, and even the magic behind run()
.
I hope this deep dive helped demystify the inner workings of Brian2 and gave you a stronger intuition for what's really happening when you write those elegant one-liners.
🧩 What’s Left? (A Few Extras)
While this blog covered most of the core ideas, there are still a few important components we haven’t touched on — like:
Preferences: Configure default settings that influence Brian’s behavior.
Monitors: Tools like
StateMonitor
,SpikeMonitor
, andPopulationRateMonitor
to record what’s happening inside your simulation.Other stuff …
These aren’t strictly necessary to understand how the system runs, but they’re incredibly useful once you start building larger and more complex simulations.
You can think of this section as the “Miscellaneous But Mighty” tools — and I’ve added their quick explanations below if you’d like to explore them.
Putting It All Together: End-to-End Example
To close it all out, here’s a final simulation that ties everything together — a full example showing how Brian2 works from start to finish, using the concepts we’ve covered so far.
“Miscellaneous But Mighty” tools
The Monitor Family Tree
Brian2 has four main types of monitors, each designed for different kinds of observations:
EventMonitor → The foundation class that handles any kind of discrete event SpikeMonitor → Specialized for recording spikes (inherits from EventMonitor)
StateMonitor → Records the continuous evolution of state variables over time PopulationRateMonitor → Tracks the instantaneous firing rate of neuron populations
Let me walk you through each one with building complexity.
1. EventMonitor: The Foundation
The EventMonitor is the most general monitoring tool. It can record any discrete event that happens in your simulation, not just spikes.
import numpy as np
from brian2 import *
# Create a simple neuron group
G = NeuronGroup(3, '''
dv/dt = (I - v)/tau : volt
I : volt
tau : second
''', threshold='v > -50*mV', reset='v = -70*mV')
G.tau = 10*ms
G.v = [-70, -60, -50]*mV # Different starting voltages
G.I = [0.1, 0.2, 0.3]*volt # Different input currents
# Create an EventMonitor to watch for threshold crossings
event_mon = EventMonitor(G, 'spike', variables=['v', 'I'])
run(50*ms)
print("Spike times:", event_mon.t)
print("Neuron indices that spiked:", event_mon.i)
print("Voltages at spike times:", event_mon.v)
How EventMonitor Works Internally
When you create an EventMonitor, several important things happen behind the scenes:
Connection to Source: The monitor establishes a connection to your neuron group and specifically watches for the 'spike' event.
Variable References: It creates references to any variables you want to record (like 'v' and 'I' in our example).
CodeRunner Creation: The monitor becomes a CodeRunner that executes code every time the watched event occurs.
Dynamic Arrays: It creates dynamic arrays to store the recorded data - these can grow as needed during the simulation.
Looking at the source code, the EventMonitor sets up these key components:# From the EventMonitor.__init__ method self.variables.add_dynamic_array('i', size=0, ...) # neuron indices self.variables.add_dynamic_array('t', size=0, ...) # spike times # Plus arrays for any additional variables you specify
2. SpikeMonitor: The Specialist
SpikeMonitor is essentially an EventMonitor that's hardcoded to watch for 'spike' events. It's the most commonly used monitor because spikes are the primary currency of neural communication.
# Create a more interesting neuron group
neurons = NeuronGroup(10, '''
dv/dt = (I - v + xi*tau**-0.5)/tau : volt
I : volt
tau : second
''', threshold='v > -50*mV', reset='v = -70*mV')
neurons.tau = 10*ms
neurons.v = -70*mV
neurons.I = '0.5*volt + 0.3*volt*rand()' # Random input currents
# Monitor spikes and membrane voltage at spike times
spike_mon = SpikeMonitor(neurons, variables=['v'])
run(100*ms)
# Access the recorded data
print(f"Total spikes: {spike_mon.num_spikes}")
print(f"Spike trains for neuron 0: {spike_mon.spike_trains()[0]}")
# Plot raster plot
figure(figsize=(10, 4))
plot(spike_mon.t/ms, spike_mon.i, '.k')
xlabel('Time (ms)')
ylabel('Neuron index')
title('Raster plot')
show()
The Spike Train Dictionary
One of SpikeMonitor's most useful features is the spike_trains()
method, which returns a dictionary where each key is a neuron index and each value is an array of that neuron's spike times:
trains = spike_mon.spike_trains()
for neuron_id, spike_times in trains.items():
print(f"Neuron {neuron_id} spiked at: {spike_times}")
This organization makes it easy to analyze individual neuron firing patterns or compute statistics like interspike intervals.
3. StateMonitor: Watching Continuous Evolution
While spike monitors capture discrete events, StateMonitor records the continuous evolution of state variables over time. This is crucial for understanding the dynamics that lead to spikes.
# Create neurons with more complex dynamics
eqs = '''
dv/dt = (I - v + g_e*(E_e - v) + g_i*(E_i - v))/tau : volt
dg_e/dt = -g_e/tau_e : siemens
dg_i/dt = -g_i/tau_i : siemens
I : volt
tau : second
E_e : volt
E_i : volt
tau_e : second
tau_i : second
'''
neurons = NeuronGroup(5, eqs, threshold='v > -50*mV', reset='v = -70*mV')
# Set parameters
neurons.tau = 10*ms
neurons.E_e = 0*mV
neurons.E_i = -80*mV
neurons.tau_e = 5*ms
neurons.tau_i = 10*ms
neurons.v = -70*mV
neurons.I = '0.1*volt*rand()'
# Monitor multiple state variables for first 3 neurons
state_mon = StateMonitor(neurons, ['v', 'g_e', 'g_i'], record=[0, 1, 2])
# Also monitor spikes
spike_mon = SpikeMonitor(neurons)
run(50*ms)
# Plot the membrane voltage evolution
figure(figsize=(12, 8))
for i in range(3):
subplot(3, 1, i+1)
plot(state_mon.t/ms, state_mon.v[i]/mV, label=f'Neuron {i}')
ylabel('Voltage (mV)')
legend()
xlabel('Time (ms)')
title('Membrane Voltage Evolution')
show()
StateMonitor's Indexing System
StateMonitor has a sophisticated indexing system. The record
parameter determines which neurons to monitor:
record=True
: Record all neurons (memory intensive!)record=[0, 2, 4]
: Record specific neuronsrecord=slice(0, 10, 2)
: Record every other neuron from 0 to 10
The recorded data is stored in 2D arrays where:
First dimension: time steps
Second dimension: recorded neurons
print(f"Shape of voltage data: {state_mon.v.shape}")
print(f"Voltage of neuron 1 at all times: {state_mon.v[1]}")
print(f"Voltage of all neurons at time index 100: {state_mon.v[:, 100]}")
StateMonitorView: Convenient Data Access
When you index a StateMonitor (like state_mon[0]
), you get a StateMonitorView that provides convenient access to that neuron's data:
neuron_0_data = state_mon[0]
plot(neuron_0_data.t/ms, neuron_0_data.v/mV)
4. PopulationRateMonitor: The Big Picture View
PopulationRateMonitor gives you a population-level view by computing the instantaneous firing rate of a group of neurons at each time step.
# Create a larger population
N = 100
neurons = NeuronGroup(N, '''
dv/dt = (I - v + sigma*xi*tau**-0.5)/tau : volt
I : volt
sigma : volt
tau : second
''', threshold='v > -50*mV', reset='v = -70*mV')
neurons.tau = 10*ms
neurons.v = -70*mV
neurons.I = '(0.5 + 0.2*rand())*volt' # Variable input
neurons.sigma = 0.1*volt # Noise strength
# Monitor population rate
rate_mon = PopulationRateMonitor(neurons)
spike_mon = SpikeMonitor(neurons)
run(200*ms)
# Plot population rate over time
figure(figsize=(12, 6))
subplot(2, 1, 1)
plot(rate_mon.t/ms, rate_mon.rate/Hz)
ylabel('Population Rate (Hz)')
title('Population Firing Rate')
subplot(2, 1, 2)
plot(spike_mon.t/ms, spike_mon.i, '.k', markersize=1)
ylabel('Neuron Index')
xlabel('Time (ms)')
title('Raster Plot')
tight_layout()
show()
Rate Calculation and Smoothing
The raw population rate can be quite noisy. PopulationRateMonitor provides a smooth_rate()
method for cleaner analysis:
# Smooth the rate with a Gaussian window
smooth_rate = rate_mon.smooth_rate(window='gaussian', width=5*ms)
plot(rate_mon.t/ms, rate_mon.rate/Hz, alpha=0.5, label='Raw rate')
plot(rate_mon.t/ms, smooth_rate/Hz, label='Smoothed rate')
legend()
The Technical Infrastructure: How Monitors Really Work
Now let's dive into the technical details that make all this possible.
CodeRunner Architecture
All monitors inherit from CodeRunner, which is Brian2's system for executing code during simulations. When you create a monitor, you're actually creating a specialized CodeRunner that:
Runs at specific times: Monitors typically run in the 'start' slot for StateMonitor or immediately after spike detection for spike monitors.
Has access to group variables: Through the Variables system, monitors can read any variable from the groups they're monitoring.
Executes template-based code: Brian2 uses templates to generate efficient C++ or Python code for different backends.
The Variables System Integration
Monitors integrate deeply with Brian2's Variables system. Let's see what happens when you create a StateMonitor:
# When you create: StateMonitor(neurons, ['v', 'g_e'], record=[0, 1, 2])
# The monitor creates these internal variables:
self.variables.add_dynamic_array('v', size=(0, 3), ...) # Time x neurons
self.variables.add_dynamic_array('g_e', size=(0, 3), ...)
self.variables.add_dynamic_array('t', size=0, ...) # Time points
self.variables.add_reference('_source_v', neurons, 'v') # Reference to source
Memory Management and Dynamic Arrays
Monitors use dynamic arrays that can grow during simulation. This is more complex than it sounds:
# From StateMonitor.resize() method - called automatically
def resize(self, new_size):
self.variables['N'].set_value(new_size) # Update size counter
self.variables['t'].resize(new_size) # Resize time array
# Resize all monitored variable arrays
for var in self.recorded_variables.values():
var.resize((new_size, self.n_indices)) # 2D: time x neurons
Input Generation
In real neuroscience experiments, researchers use various devices to stimulate neurons: electrodes that inject current, optogenetic light pulses, sensory stimuli, or pharmacological agents. In computational neuroscience, we need digital equivalents that can generate realistic input patterns.
Brian2's input system provides several specialized tools, each designed for different types of stimulation:
SpikeGeneratorGroup: Like a programmable electrode that fires at exact predetermined times
PoissonGroup: Like a noisy input source that generates random spikes with controlled statistics
PoissonInput: Like multiple independent noisy inputs efficiently bundled together
TimedArray: Like a signal generator that can produce any time-varying waveform
BinomialFunction: The mathematical engine that efficiently generates random spike patterns
Let me walk you through each one, building your understanding from the ground up.
1. SpikeGeneratorGroup: The Precision Stimulator
SpikeGeneratorGroup is like having a computer-controlled electrode that can fire spikes at exactly the times you specify. This is perfect when you want precise control over timing.
from brian2 import *
import numpy as np
import matplotlib.pyplot as plt
# Create a SpikeGeneratorGroup with predetermined spike times
spike_indices = [0, 0, 1, 1, 2] # Which neurons spike
spike_times = [10, 30, 15, 35, 25] * ms # When they spike
# Create the generator
stimulus = SpikeGeneratorGroup(3, spike_indices, spike_times)
# Create target neurons to receive the input
target_neurons = NeuronGroup(3, '''
dv/dt = -v/(10*ms) : volt
''', threshold='v > 0.8*volt', reset='v = 0*volt')
# Connect stimulus to targets
connections = Synapses(stimulus, target_neurons, 'w : volt', on_pre='v += w')
connections.connect(j='i') # One-to-one connection
connections.w = 0.5*volt
# Monitor both stimulus and target activity
stimulus_mon = SpikeMonitor(stimulus)
target_mon = SpikeMonitor(target_neurons)
voltage_mon = StateMonitor(target_neurons, 'v', record=True)
run(50*ms)
# Visualize the precise timing
figure(figsize=(12, 8))
subplot(3, 1, 1)
plot(stimulus_mon.t/ms, stimulus_mon.i, 'ro', markersize=8, label='Stimulus spikes')
ylabel('Stimulus\nNeuron')
title('SpikeGeneratorGroup: Precise Timing Control')
legend()
subplot(3, 1, 2)
plot(target_mon.t/ms, target_mon.i, 'bs', markersize=6, label='Target spikes')
ylabel('Target\nNeuron')
legend()
subplot(3, 1, 3)
for i in range(3):
plot(voltage_mon.t/ms, voltage_mon.v[i]/volt, label=f'Neuron {i}')
xlabel('Time (ms)')
ylabel('Voltage (V)')
legend()
show()
How SpikeGeneratorGroup Works Internally
When you create a SpikeGeneratorGroup, several sophisticated things happen behind the scenes:
Spike Storage and Sorting: Brian2 takes your spike times and indices and sorts them chronologically. This is crucial for efficient playback during simulation.
Time Binning: The system converts your continuous spike times into discrete time bins that match the simulation timestep. This involves careful rounding to ensure spikes occur at the right simulation steps.
Efficient Playback: During simulation, the group uses a pointer system to efficiently deliver spikes without searching through the entire spike list each timestep.
Looking at the internal structure from the source code:
# Key internal variables that SpikeGeneratorGroup creates:
self.variables.add_dynamic_array('neuron_index', values=indices, ...) # Which neurons
self.variables.add_dynamic_array('spike_time', values=times, ...) # When they spike
self.variables.add_dynamic_array('_timebins', ...) # Discretized times
self.variables.add_array('_lastindex', ...) # Current position
2. PoissonGroup: The Random Noise Generator
PoissonGroup generates spikes randomly according to Poisson statistics. Think of it like having many independent random spike generators, each with its own controllable firing rate.
# Create neurons that fire randomly at different rates
N_poisson = 10
poisson_neurons = PoissonGroup(N_poisson, rates=np.linspace(1, 20, N_poisson)*Hz)
# Create target network
target_network = NeuronGroup(5, '''
dv/dt = (-v + I_syn)/(10*ms) : volt
dI_syn/dt = -I_syn/(5*ms) : volt
''', threshold='v > 1*volt', reset='v = 0*volt')
# Connect Poisson input to network
poisson_synapses = Synapses(poisson_neurons, target_network,
'w : volt', on_pre='I_syn += w')
poisson_synapses.connect(p=0.3) # Random connectivity
poisson_synapses.w = '0.1*volt + 0.05*volt*rand()'
# Monitor everything
poisson_mon = SpikeMonitor(poisson_neurons)
target_mon = SpikeMonitor(target_network)
current_mon = StateMonitor(target_network, 'I_syn', record=True)
run(500*ms)
PoissonGroup's Mathematical Foundation
The beauty of PoissonGroup lies in its mathematical properties. A Poisson process has these key characteristics:
Independence: Each spike is independent of all others
Exponential ISIs: Time between spikes follows exponential distribution
Rate Parameter: Higher rates mean more frequent spikes
Memoryless: Past activity doesn't influence future spikes
3. PoissonInput: The Efficient Background Noise
PoissonInput is like having thousands of weak, independent Poisson inputs efficiently bundled together. Instead of creating individual synapses for each input (which would be computationally expensive), PoissonInput uses statistical sampling to achieve the same effect.
How PoissonInput Achieves Efficiency
The genius of PoissonInput lies in its use of the BinomialFunction under the hood. Instead of simulating 1000 individual synapses, it asks: "Given 1000 inputs each with probability p*dt of firing this timestep, how many actually fire?" This is a binomial distribution, which can be sampled much more efficiently.
4. BinomialFunction: The Mathematical Engine
BinomialFunction is the mathematical powerhouse that makes PoissonInput efficient. It's like having a sophisticated random number generator that can efficiently sample from binomial distributions.
5. TimedArray: The Signal Generator
TimedArray is like having a programmable signal generator that can produce any time-varying waveform you can imagine. It's perfect for creating complex, predetermined input patterns.
Namespaces
The Fundamental Problem: From Strings to Executable Code
Before we dive into Brian's solution, let's understand the core challenge. When you write a neural equation in Brian, you're writing it as a string:
# This is a string, not executable Python code
eqs = '''
dv/dt = (-v + I*R)/tau : volt
I = I_base + I_noise*xi : amp
'''
But for simulation, Brian needs to convert this into fast, executable code (Python, C++, or other languages). The challenge is: what do I_base
, I_noise
, R
, tau
, and xi
refer to? How does Brian know that tau
means your time constant variable, xi
is Gaussian white noise, and exp
is the exponential function?
This is where namespaces come in. Think of a namespace as Brian's "address book" that maps names to actual objects.
Understanding Namespaces: The Address Book Metaphor
A namespace is fundamentally a dictionary that maps names (strings) to Python objects:
# A simple namespace might look like this
namespace = {
'tau': 10*ms, # Your time constant
'I_base': 100*pA, # Your base current
'exp': numpy.exp, # The exponential function
'ms': brian2.ms, # The millisecond unit
'xi': '<special noise function>' # Gaussian white noise
}
When Brian encounters the name tau
in your equation, it looks it up in this address book to find the actual value 10*ms
. This process is called "namespace resolution."
The Namespace Hierarchy: Who Gets Priority?
Brian uses a sophisticated hierarchy of namespaces, like a series of address books that it checks in order. Understanding this hierarchy is crucial for predicting how Brian will interpret your equations.
The Resolution Order
When Brian encounters a name in an equation, it searches through namespaces in this order:
Group-specific namespace (highest priority)
Run namespace (passed to run() function)
Local namespace (variables in your current function/scope)
Global namespace (variables in your module/script)
Default namespaces (lowest priority)
Brian's built-in functions (exp, sin, cos, etc.)
Brian's units (ms, mV, nS, etc.)
Brian's constants (pi, e, inf, etc.)
Think of this as asking different people for a phone number: first you check your personal contacts, then your work contacts, then the phone book, and finally directory assistance.
Now let's examine the technical implementation of how Brian actually gathers and manages namespaces. This happens through sophisticated Python introspection.
Stack Frame Inspection: Reading Your Code's Context
The most fascinating part of Brian's namespace system is how it automatically captures your local variables. It does this through "stack frame inspection" - essentially looking at the state of your Python program when you create Brian objects:
def get_local_namespace(level):
"""Get the surrounding namespace by inspecting the call stack"""
# Get the locals and globals from the stack frame
frame = inspect.currentframe()
# Go back 'level' frames in the call stack
for _ in range(level + 1):
frame = frame.f_back
# Combine local and global variables from that frame
return dict(itertools.chain(frame.f_globals.items(), frame.f_locals.items()))
This function performs a kind of "time travel" - it looks back up the chain of function calls to see what variables were available when you created your Brian objects.
Functions in Brian2
The Fundamental Challenge: One Function, Many Languages
Before diving into the technical details, let's understand the core problem that Brian's function system solves. When you write a neural equation like this:
equations = '''
dv/dt = (-v + I*exp(-t/tau))/C : volt
firing_rate = clip(v/mV, 0, 100)*Hz : Hz
noise = 10*randn() : 1
'''
You're using several mathematical functions: exp()
, clip()
, and randn()
. But here's the challenge - Brian needs to make these work in multiple execution environments:
Python/NumPy runtime: For interactive development and small simulations
C++ standalone: For high-performance compiled simulations
Other targets: GPU, specialized hardware, etc.
Each environment has different syntax, different function names, and different capabilities. The exp()
function might be math.exp()
in Python, exp()
in C++, and expf()
for single precision. Brian's function system provides a unified interface that automatically adapts to each target.
Think of it like having a single remote control that works with any TV brand - the interface stays the same, but the underlying implementation adapts to each specific device.
The Function Class: Universal Function Abstraction
The Function
class is the cornerstone of Brian's function system. It represents a mathematical function in an abstract way, storing everything needed to use that function across different execution environments.
Anatomy of a Function Object
Let's examine how a Function
object is constructed:
class Function:
def __init__(self, pyfunc, sympy_func=None, arg_units=None,
return_unit=None, arg_types=None, return_type=None,
stateless=True, auto_vectorise=False):
self.pyfunc = pyfunc # Python implementation
self.sympy_func = sympy_func # Symbolic mathematics representation
self._arg_units = arg_units # Expected units for arguments
self._return_unit = return_unit # Units of return value
self._arg_types = arg_types # Data types of arguments
self._return_type = return_type # Data type of return value
self.stateless = stateless # Whether function has internal state
self.auto_vectorise = auto_vectorise # Whether to add vectorization info
Each Function
object knows several crucial things about the mathematical function it represents:
Python Implementation: The actual Python function that computes the result
Symbolic Representation: How symbolic math systems should interpret this function
Unit Information: What units the arguments should have and what units the result will have
Type Information: What data types (integer, float, boolean) are expected
Behavioral Properties: Whether the function has side effects or needs special handling
Structure of DEFAULT_FUNCTIONS
# This is how Brian structures its default function library
DEFAULT_FUNCTIONS = {
# Trigonometric functions
'sin': Function(unitsafe.sin, sympy_func=sympy.sin),
'cos': Function(unitsafe.cos, sympy_func=sympy.cos),
'tan': Function(unitsafe.tan, sympy_func=sympy.tan),
# Exponential and logarithmic functions
'exp': Function(unitsafe.exp, sympy_func=sympy.exp),
'log': Function(unitsafe.log, sympy_func=sympy.log),
'log10': Function(unitsafe.log10, sympy_func=sympy.log10),
# Power and root functions
'sqrt': Function(np.sqrt, sympy_func=sympy.sqrt,
arg_units=[None], return_unit=lambda u: u**0.5),
# Random number functions
'rand': Function(pyfunc=rand, arg_units=[], return_unit=1,
stateless=False, auto_vectorise=True),
'randn': Function(pyfunc=randn, arg_units=[], return_unit=1,
stateless=False, auto_vectorise=True),
# Utility functions
'clip': Function(pyfunc=np.clip, arg_units=[None, "a", "a"],
arg_names=["a", "a_min", "a_max"],
return_type="highest", return_unit=lambda u1, u2, u3: u1),
}
Why Preference Systems Matter
Before diving into the technical details, let's understand why Brian needs such a sophisticated preference system. Imagine Brian as a highly customizable sports car - you might want different settings for different situations:
Performance mode: Fast floating-point operations, aggressive optimizations
Debug mode: Slower but more careful error checking
Research mode: Different numerical precision for reproducibility
Teaching mode: Simpler defaults for educational use
Brian's preference system manages hundreds of such settings across the entire framework, ensuring they're consistent, validated, and easily accessible.
BrianPreference: The Building Block of Configuration
Let's start with BrianPreference
, which is like a smart container for a single setting. Think of it as a sophisticated variable that knows not just its value, but also its documentation, validation rules, and how to represent itself.
Anatomy of a BrianPreference
class BrianPreference:
def __init__(self, default, docs, validator=None, representor=repr):
self.default = default # The default value
self.docs = docs # Human-readable documentation
self.validator = validator # Function to check if values are valid
self.representor = representor # Function to convert values to strings
Each BrianPreference
is like a specialized container that knows four crucial things:
1. Default Value: What should this setting be if the user doesn't specify anything?
2. Documentation: What does this setting do and how should it be used?
3. Validation: What values are acceptable for this setting?
4. Representation: How should this setting be displayed or saved to a file?
# This is how Brian defines the default float type preference
default_float_dtype = BrianPreference(
default=float64, # Use 64-bit floats by default
docs="""
Default dtype for all arrays of scalars (state variables, weights, etc.).
This affects memory usage and numerical precision throughout Brian.
""",
representor=dtype_repr, # Custom function to display numpy dtypes nicely
validator=default_float_dtype_validator # Ensures only float32/float64 allowed
)
This preference encapsulates a crucial decision: should Brian use 32-bit floats (faster, less memory, less precise) or 64-bit floats (slower, more memory, more precise)? The preference system makes this choice explicit and modifiable.
The Validation System: Keeping Things Safe
The validator is like a bouncer at a club - it decides what gets in and what gets rejected. Here's how Brian implements the float type validator:
def default_float_dtype_validator(dtype):
"""Only allow float32 or float64"""
return dtype in [float32, float64]
If you try to set the float type to something invalid, the validator catches it:
# This would work
prefs.core.default_float_dtype = np.float32 ✓
# This would raise PreferenceError
prefs.core.default_float_dtype = np.int32 ✗ - PreferenceError!
The Default Validator: Smart Type Checking
When you don't provide a custom validator, Brian uses the DefaultValidator
, which is surprisingly intelligent:
class DefaultValidator:
def __init__(self, value):
self.value = value # Remember the default value
def __call__(self, value):
# Must be the same type as default
if not isinstance(value, self.value.__class__):
return False
# For Quantity objects, units must match
if isinstance(self.value, Quantity):
if not have_same_dimensions(self.value, value):
return False
return True
This validator automatically ensures type safety and unit consistency. If the default is 10*ms
, you can only set values with time dimensions - no accidentally setting a voltage where time is expected!
The Global prefs Object: Mission Control
Now let's explore prefs
, which is Brian's global preference manager. Think of it as the master control panel for the entire Brian framework - every preference in the system flows through this single object.
BrianGlobalPreferences: The Orchestrator
class BrianGlobalPreferences(MutableMapping):
def __init__(self):
self.prefs = {} # Actual preference values
self.backup_prefs = {} # For restore functionality
self.prefs_unvalidated = {} # Preferences waiting for validation
self.pref_register = {} # Registry of all possible preferences
self.eval_namespace = {} # Namespace for evaluating string values
This class is fascinating because it acts like a dictionary but with superpowers. You can access preferences in multiple ways:
python# Dictionary-style access
prefs['core.default_float_dtype'] = np.float32
# Attribute-style access (more Pythonic)
prefs.core.default_float_dtype = np.float32
# Both do exactly the same thing!
Preference Registration: The Foundation
Before any preference can be used, it must be registered. This is like declaring all possible settings before the system starts:
def register_preferences(self, prefbasename, prefbasedoc, **prefs):
"""Register a category of preferences with their definitions"""
# Example call:
# prefs.register_preferences(
# 'core', # Category name
# 'Core Brian preferences', # Category documentation
# default_float_dtype=BrianPreference(...),
# default_integer_dtype=BrianPreference(...),
# )
self.pref_register[prefbasename] = (prefs, prefbasedoc)
# Set all preferences to their default values
for k, v in prefs.items():
fullname = prefbasename + '.' + k
self.prefs_unvalidated[fullname] = v.default
Registration Example
Here's how Brian registers its core preferences:
prefs.register_preferences(
'core',
'Core Brian preferences',
default_float_dtype=BrianPreference(
default=float64,
docs='Default dtype for all arrays of scalars.',
representor=dtype_repr,
validator=default_float_dtype_validator,
),
default_integer_dtype=BrianPreference(
default=int32,
docs='Default dtype for all arrays of integer scalars.',
representor=dtype_repr,
),
outdated_dependency_error=BrianPreference(
default=True,
docs='Whether to raise an error for outdated dependencies.',
),
)
This single call creates multiple preferences under the 'core' category, each with its own validation and documentation.
Subscribe to my newsletter
Read articles from Mrigesh Thakur directly inside your inbox. Subscribe to the newsletter, and don't miss out.
Written by
