Devlog #83 – Wave 62: Neuron Network, Cell Growth, Gyroid, Generative City & Computational Photo

New Simulations

🧠 Neuron Network — Leaky Integrate-and-Fire Dynamics

The LIF model

The Leaky Integrate-and-Fire (LIF) neuron is the workhorse of computational neuroscience. Despite its simplicity, it reproduces the essential features of real neurons: integration of synaptic inputs, a threshold-crossing spike, and subthreshold voltage decay back to rest. The membrane potential V of each neuron obeys:

τ_m · dV/dt = -(V - V_rest) + R · I(t)

where:
  τ_m   = membrane time constant (typically 10–20 ms)
  V_rest = resting potential (−65 mV)
  R     = membrane resistance
  I(t)  = total synaptic current at time t

Spike rule:
  if V ≥ V_thresh:
      record spike; V ← V_reset; enter refractory for t_ref ms

The simulation discretises this with forward Euler at dt = 1 ms — coarse enough to run in real time but fine enough to capture inter-spike intervals accurately. With τ_m = 10 ms, the free decay per step is exp(−dt/τ_m) ≈ 0.905, so a neuron receiving no input loses about 9.5% of its excess voltage each millisecond.

Synaptic model and E/I balance

Each spike from neuron j delivers a current pulse to neuron i after a fixed synaptic delay d_ij (1–5 ms). The current decays exponentially with synaptic time constant τ_syn:

I_syn(t) = w_ij · exp(-(t - t_spike - d_ij) / τ_syn)   [for t > t_spike + d_ij]

Excitatory synapse: w_ij > 0  (AMPA-like, τ_syn ≈ 5 ms)
Inhibitory synapse: w_ij < 0  (GABA-like, τ_syn ≈ 10 ms)

E/I ratio: 80% excitatory neurons, 20% inhibitory (Dale's law)

The E/I balance is crucial. When excitation dominates, the network enters a run-away firing state (epileptiform bursts). When inhibition dominates, activity quenches. The balanced preset (roughly equal total excitatory and inhibitory drive) produces the asynchronous irregular (AI) state that resembles cortex during active wakefulness.

Voltage heatmap and spike raster

The left panel renders a 2D heatmap where each pixel corresponds to a neuron (laid out in a grid), coloured by membrane voltage: dark blue at V_rest, transitioning through cyan and yellow to bright white at threshold. The right panel shows a scrolling spike raster — each dot is one spike from one neuron, plotted at its time of occurrence. The raster allows immediate visual diagnosis of network state: sparse random dots indicate AI state, horizontal bands of synchronous dots indicate oscillatory or burst activity.

Synchrony and frequency statistics

Two live statistics are computed each frame. Population firing rate r(t) is the number of spikes in the last 50 ms divided by N × 0.050 s, giving a rate in Hz. Synchrony is measured as the coefficient of variation of the inter-spike interval (CV_ISI) averaged across neurons: CV = σ_ISI / μ_ISI. CV ≈ 1 indicates Poisson-like (AI) spiking; CV ≪ 1 is regular (clock-like); CV ≫ 1 indicates bursting. The oscillatory preset drives the network into a ~40 Hz gamma rhythm, visible as rhythmic columns on the raster and CV ≈ 0.3.

🦠 Cell Growth & Morphogenesis — Division, Packing and Pattern Formation

Force-based cell packing

Cells are modelled as soft disks with position x_i and radius r_i. At each timestep, three forces act on every cell:

F_push(i,j)     = k_push  · max(0, r_i + r_j - |x_i - x_j|) · n̂_ij   // overlap repulsion
F_adhesion(i,j) = k_adh   · max(0, |x_i - x_j| - (r_i+r_j+d_adh)) · n̂_ij  // short-range attraction
F_boundary(i)   = k_wall  · max(0, r_i - dist_to_boundary) · n̂_wall   // soft container wall

Net: x_i += (ΣF) · dt / ζ    (overdamped; no inertia)

The spatial grid accelerates the neighbour search to O(n) by hashing cells into grid buckets of size 2 r_max. Only cells in the same or adjacent buckets need force evaluation — this keeps the simulation at 60 fps even near the 600-cell cap.

Gray-Scott reaction-diffusion

In parallel, an 80×60 continuous grid evolves under the Gray-Scott equations for two chemical species U and V:

∂U/∂t = D_U · ∇²U  -  U·V²  +  f·(1 - U)
∂V/∂t = D_V · ∇²V  +  U·V²  -  (f + k)·V

Typical parameters for spots pattern:
  D_U = 0.16,  D_V = 0.08,  f = 0.035,  k = 0.065

Laplacian (5-point stencil, periodic boundaries):
  ∇²U[i,j] = U[i-1,j] + U[i+1,j] + U[i,j-1] + U[i,j+1] - 4·U[i,j]

The feed rate f and kill rate k determine which Turing pattern emerges: spots, stripes, labyrinthine mazes or spatiotemporal chaos. The simulation maps cell positions onto the continuous grid to sample the local V concentration, which drives cell fate: high V → specialised (differentiated) fate; low V → stem-like.

Cell division and colour modes

A cell grows at a constant rate dr/dt = growth_rate. When r > r_divide, it splits into two daughter cells at ±45° from a random orientation, each with r = r_divide / √2 (conserving area). The division event generates a short-range pressure burst that propagates through the colony via the force system, producing the characteristic cell-packing dynamics seen in growing epithelia.

• Fate mode — Colour encodes Gray-Scott V concentration sampled at the cell's grid position. Morphogen gradients produce spatially coherent colour domains.
• Age mode — Colour encodes cell generation number (time since last division), from young (cyan) to old (magenta).
• Size mode — Colour encodes current radius relative to division threshold, showing which cells are about to divide (bright white).

💠 Gyroid — Triply Periodic Minimal Surfaces

What is a TPMS?

A Triply Periodic Minimal Surface (TPMS) is a surface that is periodic in all three spatial directions and has zero mean curvature at every point (H = 0, the defining property of a minimal surface). Nature produces these structures spontaneously: they appear in butterfly wing scales, block copolymer membranes, and lipid bilayer systems. The gyroid in particular is a bicontinuous structure — both sides of the surface form interconnected labyrinths — making it highly prized in materials science for filtration, energy storage, and photonic applications.

Implicit surface equations

Each surface is defined by an implicit equation F(x, y, z) = t where t is an isovalue (default t = 0 gives the exact minimal surface). In scaled coordinates (u = 2πx/a, v = 2πy/a, w = 2πz/a):

Gyroid:   sin(u)cos(v) + sin(v)cos(w) + sin(w)cos(u) = t
Schwarz P: cos(u) + cos(v) + cos(w) = t
Diamond:  sin(u)sin(v)sin(w) + sin(u)cos(v)cos(w)
        + cos(u)sin(v)cos(w) + cos(u)cos(v)sin(w) = t
Neovius:  3(cos(u)+cos(v)+cos(w)) + 4cos(u)cos(v)cos(w) = t

The gyroid was discovered by Alan Schoen in 1970 during a NASA study of lightweight structures. Unlike the Schwarz P and D surfaces, it has no straight lines or planes of mirror symmetry, giving it a chiral, labyrinthine character.

Marching cubes algorithm

The simulation evaluates F on a regular grid of N³ voxels (default N = 48) and extracts the isosurface F = t using the Bourke 256-case marching cubes lookup table. Each voxel cube is assigned an 8-bit index (one bit per corner, set if F > t) — the index looks up one of 256 triangle configurations, giving up to 5 triangles per cube. Vertices are placed on cube edges by linear interpolation:

// For edge from vertex A (value fA) to vertex B (value fB):
t_lerp = (isovalue - fA) / (fB - fA)
vertex  = A + t_lerp · (B - A)

Total triangles: O(N² · surface area in voxels)
For N=48, tiling=1: typically ~80,000–120,000 triangles

Vertex normals are computed analytically from ∇F rather than by averaging face normals — this gives perfectly smooth shading at any resolution. The mesh is uploaded to a Three.js BufferGeometry and rendered with a MeshStandardMaterial under an environment map, producing the metallic/iridescent look.

Tiling and isovalue controls

The tiling slider repeats the unit cell 1–4 times in each axis, demonstrating the periodic nature of the surface. Increasing the isovalue t shifts the surface toward one side of the bicontinuous split: at t = 0 the two labyrinths are equal volume; at t > 0 one labyrinth shrinks and eventually pinches off, transitioning from a TPMS to a closed-cell foam topology.

🌃 Generative City — Procedural Urban Landscape

Recursive block subdivision

The city layout is generated by recursively subdividing a rectangular canvas. Starting from the full canvas rectangle, the algorithm splits it along either the horizontal or vertical axis, adding a small random jitter (±10% of side length) to the split position to avoid a grid-like feel. Splitting continues until the block width or height falls below a minimum lot size. The split axis alternates deterministically based on the current rectangle's aspect ratio: always split the long side.

function subdivide(rect, depth):
    if rect.w < MIN_LOT or rect.h < MIN_LOT: return [rect]
    axis = rect.w > rect.h ? 'x' : 'y'
    split = lerp(0.35, 0.65, lcg_random()) * (axis='x' ? rect.w : rect.h)
    [left, right] = split_rect(rect, axis, split + jitter())
    return subdivide(left, depth+1) + subdivide(right, depth+1)

This produces an irregular but realistic block pattern that resembles organic city growth rather than a rigid Cartesian grid.

Zoning and density gradient

Each leaf block is assigned a zone based on its distance from the city centre (normalised to [0, 1]):

• Commercial (d < 0.25) — Dense tall buildings; rectangles fill most of the block, leaving only a narrow street gap. Colours: warm greys and slate blues.
• Residential (0.25 ≤ d < 0.65) — Smaller buildings with gardens. Blocks contain 2–8 house footprints with a shared courtyard gap. Colours: warm tans.
• Park (random 12% of mid-ring blocks) — Filled with a stippled green pattern. Provides visual breathing room and realistic open space.
• Water (contiguous blobs in outer ring, ~8% of area) — Rendered with animated shimmer lines. Rivers and lakes emerge from merged adjacent water blocks.

Building height is encoded by colour brightness: taller buildings receive lighter fill and a white cap stripe, giving the top-down view an implicit 3D cue.

Seeded LCG random for reproducibility

All randomness flows through a Linear Congruential Generator seeded from the user's seed value:

// LCG parameters (Numerical Recipes)
a = 1664525;  c = 1013904223;  m = 2^32
seed = (a * seed + c) % m
random() { seed = ...; return seed / m; }

The same seed always produces the same city. The four themes (night, day, cyberpunk, sunset) change only the colour palette, not the layout — so you can compare how the same urban structure looks at different times of day.

🖼️ Computational Photography — Four Image Processing Algorithms

ASCII art via luminance mapping

The image is downsampled to a character grid (default 120 columns). Each cell's average luminance L = 0.299R + 0.587G + 0.114B (ITU-R BT.601) maps to one of 70 characters sorted by visual density from darkest to lightest:

charset = "$@B%8&WM#*oahkbdpqwmZO0QLCJUYXzcvunxrjft/\\|()1{}[]?-_+~<>i!lI;:,\"^`'. "
char_index = floor(L / 255 * (charset.length - 1))
// Render in a monospaced font; colour from original pixel for colour-ASCII mode

Photomosaic via tile colour averaging

The photomosaic effect divides the image into rectangular tiles (size adjustable from 4×4 to 32×32 pixels) and replaces each tile with a solid fill of its average colour. This is the digital equivalent of a classical photo mosaic, and at large tile sizes it functions as a pixelation / privacy filter. The averaging uses a full-scan of all pixels in each tile — not a sampled approximation — ensuring colour accuracy.

Poisson-disk stippling

Stippling converts an image into a set of dots whose density encodes brightness. The simulation uses a rejection-sampling approximation: candidate dot positions are drawn uniformly at random; a dot is placed only if the local luminance (inverted, so dark areas get more dots) exceeds a threshold. The minimum inter-dot distance enforces a Poisson-disk distribution, preventing clumping:

for each candidate (x, y):
    L = luminance(x, y)
    density = 1 - L/255           // dark pixel → high density
    if random() < density and no_dot_within_r(x, y, r_min):
        place dot at (x, y)

Floyd-Steinberg error diffusion dithering

Classic Floyd-Steinberg dithering produces a 1-bit (black/white) representation that perceptually matches the continuous original. The algorithm scans pixels left-to-right, top-to-bottom. Each pixel is quantised to black (0) or white (255); the quantisation error is distributed to four neighbours:

error = original_pixel - quantised_pixel

pixel[x+1][y  ] += error × 7/16
pixel[x-1][y+1] += error × 3/16
pixel[x  ][y+1] += error × 5/16
pixel[x+1][y+1] += error × 1/16

The error propagates forward and downward, so dark regions accumulate positive error (making later pixels more likely to round to white) and vice versa. The result is a halftone-like pattern with a high-frequency dither noise that, at normal viewing distance, appears as a grey tone matching the original luminance. Floyd-Steinberg dithering was published in 1976 and remains the standard algorithm in image editors and printers worldwide.

Up Next

Wave 63 moves into fluid dynamics, combinatorial mathematics and optics. Coming up: Stokes creeping flow with a reversibility demonstration, 2D and 3D sphere packing fractions, space-filling Hilbert and Peano curves, the Ulam spiral of primes, Edmonds-Karp max-flow with animated augmenting paths, and a Monte Carlo ice-halo ray tracer reproducing 22° halos and sun dogs.