Biological Brain Simulator

An open-source, research-oriented simulator for studying biologically inspired learning in spiking neural systems. The project combines spiking neurons, synaptic plasticity, structural growth, memory / replay, oscillations, and full-state persistence in a single experimental framework.

The simulator is an experimental sandbox for testing whether local learning rules and biologically motivated dynamics can support nontrivial behavior on controlled AI tasks. It does not claim to model a whole brain faithfully.

Modeling scope

Biological mechanism	Implementation
Neurons	Izhikevich model — reproduces 7 real firing patterns (regular spiking, bursting, chattering, fast spiking, ...)
Synapses	Variable-weight connections with axonal delay, neurotransmitters (glutamate, GABA, dopamine, serotonin, acetylcholine), short-term plasticity (vesicle depletion + facilitation)
STDP	Event-driven Spike-Timing Dependent Plasticity — neurons that fire together wire together, computed only at spike events
R-STDP	Reward-Modulated STDP — three-factor learning rule (STDP × dopamine) with per-synapse eligibility traces and per-target dopamine
Homeostatic plasticity	Synaptic scaling to keep firing rates in a healthy range
Metaplasticity	The plasticity of plasticity itself adapts (BCM theory)
Neurogenesis	Birth of new neurons in active regions
Synaptogenesis	New connections between co-active neurons
Pruning	Elimination of weak/unused synapses
Apoptosis	Programmed death of chronically inactive neurons
Memory	Activity traces, consolidation via replay (like during sleep), pattern completion
Oscillations	Theta/alpha/beta/gamma rhythms with theta-gamma cross-frequency coupling
Morphology	Multi-compartment neurons (dendrites/soma/axon) with distance-based synaptic attenuation
Persistence	Full save/restore of the entire brain state to disk

Morphology templates assign synapses to dendritic compartments and apply a static attenuation factor based on distance from the soma. Compartments do not maintain independent voltage states; Izhikevich membrane dynamics are integrated at the soma.

Research positioning

Mainstream ML pipeline	This simulator
Distinct training and inference phases	Continuous adaptation during interaction
Fixed size	Grows organically (neurogenesis)
Backpropagation	Local plasticity rules (STDP, R-STDP, homeostasis)
Abstract activations	Spiking membrane dynamics
Minimal internal structure	Regions, projections, neurotransmitters, oscillations

Status

This is an experimental research codebase, not a polished general-purpose framework. It supports controlled benchmarks and experiments with biologically inspired learning mechanisms. The benchmark scripts are the canonical reference for exact experimental settings.

The repository contains supervised R-STDP, reward-modulated navigation, and unsupervised MNIST protocols. Current validation results and their limits are documented in BENCHMARKS.md. Historical scores should not be compared directly with runs made after the protocol and dynamics corrections.

Quickstart

Requirements: Python >= 3.10, plus PyTorch, matplotlib, networkx, and scikit-learn as listed in requirements.txt.

pip install -r requirements.txt

# Run from the project root
python examples/association_demo.py
python examples/growth_demo.py
python examples/iris_benchmark.py      # classification benchmark
python examples/grid_nav_benchmark.py  # grid navigation benchmark
python examples/mnist_benchmark.py     # MNIST benchmark (10-class, 784+400+400)
python examples/mnist_diagnosis.py     # short MNIST diagnosis sweeps

Project structure

src/
├── __init__.py      # Public API re-exports
├── neuron.py        # Neuron types and Izhikevich parameters
├── integration.py   # Stable Heun dynamics, CPU loop and legacy Euler compatibility
├── synapse.py       # Neurotransmitter types and properties
├── device.py        # Compute device selection (CPU/MPS/CUDA)
├── region.py        # Brain region (vectorized PyTorch tensors)
├── brain.py         # Top-level orchestrator + inter-region projections
├── plasticity.py    # STDP, R-STDP, homeostatic, metaplasticity (BCM)
├── growth.py        # Neurogenesis, synaptogenesis, pruning, apoptosis
├── memory.py        # Memory traces, consolidation, pattern completion
├── stimulus.py      # Sensory encoding (rate/temporal/population coding)
├── oscillator.py    # Brain oscillations (theta/alpha/beta/gamma)
├── morphology.py    # Multi-compartment neuron morphology
└── persistence.py   # Full brain serialization and deserialization

examples/
├── association_demo.py       # Pavlovian conditioning with spiking neurons
├── growth_demo.py            # A brain that grows from scratch
├── iris_benchmark.py         # Iris classification with R-STDP
├── grid_nav_benchmark.py     # Grid navigation with R-STDP
├── mnist_benchmark.py        # MNIST benchmark: 10-class unsupervised STDP
├── mnist_diagnosis.py        # MNIST diagnosis: readout / inhibition / STDP / theta sweeps
├── fewshot_benchmark.py      # Few-shot learning: SNN vs MLP data efficiency
├── forgetting_benchmark.py   # Catastrophic forgetting: sequential task retention
├── degradation_benchmark.py  # Graceful degradation: robustness to neuron damage
├── sleep_benchmark.py        # Sleep/replay: consolidation benefit measurement
└── growth_benchmark.py       # Growth vs fixed: structural plasticity ablation

BENCHMARKS.md              # Concise benchmark notes, results, and caveats

Usage

from src.brain import Brain
from src.region import RegionType
from src.persistence import save_brain, load_brain

brain = Brain(dt=1.0, seed=42)

# Available region types: SENSORY, ASSOCIATION, MEMORY, MOTOR, REWARD
brain.add_region("input", RegionType.SENSORY, n_neurons=50)
brain.add_region("cortex", RegionType.ASSOCIATION, n_neurons=200)
brain.add_region("output", RegionType.MOTOR, n_neurons=10)
brain.connect_regions("input", "cortex", density=0.1)
brain.connect_regions("cortex", "output", density=0.1)

# Enable multi-compartment neuron morphology (optional)
brain.regions["cortex"].enable_morphology()

# Stimulate and let it learn
for step in range(10000):
    brain.stimulate("input", [0.8, 0.3, 0.5])
    if step == 5000:
        brain.reward(2.0, target="proj:input->cortex")
    if step == 7000:
        brain.punish(1.0, target="proj:cortex->output")
    brain.step()

# The brain has grown
print(brain.summary())

# Save to disk — full state, down to individual synapse weights
save_brain(brain, "my_brain")

# Load and continue from where it left off
brain = load_brain("my_brain")
for step in range(5000):
    brain.stimulate("input", [0.2, 0.9, 0.4])
    brain.step()

Convenience helpers

The core simulator provides helpers for benchmarks and structured training loops:

• brain.reset_traces() clears all synaptic eligibility traces and per-target dopamine
• brain.freeze_structural_plasticity() disables growth and metaplasticity for stable experiments
• brain.freeze_homeostatic_scaling() freezes synaptic scaling while leaving adaptive thresholds active
• brain.freeze_adaptive_thresholds() freezes adaptive thresholds independently
• brain.freeze_plasticity() freezes STDP/R-STDP on existing regions and projections, clears pending eligibility/dopamine, and keeps metaplasticity from re-enabling them; synaptic scaling and growth have separate controls
• brain.enable_projection_plasticity("input", "output") re-enables STDP on one projection
• brain.disable_memory(), brain.disable_oscillations(), and brain.disable_reward_modulated_plasticity() isolate experimental subsystems without monkey-patching methods
• brain.get_projection("input", "output") fetches an inter-region projection safely
• brain.regions["output"].add_lateral_inhibition(weight=5.0) creates all-to-all lateral inhibition

Example:

brain.freeze_plasticity()
brain.enable_projection_plasticity("input", "output", A_plus=0.01, A_minus=0.012)
brain.freeze_structural_plasticity()
brain.regions["output"].add_lateral_inhibition(weight=5.0)

To reinforce only one output, restrict credit for the entire dopamine pulse:

target = Brain.projection_target("input", "output")
with brain.reward_stdp.restrict_to_posts(target, [chosen_output]):
    brain.reward(0.1, target)
    for _ in range(10):
        brain.step()

The restriction applies to both existing eligibility and new spikes. Dopamine for that target is cleared on entry and exit; overlapping restrictions on the same target are rejected. Save checkpoints after leaving the context. Saves preserve plasticity flags and restore the previous checkpoint if a handled error or interruption prevents replacement. They do not guarantee atomic replacement across process crashes or power loss.

Experimental results

The repository includes three reference benchmark protocols. Their validation status differs:

Benchmark	Script	Main question	Current status
Iris classification	examples/iris_benchmark.py	Can reward-modulated local plasticity solve a standard supervised classification task?	Revalidated at 90% spike-only test accuracy (27/30), one seed; epoch selected on separate validation data
Grid navigation	examples/grid_nav_benchmark.py	Can the simulator learn a usable control policy with reward-modulated spiking dynamics?	Calibrated readout solves all 24 starts entirely through spikes on three seeds; mean 4.25 steps across seeds in the fixed 5×5 world
MNIST	examples/mnist_benchmark.py	Can the simulator scale to a non-trivial unsupervised vision benchmark?	Canonical MNIST split, train-only preprocessing, fixed WTA and adaptive theta; revalidation pending

A smaller MNIST study trained three networks on 1,000 images and fitted each ridge readout on 100 labelled images from that training set. Centering each image's voltage features across excitatory neurons raised mean accuracy from 61.54% to 63.04%, a gain of 1.50 percentage points. Both decoders used the same 800 additional validation images for every network, excluding the source studies' training and prior validation images. The canonical test set was not used in this comparison. See the voltage-centering results for per-seed scores and checks.

examples/mnist_learning_check.py reports the centered-voltage ridge alongside its three existing decoders. The gain comes from the readout on fixed networks; it does not establish better STDP learning. The main MNIST benchmark retains its template classifier. Its full 50,000-image training run and canonical test evaluation remain to be repeated with the corrected simulator.

A subsequent fixed-budget check varied 100-label readout subsets on three frozen networks and reused 800 validation images. The STDP condition averaged 62.54% with centered-voltage ridge and 36.02% with spike-only ridge; neither comparison showed a consistent STDP benefit over its matched control. These are validation sensitivity measurements, not a new canonical test score. See the fixed-budget results.

An isolated conductance-LIF reference pilot tests the triplet-STDP variant in Diehl and Cook's released code. Its completed paired experiment found a consistent STDP benefit in spike-based classification across three seeds. A subsequent train-only decoder selection improved the combined spike/voltage readout without retraining the network. A frozen evaluation on new images confirmed the decoder gain, while finding only a small, uncertain advantage over voltage-only decoding and the initial network. Separate frozen-inference and early-training diagnostics replay matched input events across timestep grids; weight/theta swaps probe their separate effects without retraining. The linked reports separate these effects from the numerical limits. This experimental implementation leaves Iris and the production simulator unchanged and does not reproduce the paper's 87% power-law result.

A deterministic timing audit identified late-inhibition effects and a decimal-grid refractory defect. The corrected reference v2 uses integer refractory deadlines and eight complete internal steps per 0.5 ms input interval. Existing checkpoints retain v1 dynamics. A matched 50-image training check found that halving v2's step still changes learned weights and neuron-level responses substantially; v2 used about 3.19 times v1's process CPU time in that short paired comparison. The reported MNIST scores, including 69.27% on fresh validation images, belong to v1. V2 accuracy has not yet been measured.

BENCHMARKS.md records the protocols, results and limitations. The benchmark scripts define the exact hyperparameters.

Computational performance

The simulator is implemented as a vectorized research codebase rather than a neuron-by-neuron object model. All neuron and synapse state lives in dense PyTorch tensors (structure-of-arrays layout), so a single Region.step() call can advance thousands of neurons in parallel.

Key techniques:

• PyTorch tensors on CPU: profiling through 3200 neurons and 2.9M synapses showed CPU outperforming MPS. GPU kernel launch overhead dominates the tested conditional operations, torch.compile falls back because of dynamic shapes, and masked full-tensor approaches process the roughly 95-99% of inactive synapses. See BENCHMARKS.md for measurements. Use BRAIN_DEVICE=mps to profile other network sizes or densities.
• Vectorized Izhikevich integration: new networks use second-order Heun integration with stability-limited internal steps, capped at 0.1ms by default. The external timestep still controls synapses, plasticity and spike timestamps. Version 4/5 checkpoints retain their legacy Euler dynamics; version 6 saves the method and step cap for each region. See BENCHMARKS.md for the numerical checks and remaining validation work.
• CPU Heun loop: compatible float32/float64 CPU arrays use NumPy operations in the same order as the PyTorch reference. GPU, autograd, scalar and mixed-dtype inputs retain the PyTorch loop. Whole-training state/RNG and inference responses are checked for exact equality; the measured speedups are relative to Heun, not the former Euler model.
• Ring buffer for spike delays: index_add_ deposits postsynaptic currents into future slots; each timestep reads and clears the current slot
• Event-driven STDP: nearest-neighbor pairing updates synapses whose pre or post neuron fired this step. Large, sparsely active CPU pathways use checked adjacency indices; other cases retain the PyTorch scan. Indexed updates preserve the original synapse and LTP/LTD order.
• Sparse scatter/gather propagation: synapses remain in COO-style parallel arrays. CPU adjacency lookup avoids gathering every connection's firing flag, while comparing endpoint arrays with a cached snapshot to detect topology edits. This validity check still scans all endpoints; transmission updates only the selected live synapses. Small pathways, dense activity and non-CPU devices use fired[syn_pre] filtering.
• Vectorized stimulus encoding: rate, temporal, and population coding run in a single tensor operation instead of per-neuron Python loops
• Vectorized synaptogenesis: co-activity scoring uses torch.outer instead of nested Python loops
• Pre-allocated neuron arrays up to max_neurons; synapse arrays grow via capacity-doubling when needed
• Dead flags (syn_alive, neuron_alive) for pruning and apoptosis without array compaction

Reusing responses in independent evaluation reduced median Iris evaluation time from 10.27 to 1.70 seconds (6.03x) and Grid evaluation time from 11.56 to 2.92 seconds (3.96x). The checks reproduced the reference outputs and preserved the source model and random number generator state. In a separate compact MNIST comparison, the CPU Heun loop was 1.44x faster for training and 1.56x faster for inference than the PyTorch Heun loop. These measurements depend on the workload and host load; their ratios do not estimate a combined speedup. BENCHMARKS.md documents the timing protocols and subsequent CPU optimizations.

License

This project is licensed under the MIT License.

Author

Francesco Pace
Email: francesco.pace@gmail.com LinkedIn: linkedin.com/in/francescopace