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, or amp—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:

1. Checks that its input is dimensionless (like an angle in radians)

2. If so, calls the regular np.sin

3. 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:

1. Variables in Brian2

2. 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:

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 conditions

• Intermediate 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:

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 → Integrating dv/dt

• threshold.pyx → Detecting spikes

• reset.pyx → Resetting post-spike neurons

• synapses.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:

1. Devices – the where and how of execution

2. 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:

1. Film it live on stage - Everything happens in real-time, you see results immediately, but it's limited by the theater's resources

2. 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:

1. Immediate execution: The simulation loop starts right away

2. Direct array access: Code objects read/write directly to numpy arrays

3. Python/Cython speed: Each timestep executes compiled code, but coordination happens in Python

4. 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:

1. Your equation: dv/dt = -v/tau

2. Brian's analysis: This becomes vector code that updates each neuron

3. Code generation: The template gets filled with C++ code

4. 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:

1. Contains simulation objects: Neuron groups, synapse groups, monitors, etc.

2. Controls simulation flow: Determines when different components should update

3. Manages simulation time: Keeps track of the current simulation time and advances it

4. Handles scheduling: Ensures objects are updated in the correct order

5. 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 called

• Clocks are set up to cover the simulation duration

Step 2: Main Simulation Loop

The main simulation loop does the following at each step:

1. Find the clock(s) with the smallest current time

2. Set the network time to this clock's time

3. Update all objects associated with these clocks in their scheduled order

4. Advance these clocks by their time step (dt)

5. 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 called

• Profiling 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:

1. Identifying which clocks have the earliest time

2. Updating only objects associated with those clocks

3. Advancing those clocks, leaving other clocks unchanged

4. 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, and PopulationRateMonitor 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

1. 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:

1. Connection to Source: The monitor establishes a connection to your neuron group and specifically watches for the 'spike' event.

2. Variable References: It creates references to any variables you want to record (like 'v' and 'I' in our example).

3. CodeRunner Creation: The monitor becomes a CodeRunner that executes code every time the watched event occurs.

4. 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 neurons

• record=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:

1. Runs at specific times: Monitors typically run in the 'start' slot for StateMonitor or immediately after spike detection for spike monitors.

2. Has access to group variables: Through the Variables system, monitors can read any variable from the groups they're monitoring.

3. 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

2. 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:

1. Spike Storage and Sorting: Brian2 takes your spike times and indices and sorts them chronologically. This is crucial for efficient playback during simulation.

2. 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.

3. 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:

1. Independence: Each spike is independent of all others

2. Exponential ISIs: Time between spikes follows exponential distribution

3. Rate Parameter: Higher rates mean more frequent spikes

4. 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.

3. 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:

1. Group-specific namespace (highest priority)

2. Run namespace (passed to run() function)

3. Local namespace (variables in your current function/scope)

4. Global namespace (variables in your module/script)

5. 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.

4. 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),
}

5. 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.