diff --git a/CHANGELOG.md b/CHANGELOG.md
index 03d1c1ed..e300afa3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -53,6 +53,8 @@ This release changes default outputs; regenerate anything calibrated against ear
 
 ### Added
 
+- **Named profiles.** `myogen.load_profile()` loads the `fdi_semg` card: every validated setting of an FDI surface-EMG simulation in one YAML file (`myogen/profiles/fdi_semg.yaml`), each value with its source and, where it matters, a note on what else to check when it is changed. The card builds the library objects (`card.thresholds()`, `card.muscle()`, `card.force_model()`, `card.surface_emg(muscle, electrodes)`, `card.motor_neuron_pool()`) and runs the calibrated descending drive in % MVC (`card.simulate([(duration_s, start_pct_mvc, end_pct_mvc), ...])`). Keywords override card values; an override stays on that card object for every later build, and passing the card value clears it. Overriding a value the drive calibration or the λ fit depends on prints the card's note once per card and field per run (`ProfileNote` warnings), and every override is recorded (`card.record()`, saved on built objects as `profile_record` with the build's seed and passed-through keywords, and on spike and surface-EMG blocks as the `myogen_profile` annotation). Each build draws from the card's seed for that stage (or `seed=`), so it does not depend on what ran before it, and then gives the caller's generator and seed back; this differs from the rest of MyoGen, where library code never seeds and the script owns the seed. `load_profile("fdi_semg.published")` is a variant with each component as its source paper describes it (Fuglevand 1993 force-frequency, Rosenfalck λ = 1, no fibre dispersion). Class defaults are unchanged, and classes built by hand print nothing about cards. The card is validated against frozen limits by `benchmarks/profiles/validate_card.py`. PyYAML is now a declared dependency (`import myogen` loads the cards).
+- **Calibrated descending drive in the library.** `myogen.simulator.neuron.descending_drive_network` (`build_pools`, `run_descending_drive`, `excitation_to_drive_hz`) holds the drive network the FDI calibration fitted; the optimiser imports it.
 - `SurfaceEMG(fiber_diameter__um=57.9, intracellular_conductivity__S_m=1.01)` set the source strength; amplitude scales with the fibre cross-section.
 - `SurfaceEMG(iap_lambda__per_mm=1.0)` sets the length scale of the Rosenfalck IAP, `Vm = 96 (λz)³ e^{−λz} − 90` mV. The default, 1 / mm, is the published value and leaves outputs unchanged. Values below 1 stretch the IAP, lower the MUAP's spectral content and raise surface amplitude. On the FDI defaults, λ = 0.7 matched the median frequency (~175 vs ~165 Hz measured) and adjacent-pin RMS of the Delsys recordings of Botelho et al. (2019) at 25 % MVC; λ = 1 gives ~200 Hz. Botelho et al. used 0.5 in their own FE model for the same reason. Any value other than 1 is a fitted modelling assumption, from one subject. Supported on the default and `execution=` routes; `use_unified=True` refuses λ ≠ 1.
 - **Fibre dispersion in surface EMG.** `SurfaceEMG(fiber_dispersion=SurfaceFiberDispersion(...))` adds, each switched on independently:
@@ -85,6 +87,7 @@ This release changes default outputs; regenerate anything calibrated against ear
 
 ### Changed
 
+- **Basic examples use the `fdi_semg` card.** Examples 01, 03, 04, 05, 06 and 08 build their components with the card, so their outputs differ from earlier releases: 120 FDI units, the calibrated drive in % MVC (03 replaces its Poisson drive loop; 08 uses the spike trains of 03 instead of current injection), the card's muscle and surface-EMG settings.
 - **End-plate spreads in intramuscular EMG.** The two-layer end-plate model now uses the spreads printed in Akhmadeev et al. (2020, eq. 5-6) instead of the 1.5–2× wider values of their reference code: branch-point SD 1.0 + 2.5·s mm (was 1.5 + 4.0·s) and arborisation SD 0.25 + 1.0·s mm (was 0.5 + 1.5·s), with s the share of total threshold below the unit. In the sEMG limb sweep the realised within-unit end-plate SD fell from a median of 1.68 mm to 1.07 mm. iEMG MUAPs and signals change.
 - Citation corrected: the two-layer end-plate and recruitment-threshold model cited as "Konstantin et al. 2020" is Akhmadeev, Yu, Le Carpentier, Aoustin & Farina (2020); Konstantin is the first author's given name. The `konstantin` threshold mode keeps its name.
 - **Default surface fibre diameter (FDI).** `SurfaceEMG(fiber_diameter__um)` 50 → 57.9 µm, the root-mean-square FDI diameter from Polgar et al. (1973): type I 50.0 µm, type II 67.0 µm, 57.4 % type I. The RMS is the right single value because the source scales with d². Default surface amplitudes rise by (57.9/50)² = 1.34; see Breaking changes. `Muscle` fibres are unchanged unless `fiber_model="typed"`.
diff --git a/benchmarks/profiles/README.md b/benchmarks/profiles/README.md
new file mode 100644
index 00000000..ed1439bf
--- /dev/null
+++ b/benchmarks/profiles/README.md
@@ -0,0 +1,36 @@
+# Profile card validation
+
+`validate_card.py` runs a profile card end to end and checks it against limits frozen before the run.
+
+```bash
+python benchmarks/profiles/validate_card.py fdi_semg benchmarks/profiles/results/fdi_semg.json
+python benchmarks/profiles/validate_card.py --mvc-mapping fdi_semg.published benchmarks/profiles/results/fdi_semg.published_mvc.json
+python benchmarks/profiles/validate_card.py fdi_semg.published benchmarks/profiles/results/fdi_semg.published.json
+```
+
+Run at Idle priority with two logical cores left free; `--n-jobs` (default 16) sets the MUAP workers.
+
+## What it checks (`fdi_semg`)
+
+- **Protocol**: `card.simulate` at 25 % MVC with the sEMG benchmark timing (0.5 s rest, 2 s ramp, 4 s hold, 2 s ramp,
+  0.5 s rest), the calibrated drive.
+- **Recording**: `card.muscle()`, `card.surface_emg` with eight Delsys 5-pin sensors (z 3 / 9 mm, rotation 0 / 45°,
+  θ 0 / 180°; 3 × 3 monopolar grids, 2.5 mm pitch, 0.25 mm pins) at 10240 Hz, as arm F of
+  `benchmarks/semg_dispersion`.
+- **Score** (`delsys.py`): corner–centre (cc) and adjacent-corner (adj) channels, 20–2000 Hz, plateau RMS and median
+  frequency, medians over sensors × channels.
+- **Limits** (frozen by Devon 2026-11-03, from arm F re-scored with `delsys.py`): cc and adj RMS within ±10 % of
+  0.196 / 0.233 mV, median MDF within ±5 % of 162.5 / 157.5 Hz. The membrane noise of the drive is not seeded
+  (NEURON `normrand`), so a rerun differs slightly; the limits allow for that.
+
+The published variant is run without limits (its components are the published ones, not the fitted ones); its result
+records what it gives.
+
+## MVC mapping
+
+`--mvc-mapping` runs the card's drive on the MVC ramp of `examples/02_finetune/07_fdi_spikes_force.py` (0.5 s rest,
+excitation 0 → 1 at 0.1 / s, 3 s hold; MVC = mean force over the last 2 s; force smoothed over 0.25 s) with the card's
+force model and reports `drive.mvc_model_units` and `drive.excitation_at_pct_mvc`. The `fdi_semg.published` values
+were measured this way, since its force-frequency curve differs from the calibrated card's.
+
+Results are in `results/`; each card's `validation:` block points to its file.
diff --git a/benchmarks/profiles/delsys.py b/benchmarks/profiles/delsys.py
new file mode 100644
index 00000000..1157804b
--- /dev/null
+++ b/benchmarks/profiles/delsys.py
@@ -0,0 +1,38 @@
+"""Delsys 5-pin scorer: plateau RMS and median frequency of the corner-centre (cc) and adjacent-corner (adj) channels.
+
+The scorer the fdi_semg validation limits were frozen with (arm F of benchmarks/semg_dispersion, 2026-11-03).
+``mono_emg`` is monopolar EMG shaped (time, sensors, 3, 3): each sensor a 3 x 3 grid whose corners and centre are
+the Delsys pins. Each channel is band-passed (4th-order Butterworth, 20-2000 Hz, zero phase), cut to the plateau
+(0.5 + ramp + 0.5 s to 0.5 + ramp + hold - 0.5 s), and scored: RMS, and the median frequency of the Welch spectrum
+(nperseg 4096) over 20-2000 Hz. The result is the median over sensors x channels.
+"""
+from __future__ import annotations
+
+import numpy as np
+from scipy.signal import butter, sosfiltfilt, welch
+
+CORNERS = [(0, 0), (0, 2), (2, 2), (2, 0)]
+PAIRS = {"cc": [(c, (1, 1)) for c in CORNERS], "adj": [(CORNERS[i], CORNERS[(i + 1) % 4]) for i in range(4)]}
+BAND__HZ = (20.0, 2000.0)
+
+
+def score(mono_emg: np.ndarray, fs: float, ramp_s: float, hold_s: float) -> dict[str, tuple[float, float]]:
+    """``{"cc": (median plateau RMS, median MDF Hz), "adj": (...)}``; RMS in the units of ``mono_emg``."""
+    x = np.asarray(mono_emg, dtype=float)
+    if x.ndim != 4 or x.shape[2:] != (3, 3):
+        raise ValueError(f"mono_emg must be shaped (time, sensors, 3, 3), got {x.shape}")
+    sos = butter(4, BAND__HZ, btype="bandpass", fs=fs, output="sos")
+    t = np.arange(x.shape[0]) / fs
+    plateau = (t > 0.5 + ramp_s + 0.5) & (t < 0.5 + ramp_s + hold_s - 0.5)
+    out = {}
+    for kind, pairs in PAIRS.items():
+        rms, mdf = [], []
+        for k in range(x.shape[1]):
+            for p, q in pairs:
+                y = sosfiltfilt(sos, x[:, k, p[0], p[1]] - x[:, k, q[0], q[1]])[plateau]
+                f, P = welch(y, fs=fs, nperseg=4096)
+                band = (f >= BAND__HZ[0]) & (f <= BAND__HZ[1])
+                rms.append(np.sqrt(np.mean(y**2)))
+                mdf.append(f[band][np.searchsorted(np.cumsum(P[band]), P[band].sum() / 2)])
+        out[kind] = (float(np.median(rms)), float(np.median(mdf)))
+    return out
diff --git a/benchmarks/profiles/results/fdi_semg.json b/benchmarks/profiles/results/fdi_semg.json
new file mode 100644
index 00000000..59bd83a1
--- /dev/null
+++ b/benchmarks/profiles/results/fdi_semg.json
@@ -0,0 +1,117 @@
+{
+  "card": {
+    "profile": "fdi_semg",
+    "version": 1,
+    "parent": null,
+    "parent_version": null,
+    "overrides": {}
+  },
+  "date": "2026-11-07",
+  "git_commit": "a6fabd87d741a1bbfc7735d18703c0e7e683e07b",
+  "script": "benchmarks/profiles/validate_card.py",
+  "limits_frozen": "2026-11-03",
+  "protocol": [
+    [
+      0.5,
+      0,
+      0
+    ],
+    [
+      2.0,
+      0,
+      25.0
+    ],
+    [
+      4.0,
+      25.0,
+      25.0
+    ],
+    [
+      2.0,
+      25.0,
+      0
+    ],
+    [
+      0.5,
+      0,
+      0
+    ]
+  ],
+  "fs_Hz": 10240.0,
+  "sensors": [
+    [
+      3.0,
+      0.0,
+      0.0
+    ],
+    [
+      3.0,
+      0.0,
+      180.0
+    ],
+    [
+      3.0,
+      45.0,
+      0.0
+    ],
+    [
+      3.0,
+      45.0,
+      180.0
+    ],
+    [
+      9.0,
+      0.0,
+      0.0
+    ],
+    [
+      9.0,
+      0.0,
+      180.0
+    ],
+    [
+      9.0,
+      45.0,
+      0.0
+    ],
+    [
+      9.0,
+      45.0,
+      180.0
+    ]
+  ],
+  "n_units_simulated": 104,
+  "n_active": 104,
+  "hold_force_pct_mvc": 25.621956085825367,
+  "scores": {
+    "cc": {
+      "rms_mV": 0.20043793229022705,
+      "mdf_Hz": 161.25,
+      "rms_limits_mV": [
+        0.1764,
+        0.2156
+      ],
+      "mdf_limits_Hz": [
+        154.375,
+        170.625
+      ],
+      "passed": true
+    },
+    "adj": {
+      "rms_mV": 0.23404846872880228,
+      "mdf_Hz": 156.25,
+      "rms_limits_mV": [
+        0.20970000000000003,
+        0.25630000000000003
+      ],
+      "mdf_limits_Hz": [
+        149.625,
+        165.375
+      ],
+      "passed": true
+    }
+  },
+  "checked": true,
+  "passed": true,
+  "runtime_s": 1355.3
+}
diff --git a/benchmarks/profiles/results/fdi_semg.published.json b/benchmarks/profiles/results/fdi_semg.published.json
new file mode 100644
index 00000000..e04c74f5
--- /dev/null
+++ b/benchmarks/profiles/results/fdi_semg.published.json
@@ -0,0 +1,99 @@
+{
+  "card": {
+    "profile": "fdi_semg.published",
+    "version": 1,
+    "parent": "fdi_semg",
+    "parent_version": 1,
+    "overrides": {}
+  },
+  "date": "2026-11-07",
+  "git_commit": "a6fabd87d741a1bbfc7735d18703c0e7e683e07b",
+  "script": "benchmarks/profiles/validate_card.py",
+  "limits_frozen": "2026-11-03",
+  "protocol": [
+    [
+      0.5,
+      0,
+      0
+    ],
+    [
+      2.0,
+      0,
+      25.0
+    ],
+    [
+      4.0,
+      25.0,
+      25.0
+    ],
+    [
+      2.0,
+      25.0,
+      0
+    ],
+    [
+      0.5,
+      0,
+      0
+    ]
+  ],
+  "fs_Hz": 10240.0,
+  "sensors": [
+    [
+      3.0,
+      0.0,
+      0.0
+    ],
+    [
+      3.0,
+      0.0,
+      180.0
+    ],
+    [
+      3.0,
+      45.0,
+      0.0
+    ],
+    [
+      3.0,
+      45.0,
+      180.0
+    ],
+    [
+      9.0,
+      0.0,
+      0.0
+    ],
+    [
+      9.0,
+      0.0,
+      180.0
+    ],
+    [
+      9.0,
+      45.0,
+      0.0
+    ],
+    [
+      9.0,
+      45.0,
+      180.0
+    ]
+  ],
+  "n_units_simulated": 103,
+  "n_active": 103,
+  "hold_force_pct_mvc": 25.490816185125954,
+  "scores": {
+    "cc": {
+      "rms_mV": 0.2470491006474162,
+      "mdf_Hz": 227.5
+    },
+    "adj": {
+      "rms_mV": 0.30577169099806145,
+      "mdf_Hz": 220.0
+    }
+  },
+  "checked": true,
+  "passed": null,
+  "runtime_s": 1275.5
+}
diff --git a/benchmarks/profiles/results/fdi_semg.published_mvc.json b/benchmarks/profiles/results/fdi_semg.published_mvc.json
new file mode 100644
index 00000000..8a05bd8c
--- /dev/null
+++ b/benchmarks/profiles/results/fdi_semg.published_mvc.json
@@ -0,0 +1,20 @@
+{
+  "card": {
+    "profile": "fdi_semg.published",
+    "version": 1,
+    "parent": "fdi_semg",
+    "overrides": {}
+  },
+  "date": "2026-11-04",
+  "git_commit": "2d4911c15b7fa6700e6af77accad072f104fc115",
+  "script": "benchmarks/profiles/validate_card.py --mvc-mapping",
+  "mvc_model_units": 13470.808768515415,
+  "excitation_at_pct_mvc": {
+    "5": 0.05927734375,
+    "10": 0.09799804687500001,
+    "20": 0.153955078125,
+    "25": 0.184521484375,
+    "30": 0.212890625,
+    "50": 0.31396484375
+  }
+}
diff --git a/benchmarks/profiles/validate_card.py b/benchmarks/profiles/validate_card.py
new file mode 100644
index 00000000..44f93b83
--- /dev/null
+++ b/benchmarks/profiles/validate_card.py
@@ -0,0 +1,157 @@
+"""Validate a profile card against its frozen limits, or measure a card's MVC mapping.
+
+    python validate_card.py fdi_semg results/fdi_semg.json [--n-jobs 16]
+    python validate_card.py --mvc-mapping fdi_semg.published results/fdi_semg.published_mvc.json
+
+Validation runs the card end to end, as a user would: ``card.simulate`` (the calibrated drive, a 25 % MVC trapezoid
+with the sEMG benchmark timing: 0.5 s rest, 2 s ramp, 4 s hold, 2 s ramp, 0.5 s rest), ``card.muscle()``,
+``card.surface_emg`` with eight Delsys 5-pin sensors (z 3 / 9 mm, rotation 0 / 45 deg, theta 0 / 180 deg; 3 x 3
+monopolar grids, 2.5 mm pitch, 0.25 mm pins, as Pereira Botelho et al. 2019 and arm F of
+benchmarks/semg_dispersion) at 10240 Hz, then scores cc / adj plateau RMS and median frequency (delsys.py). The
+result JSON holds the values, the limits, pass / fail, the card record and the git commit.
+
+10240 Hz is an analysis setting of this comparison, not a recording rate to copy: it lets the scorer use the
+20-2000 Hz band of the Delsys comparator, on which the limits were frozen. Generated sEMG is normally sampled at
+2048 Hz (the ``SurfaceEMG`` default; iEMG at 10240 Hz).
+
+The MVC mapping runs the card's drive on the MVC ramp of examples/02_finetune/07_fdi_spikes_force.py (0.5 s rest,
+excitation 0 -> 1 at 0.1 / s, 3 s hold; MVC = mean force over the last 2 s; force smoothed over 0.25 s) with the
+card's force model, and reports the MVC and the excitation at which the ramp force first reaches 5, 10, 20, 25, 30
+and 50 % MVC: the values of ``drive.mvc_model_units`` and ``drive.excitation_at_pct_mvc``.
+"""
+from __future__ import annotations
+
+import argparse
+import datetime
+import json
+import subprocess
+import sys
+import time
+from pathlib import Path
+
+import numpy as np
+import quantities as pq
+
+HERE = Path(__file__).resolve().parent
+sys.path.insert(0, str(HERE))
+import delsys  # noqa: E402
+
+from myogen import load_profile, set_random_seed  # noqa: E402
+from myogen.simulator import SurfaceElectrodeArray  # noqa: E402
+from myogen.utils.neo import signal_to_grid  # noqa: E402
+from myogen.utils.nmodl import load_nmodl_mechanisms  # noqa: E402
+
+# Frozen by Devon 2026-11-03, before this script existed: arm F of benchmarks/semg_dispersion re-scored with
+# delsys.py; cc / adj plateau RMS within +-10 %, median MDF within +-5 %.
+LIMITS_FROZEN = "2026-11-03"
+REFERENCE = {"cc": (0.196, 162.5), "adj": (0.233, 157.5)}
+RMS_TOL, MDF_TOL = 0.10, 0.05
+LIMITED_CARDS = {"fdi_semg"}  # other cards (e.g. the published variant) are run and recorded, not judged
+
+FS__HZ = 10240.0  # for the comparator's 20-2000 Hz band only; generated sEMG is normally 2048 Hz (see above)
+RAMP_S, HOLD_S, LEVEL_PCT_MVC = 2.0, 4.0, 25.0
+SENSORS = [(z, rot, th) for z in (3.0, 9.0) for rot in (0.0, 45.0) for th in (0.0, 180.0)]
+
+MVC_REST_S, MVC_RAMP_PER_S, MVC_HOLD_S, MVC_LAST_S = 0.5, 0.1, 3.0, 2.0
+FORCE_SMOOTH_S = 0.25
+MVC_LEVELS = (5, 10, 20, 25, 30, 50)
+
+
+def git_commit() -> str:
+    try:
+        return subprocess.run(["git", "-C", str(HERE), "rev-parse", "HEAD"], capture_output=True, text=True,
+                              check=True).stdout.strip()
+    except (OSError, subprocess.CalledProcessError):
+        return "unknown"
+
+
+def validate(name: str, n_jobs: int) -> dict:
+    card = load_profile(name)
+    protocol = [(0.5, 0, 0), (RAMP_S, 0, LEVEL_PCT_MVC), (HOLD_S, LEVEL_PCT_MVC, LEVEL_PCT_MVC),
+                (RAMP_S, LEVEL_PCT_MVC, 0), (0.5, 0, 0)]
+    t0 = time.time()
+    run = card.simulate(protocol)
+    trains = run.spikes.segments[0].spiketrains
+    muscle = card.muscle(n_jobs=n_jobs)
+    counts = np.asarray(muscle.resulting_number_of_innervated_fibers)
+    mus = [m for m, st in enumerate(trains) if len(st) > 0 and counts[m] > 0]
+    bend = (card.value("muscle.radius__mm") + card.value("muscle.fat_thickness__mm")
+            + card.value("muscle.skin_thickness__mm")) * pq.mm
+    arrays = [SurfaceElectrodeArray(3, 3, 2.5 * pq.mm, 0.25 * pq.mm, center_point__mm_deg=(z * pq.mm, th * pq.deg),
+                                    bending_radius__mm=bend, rotation_angle__deg=rot * pq.deg,
+                                    differentiation_mode="monopolar") for z, rot, th in SENSORS]
+    semg = card.surface_emg(muscle, arrays, sampling_frequency__Hz=FS__HZ * pq.Hz,
+                            internal_sampling_frequency__Hz=FS__HZ * pq.Hz, MUs_to_simulate=mus)  # as arm F
+    semg.simulate_muaps(n_jobs=n_jobs, verbose=False, use_gpu=False)
+    emg = semg.simulate_surface_emg(run.spikes, verbose=False)
+    mono = np.stack([np.asarray(signal_to_grid(g.segments[0].analogsignals[0].rescale(pq.mV)), dtype=float)
+                     for g in emg.groups], axis=1)
+    scores = delsys.score(mono, FS__HZ, RAMP_S, HOLD_S)
+    hold = (np.arange(run.force_pct_mvc.size) / run.fs_force__Hz > 0.5 + RAMP_S + 0.5) & \
+           (np.arange(run.force_pct_mvc.size) / run.fs_force__Hz < 0.5 + RAMP_S + HOLD_S - 0.5)
+    checks = {}
+    for kind, (rms, mdf) in scores.items():
+        ref_rms, ref_mdf = REFERENCE[kind]
+        checks[kind] = {"rms_mV": rms, "mdf_Hz": mdf}
+        if name in LIMITED_CARDS:
+            checks[kind]["rms_limits_mV"] = [ref_rms * (1 - RMS_TOL), ref_rms * (1 + RMS_TOL)]
+            checks[kind]["mdf_limits_Hz"] = [ref_mdf * (1 - MDF_TOL), ref_mdf * (1 + MDF_TOL)]
+            checks[kind]["passed"] = bool(
+                checks[kind]["rms_limits_mV"][0] <= rms <= checks[kind]["rms_limits_mV"][1]
+                and checks[kind]["mdf_limits_Hz"][0] <= mdf <= checks[kind]["mdf_limits_Hz"][1])
+    return {"card": card.record(), "date": datetime.date.today().isoformat(), "git_commit": git_commit(),
+            "script": "benchmarks/profiles/validate_card.py", "limits_frozen": LIMITS_FROZEN,
+            "protocol": protocol, "fs_Hz": FS__HZ, "sensors": SENSORS, "n_units_simulated": len(mus),
+            "n_active": int(sum(len(st) > 0 for st in trains)),
+            "hold_force_pct_mvc": float(np.mean(run.force_pct_mvc[hold])),
+            "scores": checks, "checked": True,
+            "passed": all(c["passed"] for c in checks.values()) if name in LIMITED_CARDS else None,
+            "runtime_s": round(time.time() - t0, 1)}
+
+
+def mvc_mapping(name: str) -> dict:
+    from myogen.simulator.neuron.descending_drive_network import excitation_to_drive_hz, run_descending_drive
+
+    card = load_profile(name)
+    thr = card.thresholds()
+    params = card.dd_params()
+    dt = float(card.value("drive.timestep_ms"))
+    segments = [(MVC_REST_S, 0.0, 0.0), (1.0 / MVC_RAMP_PER_S, 0.0, 1.0), (MVC_HOLD_S, 1.0, 1.0)]
+    e, hz = excitation_to_drive_hz(segments, float(card.value("drive.threshold_hz")), params["dd_drive__Hz"],
+                                   timestep_ms=dt, noise_seed=int(card.seeds["drive"]))
+    set_random_seed(int(card.seeds["drive"]))
+    block = run_descending_drive(thr, hz, params, len(e) * dt, timestep_ms=dt,
+                                 config_file=card.value("motor_neurons.config_file"))
+    fm = card.force_model()
+    fs = float(card.value("force.recording_frequency__Hz"))
+    force = np.asarray(fm.generate_force(block, verbose=False).magnitude, dtype=float)[:, 0]
+    t = np.arange(force.size) / fs
+    mvc = float(np.mean(force[t >= t[-1] - MVC_LAST_S]))
+    w = max(int(FORCE_SMOOTH_S * fs), 1)
+    smooth = np.convolve(force, np.ones(w) / w, mode="same")
+    e_t = np.interp(t, np.arange(e.size) * dt / 1000.0, e)
+    ramp = (t > MVC_REST_S) & (t < MVC_REST_S + 1.0 / MVC_RAMP_PER_S)
+    curve = {p: float(e_t[ramp][np.argmax(smooth[ramp] >= p / 100 * mvc)]) for p in MVC_LEVELS}
+    return {"card": card.record(), "date": datetime.date.today().isoformat(), "git_commit": git_commit(),
+            "script": "benchmarks/profiles/validate_card.py --mvc-mapping", "mvc_model_units": mvc,
+            "excitation_at_pct_mvc": curve}
+
+
+def main() -> None:
+    ap = argparse.ArgumentParser(description=__doc__.split("\n\n")[0])
+    ap.add_argument("card")
+    ap.add_argument("out", type=Path)
+    ap.add_argument("--mvc-mapping", action="store_true", help="measure the MVC mapping instead of validating")
+    ap.add_argument("--n-jobs", type=int, default=16)
+    args = ap.parse_args()
+    load_nmodl_mechanisms()
+    result = mvc_mapping(args.card) if args.mvc_mapping else validate(args.card, args.n_jobs)
+    args.out.parent.mkdir(parents=True, exist_ok=True)
+    args.out.write_text(json.dumps(result, indent=2) + "\n", encoding="utf-8")
+    print(json.dumps({k: v for k, v in result.items() if k in ("scores", "passed", "mvc_model_units",
+                                                               "excitation_at_pct_mvc", "hold_force_pct_mvc")},
+                     indent=2))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/docs/literature.md b/docs/literature.md
index f213b802..aaf98ffe 100644
--- a/docs/literature.md
+++ b/docs/literature.md
@@ -43,20 +43,20 @@ To add a source: add it to the table in `scripts/build_literature.py`, link it f
 
 | Source | Status | Use | What it supports | Used in |
 |---|---|---|---|---|
-| <a id="farina2004"></a>**farina2004**: Farina D, Mesin L, Martina S, Merletti R (2004). A surface EMG generation model with multilayer cylindrical description of the volume conductor. IEEE Trans Biomed Eng 51(3):415–426. [doi:10.1109/TBME.2003.820998](https://doi.org/10.1109/TBME.2003.820998) | FT | ADOPTED | The sEMG volume conductor; reference for the notebook fixes (bipolar B(k_z), MUAP time axis, IAP z axis, fibre semi-lengths, amplitude). Layer σ (skin 1, fat 0.05, muscle 0.5/0.1 S/m) given without source. | [SurfaceEMG](../myogen/simulator/core/emg/surface/surface_emg.py), [simulate_fiber](../myogen/simulator/core/emg/surface/simulate_fiber.py), [Muscle](../myogen/simulator/core/muscle/muscle.py), [test_surface_amplitude](../tests/test_surface_amplitude.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [electrodes.py:189](../myogen/simulator/core/emg/electrodes.py#L189), [electrodes.py:323](../myogen/simulator/core/emg/electrodes.py#L323), [electrodes.py:328](../myogen/simulator/core/emg/electrodes.py#L328), [electrodes.py:335](../myogen/simulator/core/emg/electrodes.py#L335), [electrodes.py:340](../myogen/simulator/core/emg/electrodes.py#L340), [fiber_simulation.py:11](../myogen/simulator/core/emg/fiber_simulation.py#L11), [fiber_simulation.py:115](../myogen/simulator/core/emg/fiber_simulation.py#L115), [fiber_simulation.py:418](../myogen/simulator/core/emg/fiber_simulation.py#L418), [bioelectric.py:7](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L7), [simulate_fiber.py:107](../myogen/simulator/core/emg/surface/simulate_fiber.py#L107), [simulate_fiber.py:184](../myogen/simulator/core/emg/surface/simulate_fiber.py#L184), [simulate_fiber.py:452](../myogen/simulator/core/emg/surface/simulate_fiber.py#L452), [surface_emg.py:53](../myogen/simulator/core/emg/surface/surface_emg.py#L53), [surface_emg.py:148](../myogen/simulator/core/emg/surface/surface_emg.py#L148), [volume_conductor.py:1](../myogen/simulator/core/emg/surface/volume_conductor.py#L1), [volume_conductor.py:97](../myogen/simulator/core/emg/surface/volume_conductor.py#L97), [volume_conductor.py:140](../myogen/simulator/core/emg/surface/volume_conductor.py#L140), [muscle.py:145](../myogen/simulator/core/muscle/muscle.py#L145), [muscle.py:243](../myogen/simulator/core/muscle/muscle.py#L243), [muscle.py:295](../myogen/simulator/core/muscle/muscle.py#L295) |
-| <a id="roeleveld1997"></a>**roeleveld1997**: Roeleveld K, Blok JH, Stegeman DF, van Oosterom A (1997). Volume conduction models for surface EMG; confrontation with measurements. J Electromyogr Kinesiol 7(4):221–232. [doi:10.1016/S1050-6411(97)00009-6](https://doi.org/10.1016/S1050-6411(97)00009-6) | FT | ADOPTED | Provenance of skin σ = 1 S/m: a fitted effective value reproducing biceps MUAP fall-off (52 MUs), not a measurement. Replicated in the Roeleveld bench. | [Muscle](../myogen/simulator/core/muscle/muscle.py), [Roeleveld bench](../benchmarks/semg_roeleveld1997/README.md), [test_surface_amplitude](../tests/test_surface_amplitude.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [muscle.py:244](../myogen/simulator/core/muscle/muscle.py#L244), [muscle.py:246](../myogen/simulator/core/muscle/muscle.py#L246), [muscle.py:302](../myogen/simulator/core/muscle/muscle.py#L302) |
+| <a id="farina2004"></a>**farina2004**: Farina D, Mesin L, Martina S, Merletti R (2004). A surface EMG generation model with multilayer cylindrical description of the volume conductor. IEEE Trans Biomed Eng 51(3):415–426. [doi:10.1109/TBME.2003.820998](https://doi.org/10.1109/TBME.2003.820998) | FT | ADOPTED | The sEMG volume conductor; reference for the notebook fixes (bipolar B(k_z), MUAP time axis, IAP z axis, fibre semi-lengths, amplitude). Layer σ (skin 1, fat 0.05, muscle 0.5/0.1 S/m) given without source. | [SurfaceEMG](../myogen/simulator/core/emg/surface/surface_emg.py), [simulate_fiber](../myogen/simulator/core/emg/surface/simulate_fiber.py), [Muscle](../myogen/simulator/core/muscle/muscle.py), [test_surface_amplitude](../tests/test_surface_amplitude.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [fdi_semg.published.yaml:23](../myogen/profiles/fdi_semg.published.yaml#L23), [electrodes.py:189](../myogen/simulator/core/emg/electrodes.py#L189), [electrodes.py:323](../myogen/simulator/core/emg/electrodes.py#L323), [electrodes.py:328](../myogen/simulator/core/emg/electrodes.py#L328), [electrodes.py:335](../myogen/simulator/core/emg/electrodes.py#L335), [electrodes.py:340](../myogen/simulator/core/emg/electrodes.py#L340), [fiber_simulation.py:11](../myogen/simulator/core/emg/fiber_simulation.py#L11), [fiber_simulation.py:115](../myogen/simulator/core/emg/fiber_simulation.py#L115), [fiber_simulation.py:418](../myogen/simulator/core/emg/fiber_simulation.py#L418), [bioelectric.py:7](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L7), [simulate_fiber.py:107](../myogen/simulator/core/emg/surface/simulate_fiber.py#L107), [simulate_fiber.py:184](../myogen/simulator/core/emg/surface/simulate_fiber.py#L184), [simulate_fiber.py:452](../myogen/simulator/core/emg/surface/simulate_fiber.py#L452), [surface_emg.py:53](../myogen/simulator/core/emg/surface/surface_emg.py#L53), [surface_emg.py:148](../myogen/simulator/core/emg/surface/surface_emg.py#L148), [volume_conductor.py:1](../myogen/simulator/core/emg/surface/volume_conductor.py#L1), [volume_conductor.py:97](../myogen/simulator/core/emg/surface/volume_conductor.py#L97), [volume_conductor.py:140](../myogen/simulator/core/emg/surface/volume_conductor.py#L140), [muscle.py:145](../myogen/simulator/core/muscle/muscle.py#L145), [muscle.py:243](../myogen/simulator/core/muscle/muscle.py#L243), [muscle.py:295](../myogen/simulator/core/muscle/muscle.py#L295) |
+| <a id="roeleveld1997"></a>**roeleveld1997**: Roeleveld K, Blok JH, Stegeman DF, van Oosterom A (1997). Volume conduction models for surface EMG; confrontation with measurements. J Electromyogr Kinesiol 7(4):221–232. [doi:10.1016/S1050-6411(97)00009-6](https://doi.org/10.1016/S1050-6411(97)00009-6) | FT | ADOPTED | Provenance of skin σ = 1 S/m: a fitted effective value reproducing biceps MUAP fall-off (52 MUs), not a measurement. Replicated in the Roeleveld bench. | [Muscle](../myogen/simulator/core/muscle/muscle.py), [Roeleveld bench](../benchmarks/semg_roeleveld1997/README.md), [test_surface_amplitude](../tests/test_surface_amplitude.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [fdi_semg.yaml:55](../myogen/profiles/fdi_semg.yaml#L55), [muscle.py:244](../myogen/simulator/core/muscle/muscle.py#L244), [muscle.py:246](../myogen/simulator/core/muscle/muscle.py#L246), [muscle.py:302](../myogen/simulator/core/muscle/muscle.py#L302) |
 | <a id="roeleveld1997b"></a>**roeleveld1997b**: Roeleveld K et al. (1997), companion JEK paper on MUAP spread / exponents (details unverified). | LINK | CONTEXT | Source of the depth/lateral fall-off regression used in the bench; not read. | [Roeleveld bench](../benchmarks/semg_roeleveld1997/README.md) |
 | <a id="blok2002"></a>**blok2002**: Blok JH, Stegeman DF, van Oosterom A (2002). Three-layer volume conductor model and software package for applications in surface electromyography. Ann Biomed Eng 30:566–577. [doi:10.1114/1.1475345](https://doi.org/10.1114/1.1475345) | FT | CONTEXT | Skin 1.0 S/m 'taken from Roeleveld'; calls it 'somewhat ambiguous'. σi 0.5 S/m, d 50 µm (source-strength comparison). | [Muscle](../myogen/simulator/core/muscle/muscle.py)<br>Cited at: [muscle.py:244](../myogen/simulator/core/muscle/muscle.py#L244), [muscle.py:247](../myogen/simulator/core/muscle/muscle.py#L247), [muscle.py:303](../myogen/simulator/core/muscle/muscle.py#L303) |
 | <a id="lowery2002"></a>**lowery2002**: Lowery MM, Stoykov NS, Taflove A, Kuiken TA (2002). A multiple-layer finite-element model of the surface EMG signal. IEEE Trans Biomed Eng 49(5):446–454. [doi:10.1109/10.995683](https://doi.org/10.1109/10.995683) | FT | CONTEXT | Counterpoint: measured skin 4.55e-4 S/m (Gabriel 1996); Roeleveld's skin 'considerably higher than measured'. σi 1.01 S/m (Table I). | [Muscle](../myogen/simulator/core/muscle/muscle.py), [SurfaceEMG](../myogen/simulator/core/emg/surface/surface_emg.py)<br>Cited at: [muscle.py:249](../myogen/simulator/core/muscle/muscle.py#L249), [muscle.py:304](../myogen/simulator/core/muscle/muscle.py#L304) |
 | <a id="farinarainoldi1999"></a>**farinarainoldi1999**: Farina D, Rainoldi A (1999). Compensation of the effect of sub-cutaneous tissue layers on surface EMG: a simulation study. Med Eng Phys 21:487–496. [doi:10.1016/S1350-4533(99)00075-2](https://doi.org/10.1016/S1350-4533(99)00075-2) | FT | CONTEXT | Skin/fat ratio 20 (7 also tested), fat/muscle 0.5, cited to Roeleveld (EEG abstract), Geddes & Baker, Gootzen. 20 × 0.05 = 1 S/m. | — |
 | <a id="farinamerletti2001"></a>**farinamerletti2001**: Farina D, Merletti R (2001). A novel approach for precise simulation of the EMG signal detected by surface electrodes. IEEE Trans Biomed Eng 48(6):637–646. [doi:10.1109/10.923782](https://doi.org/10.1109/10.923782) · [PMID 11396594](https://pubmed.ncbi.nlm.nih.gov/11396594/) | FT | CONTEXT | Same conductivity ratios; IAP 96 mV/mm³ with λ = 1 / mm. | [simulate_fiber](../myogen/simulator/core/emg/surface/simulate_fiber.py)<br>Cited at: [bioelectric.py:12](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L12), [bioelectric.py:232](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L232), [simulate_fiber.py:463](../myogen/simulator/core/emg/surface/simulate_fiber.py#L463) |
 | <a id="merlettifarina2016"></a>**merlettifarina2016**: Merletti R, Farina D (eds) (2016). Surface Electromyography: Physiology, Engineering and Applications. Wiley-IEEE Press (ch. 8, Lowery). | FT | CONTEXT | §8.6.3 gives muscle 0.093–0.8 and fat 0.02–0.1 S/m with sources, but no skin value. | — |
-| <a id="botelho2019"></a>**botelho2019**: Pereira Botelho D, Curran K, Lowery MM (2019). Anatomically accurate model of EMG during index finger flexion and abduction derived from diffusion tensor imaging. PLoS Comput Biol 15(8):e1007267. [doi:10.1371/journal.pcbi.1007267](https://doi.org/10.1371/journal.pcbi.1007267) · [PMC6738720](https://pmc.ncbi.nlm.nih.gov/articles/PMC6738720/) | FT + FIG | CONTEXT | Main FDI comparator, measured: RMS 0.21 mV at 25 % MVC (IQR 0.19–0.26), MDF ~165 Hz (IQR 152–185), one subject, Delsys 5-pin. Existing muscle/fat σ defaults. Rejected as defaults: dry skin 4.88e-4 S/m; Rosenfalck α 0.5 (precedent for the optional λ). | [Muscle](../myogen/simulator/core/muscle/muscle.py), [SurfaceEMG](../myogen/simulator/core/emg/surface/surface_emg.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [surface_emg.py:108](../myogen/simulator/core/emg/surface/surface_emg.py#L108), [surface_emg.py:112](../myogen/simulator/core/emg/surface/surface_emg.py#L112), [surface_emg.py:118](../myogen/simulator/core/emg/surface/surface_emg.py#L118), [surface_emg.py:335](../myogen/simulator/core/emg/surface/surface_emg.py#L335), [muscle.py:160](../myogen/simulator/core/muscle/muscle.py#L160), [muscle.py:237](../myogen/simulator/core/muscle/muscle.py#L237), [muscle.py:239](../myogen/simulator/core/muscle/muscle.py#L239), [muscle.py:241](../myogen/simulator/core/muscle/muscle.py#L241), [muscle.py:250](../myogen/simulator/core/muscle/muscle.py#L250), [muscle.py:301](../myogen/simulator/core/muscle/muscle.py#L301) |
+| <a id="botelho2019"></a>**botelho2019**: Pereira Botelho D, Curran K, Lowery MM (2019). Anatomically accurate model of EMG during index finger flexion and abduction derived from diffusion tensor imaging. PLoS Comput Biol 15(8):e1007267. [doi:10.1371/journal.pcbi.1007267](https://doi.org/10.1371/journal.pcbi.1007267) · [PMC6738720](https://pmc.ncbi.nlm.nih.gov/articles/PMC6738720/) | FT + FIG | CONTEXT | Main FDI comparator, measured: RMS 0.21 mV at 25 % MVC (IQR 0.19–0.26), MDF ~165 Hz (IQR 152–185), one subject, Delsys 5-pin. Existing muscle/fat σ defaults. Rejected as defaults: dry skin 4.88e-4 S/m; Rosenfalck α 0.5 (precedent for the optional λ). | [Muscle](../myogen/simulator/core/muscle/muscle.py), [SurfaceEMG](../myogen/simulator/core/emg/surface/surface_emg.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [fdi_semg.yaml:61](../myogen/profiles/fdi_semg.yaml#L61), [surface_emg.py:108](../myogen/simulator/core/emg/surface/surface_emg.py#L108), [surface_emg.py:112](../myogen/simulator/core/emg/surface/surface_emg.py#L112), [surface_emg.py:118](../myogen/simulator/core/emg/surface/surface_emg.py#L118), [surface_emg.py:343](../myogen/simulator/core/emg/surface/surface_emg.py#L343), [muscle.py:160](../myogen/simulator/core/muscle/muscle.py#L160), [muscle.py:237](../myogen/simulator/core/muscle/muscle.py#L237), [muscle.py:239](../myogen/simulator/core/muscle/muscle.py#L239), [muscle.py:241](../myogen/simulator/core/muscle/muscle.py#L241), [muscle.py:250](../myogen/simulator/core/muscle/muscle.py#L250), [muscle.py:301](../myogen/simulator/core/muscle/muscle.py#L301) |
 | <a id="gabriel1996"></a>**gabriel1996**: Gabriel S, Lau RW, Gabriel C (1996). The dielectric properties of biological tissues (series). Phys Med Biol 41. | SEC | CONTEXT | Dry skin ~2–5e-4 S/m; source of Lowery's and Botelho's skin values. | [Muscle](../myogen/simulator/core/muscle/muscle.py)<br>Cited at: [muscle.py:249](../myogen/simulator/core/muscle/muscle.py#L249) |
 | <a id="itis"></a>**itis**: Baumgartner C, Hasgall PA, Di Gennaro F, Neufeld E, Lloyd B, Gosselin MC, Payne D, Klingenböck A, Kuster N (2025). IT'IS Database for thermal and electromagnetic parameters of biological tissues, Version 5.0 (low-frequency conductivity table). [doi:10.13099/VIP21000-05-0](https://doi.org/10.13099/VIP21000-05-0) | FT | CONTEXT | Low-frequency skin 0.148 ± 0.042 S/m (7 studies, 0.057–0.25; the entry pools dry skin and epidermis/dermis, and the site advises caution). The 0.15 S/m of the calibration argument for 1 S/m; 0.15 S/m in the cylinder gives about 2.5–4× too much amplitude. | [Muscle](../myogen/simulator/core/muscle/muscle.py) |
 | <a id="dideriksen2019"></a>**dideriksen2019**: Dideriksen JL, Farina D (2019). PLoS Comput Biol (amplitude cancellation and neural drive; title unverified). | FT | CONTEXT | Uses skin 2 mm / 1 S/m, attributed to 'a previously proposed model'. | — |
 | <a id="hodgkinnakajima1972"></a>**hodgkinnakajima1972**: Hodgkin AL, Nakajima S (1972). The effect of diameter on the electrical constants of frog skeletal muscle fibres. J Physiol 221:105. [doi:10.1113/jphysiol.1972.sp009742](https://doi.org/10.1113/jphysiol.1972.sp009742) | ABS | CONTEXT | Gi 0.59 S/m at 20 °C (≈1.0 S/m at 37 °C, Q10 1.37, our arithmetic): supports σi = 1.01 S/m in the line-source strength. | [SurfaceEMG](../myogen/simulator/core/emg/surface/surface_emg.py) |
-| <a id="rosenfalck1969"></a>**rosenfalck1969**: Rosenfalck P (1969). Intra- and extracellular potential fields of active nerve and muscle fibres. Acta Physiol Scand Suppl 321. [PMID 5383732](https://pubmed.ncbi.nlm.nih.gov/5383732/) | SEC | ADOPTED | IAP Vm = 96 (λz)³ e^(−λz) − 90 mV, λ = 1 / mm (default; Devon 2026-11-01). λ exposed as SurfaceEMG(iap_lambda__per_mm). FDI fit λ 0.7 documented, not adopted. | [simulate_fiber](../myogen/simulator/core/emg/surface/simulate_fiber.py), [SurfaceEMG](../myogen/simulator/core/emg/surface/surface_emg.py), [test_surface_amplitude](../tests/test_surface_amplitude.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [fiber_simulation.py:5](../myogen/simulator/core/emg/fiber_simulation.py#L5), [fiber_simulation.py:10](../myogen/simulator/core/emg/fiber_simulation.py#L10), [fiber_simulation.py:40](../myogen/simulator/core/emg/fiber_simulation.py#L40), [fiber_simulation.py:287](../myogen/simulator/core/emg/fiber_simulation.py#L287), [fiber_simulation.py:330](../myogen/simulator/core/emg/fiber_simulation.py#L330), [fiber_simulation.py:389](../myogen/simulator/core/emg/fiber_simulation.py#L389), [fiber_simulation.py:553](../myogen/simulator/core/emg/fiber_simulation.py#L553), [fiber_simulation.py:568](../myogen/simulator/core/emg/fiber_simulation.py#L568), [bioelectric.py:8](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L8), [bioelectric.py:13](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L13), [bioelectric.py:25](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L25), [bioelectric.py:48](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L48), [bioelectric.py:76](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L76), [simulate_fiber.py:381](../myogen/simulator/core/emg/surface/simulate_fiber.py#L381), [simulate_fiber.py:383](../myogen/simulator/core/emg/surface/simulate_fiber.py#L383), [simulate_fiber.py:461](../myogen/simulator/core/emg/surface/simulate_fiber.py#L461), [simulate_fiber.py:478](../myogen/simulator/core/emg/surface/simulate_fiber.py#L478), [simulate_fiber.py:783](../myogen/simulator/core/emg/surface/simulate_fiber.py#L783), [simulate_fiber.py:785](../myogen/simulator/core/emg/surface/simulate_fiber.py#L785), [surface_emg.py:104](../myogen/simulator/core/emg/surface/surface_emg.py#L104), [surface_emg.py:106](../myogen/simulator/core/emg/surface/surface_emg.py#L106), [surface_emg.py:571](../myogen/simulator/core/emg/surface/surface_emg.py#L571), [test_iemg_sfap_units.py:4](../tests/test_iemg_sfap_units.py#L4), [test_iemg_sfap_units.py:5](../tests/test_iemg_sfap_units.py#L5) |
+| <a id="rosenfalck1969"></a>**rosenfalck1969**: Rosenfalck P (1969). Intra- and extracellular potential fields of active nerve and muscle fibres. Acta Physiol Scand Suppl 321. [PMID 5383732](https://pubmed.ncbi.nlm.nih.gov/5383732/) | SEC | ADOPTED | IAP Vm = 96 (λz)³ e^(−λz) − 90 mV, λ = 1 / mm (default; Devon 2026-11-01). λ exposed as SurfaceEMG(iap_lambda__per_mm). FDI fit λ 0.7 documented, not adopted. | [simulate_fiber](../myogen/simulator/core/emg/surface/simulate_fiber.py), [SurfaceEMG](../myogen/simulator/core/emg/surface/surface_emg.py), [test_surface_amplitude](../tests/test_surface_amplitude.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [fdi_semg.published.yaml:24](../myogen/profiles/fdi_semg.published.yaml#L24), [fiber_simulation.py:5](../myogen/simulator/core/emg/fiber_simulation.py#L5), [fiber_simulation.py:10](../myogen/simulator/core/emg/fiber_simulation.py#L10), [fiber_simulation.py:40](../myogen/simulator/core/emg/fiber_simulation.py#L40), [fiber_simulation.py:287](../myogen/simulator/core/emg/fiber_simulation.py#L287), [fiber_simulation.py:330](../myogen/simulator/core/emg/fiber_simulation.py#L330), [fiber_simulation.py:389](../myogen/simulator/core/emg/fiber_simulation.py#L389), [fiber_simulation.py:553](../myogen/simulator/core/emg/fiber_simulation.py#L553), [fiber_simulation.py:568](../myogen/simulator/core/emg/fiber_simulation.py#L568), [bioelectric.py:8](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L8), [bioelectric.py:13](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L13), [bioelectric.py:25](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L25), [bioelectric.py:48](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L48), [bioelectric.py:76](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L76), [simulate_fiber.py:381](../myogen/simulator/core/emg/surface/simulate_fiber.py#L381), [simulate_fiber.py:383](../myogen/simulator/core/emg/surface/simulate_fiber.py#L383), [simulate_fiber.py:461](../myogen/simulator/core/emg/surface/simulate_fiber.py#L461), [simulate_fiber.py:478](../myogen/simulator/core/emg/surface/simulate_fiber.py#L478), [simulate_fiber.py:783](../myogen/simulator/core/emg/surface/simulate_fiber.py#L783), [simulate_fiber.py:785](../myogen/simulator/core/emg/surface/simulate_fiber.py#L785), [surface_emg.py:104](../myogen/simulator/core/emg/surface/surface_emg.py#L104), [surface_emg.py:106](../myogen/simulator/core/emg/surface/surface_emg.py#L106), [surface_emg.py:579](../myogen/simulator/core/emg/surface/surface_emg.py#L579), [test_iemg_sfap_units.py:4](../tests/test_iemg_sfap_units.py#L4), [test_iemg_sfap_units.py:5](../tests/test_iemg_sfap_units.py#L5) |
 | <a id="dimitrova2006"></a>**dimitrova2006**: Dimitrova NA, Dimitrov GV (2006). Electromyography (EMG) modeling. Wiley Encyclopedia of Biomedical Engineering. | FT | CONTEXT | IAP profile length = duration × CV; spectral content scales with IAP duration (MDF gap). | — |
 | <a id="dimitrov1998"></a>**dimitrov1998**: Dimitrov GV, Dimitrova NA et al. (1998). Three papers, Med Eng Phys 20:374, 580, 702 (titles unverified). | FT | CONTEXT | Spectra vs IAP duration and CV. No measured human IAP duration found. | — |
 | <a id="rodriguezfalces2013"></a>**rodriguezfalces2013**: Rodriguez-Falces J et al. (2013). Med Biol Eng Comput 50:447 (IAP / dipole metrics; title unverified). [doi:10.1007/s11517-012-0879-7](https://doi.org/10.1007/s11517-012-0879-7) | FT | CONTEXT | IAP duration metrics for the λ question. | — |
@@ -75,34 +75,34 @@ To add a source: add it to the table in `scripts/build_literature.py`, link it f
 | <a id="kouyoumdjian2023"></a>**kouyoumdjian2023**: Kouyoumdjian JA et al. (2023). Muscle fiber conduction velocity in situ revisited. Front Neurol 14:1118510. [doi:10.3389/fneur.2023.1118510](https://doi.org/10.3389/fneur.2023.1118510) · [PMC9996126](https://pmc.ncbi.nlm.nih.gov/articles/PMC9996126/) | FT | CONTEXT | TA single-fibre CV 4.10 ± 0.66 m/s. | — |
 | <a id="zalewska2012"></a>**zalewska2012**: Zalewska E et al. (2012). A method for determination of muscle fiber diameter using single fiber potential (SFP) analysis. Med Biol Eng Comput. [doi:10.1007/s11517-012-0965-x](https://doi.org/10.1007/s11517-012-0965-x) · [PMC3514689](https://pmc.ncbi.nlm.nih.gov/articles/PMC3514689/) | FT | CONTEXT | Alternative CV–diameter line v = 2.2 + 0.05(d − 25). | — |
 | <a id="blijham2006"></a>**blijham2006**: Blijham PJ, ter Laak HJ, Schelhaas HJ, van Engelen BGM, Stegeman DF, Zwarts MJ (2006). Relation between muscle fiber conduction velocity and fiber size in neuromuscular disorders. J Appl Physiol 100:1837–1841. [doi:10.1152/japplphysiol.01009.2005](https://doi.org/10.1152/japplphysiol.01009.2005) | FT | CONTEXT | Read 2026-11-02. Invasive MFCV (direct stimulation, biceps) against needle-biopsy diameters (mostly quadriceps; 53 of 96 in the same muscle) in 96 patients, 24 of whom had no neuromuscular disorder (pp. 1837–1838). Only extremes were paired: smallest fibre with slowest MFCV, largest with fastest; R 0.58 / 0.46 (p. 1838 Table 1). Overall fit MFCV = 0.043·d[µm] + 0.83 m/s (p. 1840, Fig. 5), 'remarkably similar' to nandedkar1983 (p. 1840). Control ranges biceps 30–70, quadriceps 50–95 µm, citing polgar1973 (p. 1840). Cited with nandedkar1983 by delvecchio2018 (p. 7, ref. 6). Supports the slope and intercept of the Nandedkar line; patient data, so not adopted as its replacement. Cited in code for the equal relative spread (within 2 %) of the two lines and for the quadriceps range against the removed 80 µm clip. | —<br>Cited at: [fiber_dispersion.py:92](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L92) |
-| <a id="delvecchio2018"></a>**delvecchio2018**: Del Vecchio A, Negro F, Felici F, Farina D (2018). Distribution of muscle fibre conduction velocity for representative samples of motor units in the full recruitment range of the tibialis anterior muscle. Acta Physiol 222:e12930. [doi:10.1111/apha.12930](https://doi.org/10.1111/apha.12930) | FT | ADOPTED | TA, 10 subjects, 406 MUs: MU CV 2.78–6.21 m/s (4.32 ± 0.71), linear in recruitment threshold (R² 0.70 ± 0.09); 3.89 ± 0.50 m/s at RT 0–30 % MVC, 5.12 ± 0.45 at 50–70 % (p. 2, Fig. 1 p. 3). HD-sEMG grid of 64 electrodes. Since 2026-11-02 (Devon) 4.32 ± 0.71 sets the Muscle CV defaults (intramuscular EMG; was 3.89 ± 0.50 before fibre typing), matched over unit means since 2026-11-03 (was over fibres); R² 0.70 gives the 0.39 m/s residual scatter of cv_from_recruitment_force (derived); the linear CV-threshold relation is why fibre typing uses one diameter ladder across the pool. Within-unit fibre CV distribution unknown (p. 5); CoV 5 % is a modelling assumption (Devon). CV→diameter by MFD = (MUCV − 0.95)/0.05 (p. 7 eq. 1). Basis of fiber_cv_threshold_slope (0.63 × mean, derived). Earlier cited in error as the 2017 J Appl Physiol paper. | [SurfaceFiberDispersion](../myogen/simulator/core/emg/surface/fiber_dispersion.py), [CHANGELOG](../CHANGELOG.md), [Muscle](../myogen/simulator/core/muscle/muscle.py)<br>Cited at: [fiber_dispersion.py:64](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L64), [fiber_dispersion.py:75](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L75), [fiber_dispersion.py:91](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L91), [fiber_dispersion.py:184](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L184), [fiber_dispersion.py:204](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L204), [surface_emg.py:111](../myogen/simulator/core/emg/surface/surface_emg.py#L111), [fiber_types.py:24](../myogen/simulator/core/muscle/fiber_types.py#L24), [muscle.py:35](../myogen/simulator/core/muscle/muscle.py#L35), [muscle.py:184](../myogen/simulator/core/muscle/muscle.py#L184), [muscle.py:197](../myogen/simulator/core/muscle/muscle.py#L197), [muscle.py:592](../myogen/simulator/core/muscle/muscle.py#L592), [ramp_analysis.py:144](../myogen/utils/ramp_analysis.py#L144), [ramp_analysis.py:152](../myogen/utils/ramp_analysis.py#L152), [test_muscle_fibre_properties.py:138](../tests/test_muscle_fibre_properties.py#L138), [test_ramp_recruitment_force.py:43](../tests/test_ramp_recruitment_force.py#L43) |
+| <a id="delvecchio2018"></a>**delvecchio2018**: Del Vecchio A, Negro F, Felici F, Farina D (2018). Distribution of muscle fibre conduction velocity for representative samples of motor units in the full recruitment range of the tibialis anterior muscle. Acta Physiol 222:e12930. [doi:10.1111/apha.12930](https://doi.org/10.1111/apha.12930) | FT | ADOPTED | TA, 10 subjects, 406 MUs: MU CV 2.78–6.21 m/s (4.32 ± 0.71), linear in recruitment threshold (R² 0.70 ± 0.09); 3.89 ± 0.50 m/s at RT 0–30 % MVC, 5.12 ± 0.45 at 50–70 % (p. 2, Fig. 1 p. 3). HD-sEMG grid of 64 electrodes. Since 2026-11-02 (Devon) 4.32 ± 0.71 sets the Muscle CV defaults (intramuscular EMG; was 3.89 ± 0.50 before fibre typing), matched over unit means since 2026-11-03 (was over fibres); R² 0.70 gives the 0.39 m/s residual scatter of cv_from_recruitment_force (derived); the linear CV-threshold relation is why fibre typing uses one diameter ladder across the pool. Within-unit fibre CV distribution unknown (p. 5); CoV 5 % is a modelling assumption (Devon). CV→diameter by MFD = (MUCV − 0.95)/0.05 (p. 7 eq. 1). Basis of fiber_cv_threshold_slope (0.63 × mean, derived). Earlier cited in error as the 2017 J Appl Physiol paper. | [SurfaceFiberDispersion](../myogen/simulator/core/emg/surface/fiber_dispersion.py), [CHANGELOG](../CHANGELOG.md), [Muscle](../myogen/simulator/core/muscle/muscle.py)<br>Cited at: [fdi_semg.yaml:58](../myogen/profiles/fdi_semg.yaml#L58), [fiber_dispersion.py:64](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L64), [fiber_dispersion.py:75](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L75), [fiber_dispersion.py:91](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L91), [fiber_dispersion.py:184](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L184), [fiber_dispersion.py:204](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L204), [surface_emg.py:111](../myogen/simulator/core/emg/surface/surface_emg.py#L111), [fiber_types.py:24](../myogen/simulator/core/muscle/fiber_types.py#L24), [muscle.py:35](../myogen/simulator/core/muscle/muscle.py#L35), [muscle.py:184](../myogen/simulator/core/muscle/muscle.py#L184), [muscle.py:197](../myogen/simulator/core/muscle/muscle.py#L197), [muscle.py:592](../myogen/simulator/core/muscle/muscle.py#L592), [ramp_analysis.py:144](../myogen/utils/ramp_analysis.py#L144), [ramp_analysis.py:152](../myogen/utils/ramp_analysis.py#L152), [test_muscle_fibre_properties.py:138](../tests/test_muscle_fibre_properties.py#L138), [test_ramp_recruitment_force.py:43](../tests/test_ramp_recruitment_force.py#L43) |
 | <a id="delvecchio2017"></a>**delvecchio2017**: Del Vecchio A, Negro F, Felici F, Farina D (2017). Associations between motor unit action potential parameters and surface EMG features. J Appl Physiol 123:835–843. [doi:10.1152/japplphysiol.00482.2017](https://doi.org/10.1152/japplphysiol.00482.2017) | FT | CONTEXT | TA, 13 subjects, ramps to 70 % MVC: MU CV linear in recruitment threshold (R² 0.64 ± 0.14); intercept 3.88 ± 0.35 m/s (Table 2), mean slope 0.0224 m/s per % MVC of threshold (mean of the 13 subject slopes in Table 3). Global MDF poorly related to recruitment. Not cited in code. | — |
-| <a id="akhmadeev2020"></a>**akhmadeev2020**: Akhmadeev K, Yu T, Le Carpentier E, Aoustin Y, Farina D (2020). Simulation of motor unit action potential recordings from intramuscular multichannel scanning electrodes. IEEE Trans Biomed Eng 67(7):2005–2014. Cited in MyoGen as 'Konstantin et al. 2020' (Konstantin is the first author's given name). [doi:10.1109/TBME.2019.2953680](https://doi.org/10.1109/TBME.2019.2953680) · [code](https://github.com/smetanadvorak/iemg_simulator) | FT + CODE | ADOPTED | Two-layer end-plate / nerve-branch model (shared iEMG/sEMG). Published eq. 5–6 coefficients are the default; reference-code coefficients kept as nmj_coefficients='reference_code'. Branch CVs 5/2 m/s from the reference code (paper: 10/1). | [nerve_branches](../myogen/simulator/core/emg/nerve_branches.py), [IntramuscularEMG](../myogen/simulator/core/emg/intramuscular/intramuscular_emg.py), [SurfaceFiberDispersion](../myogen/simulator/core/emg/surface/fiber_dispersion.py), [test_surface_dispersion](../tests/test_surface_dispersion.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [intramuscular_emg.py:70](../myogen/simulator/core/emg/intramuscular/intramuscular_emg.py#L70), [intramuscular_emg.py:73](../myogen/simulator/core/emg/intramuscular/intramuscular_emg.py#L73), [intramuscular_emg.py:76](../myogen/simulator/core/emg/intramuscular/intramuscular_emg.py#L76), [intramuscular_emg.py:79](../myogen/simulator/core/emg/intramuscular/intramuscular_emg.py#L79), [intramuscular_emg.py:100](../myogen/simulator/core/emg/intramuscular/intramuscular_emg.py#L100), [intramuscular_emg.py:280](../myogen/simulator/core/emg/intramuscular/intramuscular_emg.py#L280), [motor_unit_sim.py:8](../myogen/simulator/core/emg/intramuscular/motor_unit_sim.py#L8), [nerve_branches.py:1](../myogen/simulator/core/emg/nerve_branches.py#L1), [nerve_branches.py:21](../myogen/simulator/core/emg/nerve_branches.py#L21), [nerve_branches.py:22](../myogen/simulator/core/emg/nerve_branches.py#L22), [fiber_dispersion.py:44](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L44), [fiber_dispersion.py:53](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L53), [fiber_dispersion.py:179](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L179), [fiber_dispersion.py:181](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L181), [muscle.py:145](../myogen/simulator/core/muscle/muscle.py#L145), [muscle.py:296](../myogen/simulator/core/muscle/muscle.py#L296), [physiological_distribution.py:37](../myogen/simulator/core/physiological_distribution.py#L37), [physiological_distribution.py:75](../myogen/simulator/core/physiological_distribution.py#L75), [physiological_distribution.py:91](../myogen/simulator/core/physiological_distribution.py#L91) |
+| <a id="akhmadeev2020"></a>**akhmadeev2020**: Akhmadeev K, Yu T, Le Carpentier E, Aoustin Y, Farina D (2020). Simulation of motor unit action potential recordings from intramuscular multichannel scanning electrodes. IEEE Trans Biomed Eng 67(7):2005–2014. Cited in MyoGen as 'Konstantin et al. 2020' (Konstantin is the first author's given name). [doi:10.1109/TBME.2019.2953680](https://doi.org/10.1109/TBME.2019.2953680) · [code](https://github.com/smetanadvorak/iemg_simulator) | FT + CODE | ADOPTED | Two-layer end-plate / nerve-branch model (shared iEMG/sEMG). Published eq. 5–6 coefficients are the default; reference-code coefficients kept as nmj_coefficients='reference_code'. Branch CVs 5/2 m/s from the reference code (paper: 10/1). | [nerve_branches](../myogen/simulator/core/emg/nerve_branches.py), [IntramuscularEMG](../myogen/simulator/core/emg/intramuscular/intramuscular_emg.py), [SurfaceFiberDispersion](../myogen/simulator/core/emg/surface/fiber_dispersion.py), [test_surface_dispersion](../tests/test_surface_dispersion.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [fdi_semg.yaml:60](../myogen/profiles/fdi_semg.yaml#L60), [intramuscular_emg.py:70](../myogen/simulator/core/emg/intramuscular/intramuscular_emg.py#L70), [intramuscular_emg.py:73](../myogen/simulator/core/emg/intramuscular/intramuscular_emg.py#L73), [intramuscular_emg.py:76](../myogen/simulator/core/emg/intramuscular/intramuscular_emg.py#L76), [intramuscular_emg.py:79](../myogen/simulator/core/emg/intramuscular/intramuscular_emg.py#L79), [intramuscular_emg.py:100](../myogen/simulator/core/emg/intramuscular/intramuscular_emg.py#L100), [intramuscular_emg.py:280](../myogen/simulator/core/emg/intramuscular/intramuscular_emg.py#L280), [motor_unit_sim.py:8](../myogen/simulator/core/emg/intramuscular/motor_unit_sim.py#L8), [nerve_branches.py:1](../myogen/simulator/core/emg/nerve_branches.py#L1), [nerve_branches.py:21](../myogen/simulator/core/emg/nerve_branches.py#L21), [nerve_branches.py:22](../myogen/simulator/core/emg/nerve_branches.py#L22), [fiber_dispersion.py:44](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L44), [fiber_dispersion.py:53](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L53), [fiber_dispersion.py:179](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L179), [fiber_dispersion.py:181](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L181), [muscle.py:145](../myogen/simulator/core/muscle/muscle.py#L145), [muscle.py:296](../myogen/simulator/core/muscle/muscle.py#L296), [physiological_distribution.py:37](../myogen/simulator/core/physiological_distribution.py#L37), [physiological_distribution.py:75](../myogen/simulator/core/physiological_distribution.py#L75), [physiological_distribution.py:91](../myogen/simulator/core/physiological_distribution.py#L91) |
 | <a id="stalbergkarlsson2001"></a>**stalbergkarlsson2001**: Stålberg E, Karlsson L (2001). Simulation of the normal concentric needle electromyogram by using a muscle model. Clin Neurophysiol 112:464–471. [doi:10.1016/S1388-2457(01)00459-X](https://doi.org/10.1016/S1388-2457(01)00459-X) | FT | CONTEXT | Bound: end-plate position 0 ± 1 mm; per-fibre delay 0–1 ms. Piasecki 2020 reuses this value (one chain, not two anchors). Companion paper: Simulation of EMG in pathological situations, Clin Neurophysiol 112:869–878, doi:10.1016/S1388-2457(01)00498-9 (title and pages were mixed up here before 2026-11-03; checked against Zotero). Their fibre diameter 50 ± 5 µm within each motor unit is a model setting, no measurement cited (literature check 2026-11-03). | — |
 | <a id="piasecki2020"></a>**piasecki2020**: Piasecki M et al. (2020/2021). Near-fibre EMG (medRxiv; Clin Neurophysiol 132:1089). [doi:10.1101/2020.10.07.20208348](https://doi.org/10.1101/2020.10.07.20208348) | FT | CONTEXT | Healthy end-plate SD 1 mm (from Stålberg & Karlsson). | — |
 | <a id="lateva2010"></a>**lateva2010**: Lateva ZC, McGill KC, Johanson ME (2010). J Appl Physiol 108:1530. [PMC2886675](https://pmc.ncbi.nlm.nih.gov/articles/PMC2886675/) | FT | CONTEXT | Myelinated branch ~50 m/s bounds fast branch delays; fibre-end spread. | — |
 | <a id="keenan2005"></a>**keenan2005**: Keenan KG, Farina D et al. (2005). J Appl Physiol 98:120 (FDI sEMG simulation). [doi:10.1152/japplphysiol.00894.2004](https://doi.org/10.1152/japplphysiol.00894.2004) | FT | ADOPTED | Convention for the opt-in 5 mm tendon-end band (no measurement of fibre-end stagger found). | [SurfaceFiberDispersion](../myogen/simulator/core/emg/surface/fiber_dispersion.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [fiber_dispersion.py:187](../myogen/simulator/core/emg/surface/fiber_dispersion.py#L187) |
 | <a id="farinafosci2002"></a>**farinafosci2002**: Farina D, Fosci M, Merletti R (2002). J Appl Physiol 92:235. [doi:10.1152/jappl.2002.92.1.235](https://doi.org/10.1152/jappl.2002.92.1.235) | FT | CONTEXT | IZ / tendon spread 10 mm used as the sensitivity bound. | — |
 | <a id="jahanmiri2015"></a>**jahanmiri2015**: Jahanmiri-Nezhad F et al. (2015). Front Hum Neurosci 9:239. [doi:10.3389/fnhum.2015.00239](https://doi.org/10.3389/fnhum.2015.00239) | FT | CONTEXT | Single-MU IZ spread 9–18 mm. | — |
-| <a id="martin1988"></a>**martin1988**: Martin TP, Bodine-Fowler S, Roy RR, Eldred E, Edgerton VR (1988). Metabolic and fiber size properties of cat tibialis anterior motor units. Am J Physiol 255:C43–C50. [doi:10.1152/ajpcell.1988.255.1.C43](https://doi.org/10.1152/ajpcell.1988.255.1.C43) | FT | CONTEXT | Cat TA, 14 units by glycogen depletion, 8–10 unit fibres each, single 20-µm sections (p. C44). Within-unit fibre cross-section CoV per unit (Table 1 p. C46): type S 14/21/15/16 %, FR 10/18/19/34 %, FI 33/20 %, FF 28/30/31/32 %; mean 25 % against 46 % for non-unit fibres (p. C45–C46); Bodine et al. 1987 found 28 % vs 38 % (p. C46). Within-unit variability approximates that of the unit type (p. C49). Only within-unit fibre-size data found (read 2026-11-03). Diameter CoV ≈ area CoV / 2: S ≈ 0.08, fast ≈ 0.13. Muscle default within-unit diameter log-SD is 0.166 (from type_log_sd 0.2, an assumption), above every unit here. | —<br>Cited at: [surface_emg.py:115](../myogen/simulator/core/emg/surface/surface_emg.py#L115), [fiber_types.py:15](../myogen/simulator/core/muscle/fiber_types.py#L15), [fiber_types.py:176](../myogen/simulator/core/muscle/fiber_types.py#L176), [muscle.py:215](../myogen/simulator/core/muscle/muscle.py#L215), [test_fiber_types.py:127](../tests/test_fiber_types.py#L127), [test_muscle_fibre_properties.py:131](../tests/test_muscle_fibre_properties.py#L131) |
+| <a id="martin1988"></a>**martin1988**: Martin TP, Bodine-Fowler S, Roy RR, Eldred E, Edgerton VR (1988). Metabolic and fiber size properties of cat tibialis anterior motor units. Am J Physiol 255:C43–C50. [doi:10.1152/ajpcell.1988.255.1.C43](https://doi.org/10.1152/ajpcell.1988.255.1.C43) | FT | CONTEXT | Cat TA, 14 units by glycogen depletion, 8–10 unit fibres each, single 20-µm sections (p. C44). Within-unit fibre cross-section CoV per unit (Table 1 p. C46): type S 14/21/15/16 %, FR 10/18/19/34 %, FI 33/20 %, FF 28/30/31/32 %; mean 25 % against 46 % for non-unit fibres (p. C45–C46); Bodine et al. 1987 found 28 % vs 38 % (p. C46). Within-unit variability approximates that of the unit type (p. C49). Only within-unit fibre-size data found (read 2026-11-03). Diameter CoV ≈ area CoV / 2: S ≈ 0.08, fast ≈ 0.13. Muscle default within-unit diameter log-SD is 0.166 (from type_log_sd 0.2, an assumption), above every unit here. | —<br>Cited at: [fdi_semg.yaml:57](../myogen/profiles/fdi_semg.yaml#L57), [surface_emg.py:115](../myogen/simulator/core/emg/surface/surface_emg.py#L115), [fiber_types.py:15](../myogen/simulator/core/muscle/fiber_types.py#L15), [fiber_types.py:176](../myogen/simulator/core/muscle/fiber_types.py#L176), [muscle.py:215](../myogen/simulator/core/muscle/muscle.py#L215), [test_fiber_types.py:127](../tests/test_fiber_types.py#L127), [test_muscle_fibre_properties.py:131](../tests/test_muscle_fibre_properties.py#L131) |
 | <a id="andreassen1987"></a>**andreassen1987**: Andreassen S, Arendt-Nielsen L (1987). Muscle fibre conduction velocity in motor units of the human anterior tibial muscle: a new size principle parameter. J Physiol 391:561–571. [doi:10.1113/jphysiol.1987.sp016756](https://doi.org/10.1113/jphysiol.1987.sp016756) | FT | CONTEXT | Human TA, 33 units by intramuscular axon microstimulation, surface tripolar array 15 mm: unit CV 2.6–5.3 m/s, 3.7 ± 0.7 (p. 567, Fig. 3D p. 566), a between-unit spread matching delvecchio2018 (0.71). Unit CV correlates with twitch torque r = 0.87 (p. 568 Fig. 4B). Argues (p. 569) that a wide unit-CV range needs fibres of one unit to share about the same CV and diameter; not measured. Twitch-torque ratio for FDI ~100 / 40 / 60–120 from other studies (Table 1 p. 565). | —<br>Cited at: [muscle.py:185](../myogen/simulator/core/muscle/muscle.py#L185), [muscle.py:592](../myogen/simulator/core/muscle/muscle.py#L592), [test_muscle_fibre_properties.py:139](../tests/test_muscle_fibre_properties.py#L139) |
 
 ## C. Muscle anatomy (FDI dimensions, layers, motor-unit counts)
 
 | Source | Status | Use | What it supports | Used in |
 |---|---|---|---|---|
-| <a id="jacobson1992"></a>**jacobson1992**: Jacobson MD, Raab R, Fazeli BM, Abrams RA, Botte MJ, Lieber RL (1992). Architectural design of the human intrinsic hand muscles. J Hand Surg Am 17(5):804–809. [doi:10.1016/0363-5023(92)90446-V](https://doi.org/10.1016/0363-5023(92)90446-V) | ABS | ADOPTED | FDI fibre length 31.7 mm and var_fiber_length__mm 2.8 mm (used as a uniform half-width; probably a between-specimen spread). **Table 1 unverified.** | [Muscle](../myogen/simulator/core/muscle/muscle.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [muscle.py:156](../myogen/simulator/core/muscle/muscle.py#L156), [muscle.py:219](../myogen/simulator/core/muscle/muscle.py#L219), [muscle.py:223](../myogen/simulator/core/muscle/muscle.py#L223), [muscle.py:297](../myogen/simulator/core/muscle/muscle.py#L297) |
-| <a id="li2018"></a>**li2018**: Li H, Furst DE, Jin H et al. (2018). High-frequency ultrasound of the skin in systemic sclerosis. Arthritis Res Ther 20:181. [doi:10.1186/s13075-018-1686-9](https://doi.org/10.1186/s13075-018-1686-9) | FT | ADOPTED | Skin thickness 1.2 mm: hand dorsum over MCP II–III, median of 31 controls (0.79–1.50). | [Muscle](../myogen/simulator/core/muscle/muscle.py), [test_surface_amplitude](../tests/test_surface_amplitude.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [muscle.py:234](../myogen/simulator/core/muscle/muscle.py#L234), [muscle.py:305](../myogen/simulator/core/muscle/muscle.py#L305) |
-| <a id="nishikawa2024"></a>**nishikawa2024**: Nishikawa Y, Watanabe K, Holobar A et al. (2024). Sex differences in laterality of motor unit firing behavior of the first dorsal interosseous muscle. Eur J Appl Physiol 124:1979–1990. [doi:10.1007/s00421-024-05420-7](https://doi.org/10.1007/s00421-024-05420-7) | FT | ADOPTED | FDI depth to muscle 2.54–3.30 mm; fat 1.3 mm = 2.5 − 1.2 (derived from the low end). | [Muscle](../myogen/simulator/core/muscle/muscle.py), [test_surface_amplitude](../tests/test_surface_amplitude.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [muscle.py:230](../myogen/simulator/core/muscle/muscle.py#L230), [muscle.py:306](../myogen/simulator/core/muscle/muscle.py#L306) |
+| <a id="jacobson1992"></a>**jacobson1992**: Jacobson MD, Raab R, Fazeli BM, Abrams RA, Botte MJ, Lieber RL (1992). Architectural design of the human intrinsic hand muscles. J Hand Surg Am 17(5):804–809. [doi:10.1016/0363-5023(92)90446-V](https://doi.org/10.1016/0363-5023(92)90446-V) | ABS | ADOPTED | FDI fibre length 31.7 mm and var_fiber_length__mm 2.8 mm (used as a uniform half-width; probably a between-specimen spread). **Table 1 unverified.** | [Muscle](../myogen/simulator/core/muscle/muscle.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [fdi_semg.yaml:51](../myogen/profiles/fdi_semg.yaml#L51), [muscle.py:156](../myogen/simulator/core/muscle/muscle.py#L156), [muscle.py:219](../myogen/simulator/core/muscle/muscle.py#L219), [muscle.py:223](../myogen/simulator/core/muscle/muscle.py#L223), [muscle.py:297](../myogen/simulator/core/muscle/muscle.py#L297) |
+| <a id="li2018"></a>**li2018**: Li H, Furst DE, Jin H et al. (2018). High-frequency ultrasound of the skin in systemic sclerosis. Arthritis Res Ther 20:181. [doi:10.1186/s13075-018-1686-9](https://doi.org/10.1186/s13075-018-1686-9) | FT | ADOPTED | Skin thickness 1.2 mm: hand dorsum over MCP II–III, median of 31 controls (0.79–1.50). | [Muscle](../myogen/simulator/core/muscle/muscle.py), [test_surface_amplitude](../tests/test_surface_amplitude.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [fdi_semg.yaml:53](../myogen/profiles/fdi_semg.yaml#L53), [muscle.py:234](../myogen/simulator/core/muscle/muscle.py#L234), [muscle.py:305](../myogen/simulator/core/muscle/muscle.py#L305) |
+| <a id="nishikawa2024"></a>**nishikawa2024**: Nishikawa Y, Watanabe K, Holobar A et al. (2024). Sex differences in laterality of motor unit firing behavior of the first dorsal interosseous muscle. Eur J Appl Physiol 124:1979–1990. [doi:10.1007/s00421-024-05420-7](https://doi.org/10.1007/s00421-024-05420-7) | FT | ADOPTED | FDI depth to muscle 2.54–3.30 mm; fat 1.3 mm = 2.5 − 1.2 (derived from the low end). | [Muscle](../myogen/simulator/core/muscle/muscle.py), [test_surface_amplitude](../tests/test_surface_amplitude.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [fdi_semg.yaml:54](../myogen/profiles/fdi_semg.yaml#L54), [muscle.py:230](../myogen/simulator/core/muscle/muscle.py#L230), [muscle.py:306](../myogen/simulator/core/muscle/muscle.py#L306) |
 | <a id="storchle2018"></a>**storchle2018**: Störchle P, Müller W, Sengeis M et al. (2018). Measurement of mean subcutaneous fat thickness. Sci Rep 8:16268. [doi:10.1038/s41598-018-34213-0](https://doi.org/10.1038/s41598-018-34213-0) | FT | REJECTED | Old fat default 0.3 mm (whole-hand average, excludes skin). | [Muscle](../myogen/simulator/core/muscle/muscle.py)<br>Cited at: [muscle.py:231](../myogen/simulator/core/muscle/muscle.py#L231), [muscle.py:299](../myogen/simulator/core/muscle/muscle.py#L299) |
 | <a id="brodar1960"></a>**brodar1960**: Brodar V (1960). Observations on skin thickness and subcutaneous tissue in man. Z Morphol Anthropol 50:386–395. | SEC | REJECTED | Old skin default 1.29 mm (whole-body male average). | [Muscle](../myogen/simulator/core/muscle/muscle.py)<br>Cited at: [muscle.py:235](../myogen/simulator/core/muscle/muscle.py#L235), [muscle.py:300](../myogen/simulator/core/muscle/muscle.py#L300) |
 | <a id="krammer2025"></a>**krammer2025**: Krammer et al. (2025). Diagnostics 15:1600. | FT | CONTEXT | Hand-dorsum dermis 0.52 ± 0.16 mm. | — |
 | <a id="infantolino2010"></a>**infantolino2010**: Infantolino BW, Challis JH (2010). Architectural properties of the first dorsal interosseous muscle. J Anat 216:463. [doi:10.1111/j.1469-7580.2009.01196.x](https://doi.org/10.1111/j.1469-7580.2009.01196.x) · [PMC2849523](https://pmc.ncbi.nlm.nih.gov/articles/PMC2849523/) | FT | CONTEXT | Volume 5.2 mL, fascicle 29.9 mm; derived CSA ~89 mm² vs model 150 mm². | — |
 | <a id="burketsairis1973"></a>**burketsairis1973**: Burke RE, Tsairis P (1973). Anatomy and innervation ratios in motor units of cat gastrocnemius. J Physiol 234:749–765. [doi:10.1113/jphysiol.1973.sp010370](https://doi.org/10.1113/jphysiol.1973.sp010370) | FT | CONTEXT | Cat MG: innervation ratios FF ~550–750, FR ~400–550 (p. 756; Table 2 p. 758); unit fibres spread over a wide territory, 1–2 per primary fascicle, 3–5 per 100 fibres (p. 753); mean unit fibre areas FF 5290, FR 2890, S 1734 µm² (p. 760). No within-unit fibre-area spread given. | — |
-| <a id="polgar1973"></a>**polgar1973**: Polgar J, Johnson MA, Weightman D, Appleton D (1973). Data on fibre size in thirty-six human muscles: an autopsy study. J Neurol Sci 19:307–318. [doi:10.1016/0022-510X(73)90094-4](https://doi.org/10.1016/0022-510X(73)90094-4) | FT | ADOPTED | FDI fibre diameter defaults (since 2026-11-02, Devon). FDI (autopsy): type I 50.0 µm, type II 67.0 µm [p. 310 Table 2]; type I fibres 57.4 % [p. 312 Table 4]. Weighted by that proportion: mean 57.2 µm (Muscle since 2026-11-02: fibre types by motor unit on one diameter ladder that keeps these type means, log-SD 0.2 per type an assumption; was a per-fibre two-type mixture earlier the same day), root-mean-square 57.9 µm (SurfaceEMG source default, since the source scales with d²). Replaces the 50 µm attributed to brooke1970. Type II − type I 17.02 µm, P < 0.001 (p. 312 Table 4). FDI 95 % limits for a single subject's mean (between post-mortems, not the fibre-to-fibre spread): type I 35–66, type II 47–86 µm (p. 313 Table 5; methods p. 309). Minor-axis diameters, ≥ 200 fibres per sample (p. 308). Gives no within-muscle fibre SD. | —<br>Cited at: [surface_emg.py:95](../myogen/simulator/core/emg/surface/surface_emg.py#L95), [surface_emg.py:95](../myogen/simulator/core/emg/surface/surface_emg.py#L95), [fiber_types.py:7](../myogen/simulator/core/muscle/fiber_types.py#L7), [fiber_types.py:103](../myogen/simulator/core/muscle/fiber_types.py#L103), [fiber_types.py:105](../myogen/simulator/core/muscle/fiber_types.py#L105), [muscle.py:201](../myogen/simulator/core/muscle/muscle.py#L201), [muscle.py:206](../myogen/simulator/core/muscle/muscle.py#L206), [muscle.py:554](../myogen/simulator/core/muscle/muscle.py#L554), [muscle.py:555](../myogen/simulator/core/muscle/muscle.py#L555), [test_muscle_fibre_properties.py:57](../tests/test_muscle_fibre_properties.py#L57), [test_muscle_fibre_properties.py:278](../tests/test_muscle_fibre_properties.py#L278), [test_muscle_fibre_properties.py:278](../tests/test_muscle_fibre_properties.py#L278), [test_surface_amplitude.py:26](../tests/test_surface_amplitude.py#L26) |
+| <a id="polgar1973"></a>**polgar1973**: Polgar J, Johnson MA, Weightman D, Appleton D (1973). Data on fibre size in thirty-six human muscles: an autopsy study. J Neurol Sci 19:307–318. [doi:10.1016/0022-510X(73)90094-4](https://doi.org/10.1016/0022-510X(73)90094-4) | FT | ADOPTED | FDI fibre diameter defaults (since 2026-11-02, Devon). FDI (autopsy): type I 50.0 µm, type II 67.0 µm [p. 310 Table 2]; type I fibres 57.4 % [p. 312 Table 4]. Weighted by that proportion: mean 57.2 µm (Muscle since 2026-11-02: fibre types by motor unit on one diameter ladder that keeps these type means, log-SD 0.2 per type an assumption; was a per-fibre two-type mixture earlier the same day), root-mean-square 57.9 µm (SurfaceEMG source default, since the source scales with d²). Replaces the 50 µm attributed to brooke1970. Type II − type I 17.02 µm, P < 0.001 (p. 312 Table 4). FDI 95 % limits for a single subject's mean (between post-mortems, not the fibre-to-fibre spread): type I 35–66, type II 47–86 µm (p. 313 Table 5; methods p. 309). Minor-axis diameters, ≥ 200 fibres per sample (p. 308). Gives no within-muscle fibre SD. | —<br>Cited at: [fdi_semg.yaml:56](../myogen/profiles/fdi_semg.yaml#L56), [surface_emg.py:95](../myogen/simulator/core/emg/surface/surface_emg.py#L95), [surface_emg.py:95](../myogen/simulator/core/emg/surface/surface_emg.py#L95), [fiber_types.py:7](../myogen/simulator/core/muscle/fiber_types.py#L7), [fiber_types.py:103](../myogen/simulator/core/muscle/fiber_types.py#L103), [fiber_types.py:105](../myogen/simulator/core/muscle/fiber_types.py#L105), [muscle.py:201](../myogen/simulator/core/muscle/muscle.py#L201), [muscle.py:206](../myogen/simulator/core/muscle/muscle.py#L206), [muscle.py:554](../myogen/simulator/core/muscle/muscle.py#L554), [muscle.py:555](../myogen/simulator/core/muscle/muscle.py#L555), [test_muscle_fibre_properties.py:57](../tests/test_muscle_fibre_properties.py#L57), [test_muscle_fibre_properties.py:278](../tests/test_muscle_fibre_properties.py#L278), [test_muscle_fibre_properties.py:278](../tests/test_muscle_fibre_properties.py#L278), [test_surface_amplitude.py:26](../tests/test_surface_amplitude.py#L26) |
 | <a id="garnett1979"></a>**garnett1979**: Garnett RAF, O'Donovan MJ, Stephens JA, Taylor A (1979). Motor unit organization of human medial gastrocnemius. J Physiol 287:33–43. [doi:10.1113/jphysiol.1979.sp012643](https://doi.org/10.1113/jphysiol.1979.sp012643) · [PMC1281479](https://pmc.ncbi.nlm.nih.gov/articles/PMC1281479/) | FT | CONTEXT | 57 human MG units (microstimulation, glycogen depletion, biopsy). S units small/slow/resistant, FR intermediate, FF large/fatigable [p. 33]; S units type 1 fibres, FF type 2b [p. 39–41]. Ranked by size, S and FR interleave at low force; only FF added above ~30 % cumulative twitch tension [p. 42]. Sample biased against largest/smallest units [p. 41]. Basis (with limits) for typing units by threshold (Muscle fibre typing, 2026-11-02, cited in fiber_types.py and Muscle docstring): strict type I-first order is a simplification. | [Muscle](../myogen/simulator/core/muscle/muscle.py)<br>Cited at: [fiber_types.py:9](../myogen/simulator/core/muscle/fiber_types.py#L9), [muscle.py:204](../myogen/simulator/core/muscle/muscle.py#L204) |
 | <a id="hwang2013"></a>**hwang2013**: Hwang K, Huan F, Kim DJ (2013). Muscle fibre types of the lumbrical, interossei, flexor, and extensor muscles moving the index finger. J Plast Surg Hand Surg 47:268–272. [doi:10.3109/2000656X.2012.755988](https://doi.org/10.3109/2000656X.2012.755988) | FT | REJECTED | FDI ~31 µm (elderly, paraffin, diameter undefined): judged an underestimate. Source of Botelho's 24.8–40.6 µm. | — |
-| <a id="feinstein1955"></a>**feinstein1955**: Feinstein B, Lindegård B, Nyman E, Wohlfart G (1955). Morphologic studies of motor units in normal human muscles. Acta Anat 23:127–142. | SEC | CONTEXT | FDI ~119 MUs, ~40,000 fibres, ~340 fibres/MU (via Piasecki 2016, De Luca 1982). | — |
+| <a id="feinstein1955"></a>**feinstein1955**: Feinstein B, Lindegård B, Nyman E, Wohlfart G (1955). Morphologic studies of motor units in normal human muscles. Acta Anat 23:127–142. | SEC | CONTEXT | FDI ~119 MUs, ~40,000 fibres, ~340 fibres/MU (via Piasecki 2016, De Luca 1982). | —<br>Cited at: [fdi_semg.yaml:29](../myogen/profiles/fdi_semg.yaml#L29), [fdi_semg.yaml:52](../myogen/profiles/fdi_semg.yaml#L52) |
 | <a id="zong2021"></a>**zong2021**: Zong Y et al. (2021). IEEE Trans Neural Syst Rehabil Eng 29:1199. [doi:10.1109/TNSRE.2021.3088061](https://doi.org/10.1109/TNSRE.2021.3088061) · [PMC8780215](https://pmc.ncbi.nlm.nih.gov/articles/PMC8780215/) | FT | CONTEXT | FDI MUNE 108 ± 21; CMAP 16.75 mV (ceiling). | — |
 | <a id="brooke1970"></a>**brooke1970**: Brooke MH, Kaiser KK (1970). Muscle fiber types: how many and what kind? Arch Neurol 23:369–379. [doi:10.1001/archneur.1970.00480280083010](https://doi.org/10.1001/archneur.1970.00480280083010) | LINK | CONTEXT | The code cited it for the log-normal fibre diameters (median 50 µm until 0.13.0) and for the 20–80 µm clip, removed 2026-11-02 (it held ~10 % of fibres at 80 µm). By its title a histochemical fibre-typing paper; the values are not checked against it. The diameter mean now comes from polgar1973. | — |
 
@@ -110,11 +110,11 @@ To add a source: add it to the table in `scripts/build_literature.py`, link it f
 
 | Source | Status | Use | What it supports | Used in |
 |---|---|---|---|---|
-| <a id="deluca1982"></a>**deluca1982**: De Luca CJ, LeFever RS, McCue MP, Xenakis AP (1982). Behaviour of human motor units in different muscles during linearly varying contractions. J Physiol 329:113–128. [doi:10.1113/jphysiol.1982.sp014293](https://doi.org/10.1113/jphysiol.1982.sp014293) | FT + FIG | ADOPTED | FDI rate-coding targets (onset 8.9 Hz, Table 1 peaks; fitted shape a modelling assumption) and the recruitment ceiling (~52 % MVC) behind --max-recruitment; Table 1 mean rates at 40/80 % MVC (25.3/41.4 Hz) are HIGH_FORCE_RATE_TARGETS (--high-force-rates). | [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py), [fr_optimization_targets](../myogen/utils/fr_optimization_targets.py)<br>Cited at: [06_optimize_dd_multilevel.py:78](../examples/02_finetune/06_optimize_dd_multilevel.py#L78), [06_optimize_dd_multilevel.py:555](../examples/02_finetune/06_optimize_dd_multilevel.py#L555), [06_optimize_dd_multilevel.py:1038](../examples/02_finetune/06_optimize_dd_multilevel.py#L1038) |
+| <a id="deluca1982"></a>**deluca1982**: De Luca CJ, LeFever RS, McCue MP, Xenakis AP (1982). Behaviour of human motor units in different muscles during linearly varying contractions. J Physiol 329:113–128. [doi:10.1113/jphysiol.1982.sp014293](https://doi.org/10.1113/jphysiol.1982.sp014293) | FT + FIG | ADOPTED | FDI rate-coding targets (onset 8.9 Hz, Table 1 peaks; fitted shape a modelling assumption) and the recruitment ceiling (~52 % MVC) behind --max-recruitment; Table 1 mean rates at 40/80 % MVC (25.3/41.4 Hz) are HIGH_FORCE_RATE_TARGETS (--high-force-rates). | [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py), [fr_optimization_targets](../myogen/utils/fr_optimization_targets.py)<br>Cited at: [06_optimize_dd_multilevel.py:77](../examples/02_finetune/06_optimize_dd_multilevel.py#L77), [06_optimize_dd_multilevel.py:354](../examples/02_finetune/06_optimize_dd_multilevel.py#L354), [06_optimize_dd_multilevel.py:837](../examples/02_finetune/06_optimize_dd_multilevel.py#L837) |
 | <a id="delucacontessa2012"></a>**delucacontessa2012**: De Luca CJ, Contessa P (2012). Hierarchical control of motor units in voluntary contractions. J Neurophysiol 107:178–195. | SEC | ADOPTED | Exponential onion-skin rate law used for the de_luca targets. | [fr_optimization_targets](../myogen/utils/fr_optimization_targets.py), [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py)<br>Cited at: [physiological_distribution.py:36](../myogen/simulator/core/physiological_distribution.py#L36), [physiological_distribution.py:74](../myogen/simulator/core/physiological_distribution.py#L74), [physiological_distribution.py:85](../myogen/simulator/core/physiological_distribution.py#L85), [fr_optimization_targets.py:15](../myogen/utils/fr_optimization_targets.py#L15), [fr_optimization_targets.py:108](../myogen/utils/fr_optimization_targets.py#L108), [fr_optimization_targets.py:330](../myogen/utils/fr_optimization_targets.py#L330) |
 | <a id="skarabot2023"></a>**skarabot2023**: Škarabot J et al. (2023). J Physiol (FDI/VL motor-unit decomposition). [PMC10952962](https://pmc.ncbi.nlm.nih.gov/articles/PMC10952962/) | FT | CONTEXT | FDI recruitment 0.8–59 % MVF (ceiling); VL peak rates in fr_optimization_targets. | [fr_optimization_targets](../myogen/utils/fr_optimization_targets.py)<br>Cited at: [fr_optimization_targets.py:138](../myogen/utils/fr_optimization_targets.py#L138), [fr_optimization_targets.py:193](../myogen/utils/fr_optimization_targets.py#L193), [fr_optimization_targets.py:202](../myogen/utils/fr_optimization_targets.py#L202), [fr_optimization_targets.py:209](../myogen/utils/fr_optimization_targets.py#L209), [fr_optimization_targets.py:212](../myogen/utils/fr_optimization_targets.py#L212), [fr_optimization_targets.py:218](../myogen/utils/fr_optimization_targets.py#L218) |
-| <a id="enokaduchateau2017"></a>**enokaduchateau2017**: Enoka RM, Duchateau J (2017). Rate coding and the control of muscle force. Cold Spring Harb Perspect Med 7:a029702. [doi:10.1101/cshperspect.a029702](https://doi.org/10.1101/cshperspect.a029702) | LINK | CONTEXT | Review: ISI CV 10–30 % in voluntary contractions, falling with force. With dideriksen2012, the basis of the optimiser's ISI CV target 15–22 % at 30 % MVC. Not read here. | —<br>Cited at: [06_optimize_dd_multilevel.py:61](../examples/02_finetune/06_optimize_dd_multilevel.py#L61) |
-| <a id="dideriksen2012"></a>**dideriksen2012**: Dideriksen JL, Negro F, Enoka RM, Farina D (2012). Motor unit recruitment strategies and muscle properties determine the influence of synaptic noise on force steadiness. J Neurophysiol 107:3357–3369. [doi:10.1152/jn.00938.2011](https://doi.org/10.1152/jn.00938.2011) | LINK | CONTEXT | Simulation of synaptic noise and force steadiness; cited with enokaduchateau2017 for the ISI CV target. Not read here. | —<br>Cited at: [06_optimize_dd_multilevel.py:61](../examples/02_finetune/06_optimize_dd_multilevel.py#L61) |
+| <a id="enokaduchateau2017"></a>**enokaduchateau2017**: Enoka RM, Duchateau J (2017). Rate coding and the control of muscle force. Cold Spring Harb Perspect Med 7:a029702. [doi:10.1101/cshperspect.a029702](https://doi.org/10.1101/cshperspect.a029702) | LINK | CONTEXT | Review: ISI CV 10–30 % in voluntary contractions, falling with force. With dideriksen2012, the basis of the optimiser's ISI CV target 15–22 % at 30 % MVC. Not read here. | —<br>Cited at: [06_optimize_dd_multilevel.py:60](../examples/02_finetune/06_optimize_dd_multilevel.py#L60) |
+| <a id="dideriksen2012"></a>**dideriksen2012**: Dideriksen JL, Negro F, Enoka RM, Farina D (2012). Motor unit recruitment strategies and muscle properties determine the influence of synaptic noise on force steadiness. J Neurophysiol 107:3357–3369. [doi:10.1152/jn.00938.2011](https://doi.org/10.1152/jn.00938.2011) | LINK | CONTEXT | Simulation of synaptic noise and force steadiness; cited with enokaduchateau2017 for the ISI CV target. Not read here. | —<br>Cited at: [06_optimize_dd_multilevel.py:60](../examples/02_finetune/06_optimize_dd_multilevel.py#L60) |
 | <a id="moritz2005"></a>**moritz2005**: Moritz CT, Barry BK, Pascoe MA, Enoka RM (2005). J Neurophysiol 93:2449 (FDI discharge variability; details from memory, unverified). | LINK | CONTEXT | ISI CV falls with rate — **quoted from memory, unverified**. | — |
 
 ## E. Motor-neuron model and spinal mechanics (synaptic input, PICs, inhibition, noise)
@@ -122,7 +122,7 @@ To add a source: add it to the table in `scripts/build_literature.py`, link it f
 | Source | Status | Use | What it supports | Used in |
 |---|---|---|---|---|
 | <a id="destexhe2001"></a>**destexhe2001**: Destexhe A, Rudolph M, Fellous J-M, Sejnowski TJ (2001). Fluctuating synaptic conductances recreate in vivo-like activity in neocortical neurons. Neuroscience 107:13–24. | SEC | ADOPTED | Gfluctdv conductance noise (was silent on main); scale calibrated per drive. | [test_gfluctdv_noise](../tests/test_gfluctdv_noise.py), [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py), [CHANGELOG](../CHANGELOG.md) |
-| <a id="heckmanbinder1991"></a>**heckmanbinder1991**: Heckman CJ, Binder MD (1991). Computer simulation of the steady-state input-output function of the cat medial gastrocnemius motoneuron pool. J Neurophysiol 65:952–967. [PMID 2051212](https://pubmed.ncbi.nlm.nih.gov/2051212/) | ABS | ADOPTED | Stated basis of the modelling assumption 'equal effective synaptic current per input' (DD input weights scaled by dendrite→soma transfer). Cat MG data; no data for corticospinal input to a hand muscle. | [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py), [test_dd_optimizer_network](../tests/test_dd_optimizer_network.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [06_optimize_dd_multilevel.py:195](../examples/02_finetune/06_optimize_dd_multilevel.py#L195) |
+| <a id="heckmanbinder1991"></a>**heckmanbinder1991**: Heckman CJ, Binder MD (1991). Computer simulation of the steady-state input-output function of the cat medial gastrocnemius motoneuron pool. J Neurophysiol 65:952–967. [PMID 2051212](https://pubmed.ncbi.nlm.nih.gov/2051212/) | ABS | ADOPTED | Stated basis of the modelling assumption 'equal effective synaptic current per input' (DD input weights scaled by dendrite→soma transfer). Cat MG data; no data for corticospinal input to a hand muscle. | [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py), [test_dd_optimizer_network](../tests/test_dd_optimizer_network.py), [CHANGELOG](../CHANGELOG.md)<br>Cited at: [fdi_semg.yaml:41](../myogen/profiles/fdi_semg.yaml#L41), [descending_drive_network.py:68](../myogen/simulator/neuron/descending_drive_network.py#L68) |
 | <a id="heckmanbinder1993"></a>**heckmanbinder1993**: Heckman CJ, Binder MD (1993). J Neurophysiol 70 (effects of different synaptic input systems on recruitment; details unverified). | LINK | CONTEXT | Likelier source for input-system differences; not read. | — |
 | <a id="beauchamp2023"></a>**beauchamp2023**: Beauchamp JA, Pearcey GEP, Khurram OU, Chardon M, Wang YC, Powers RK, Dewald JPA, Heckman CJ (2023). A geometric approach to quantifying the neuromodulatory effects of persistent inward currents on individual motor unit discharge patterns. J Neural Eng 20(1):016034. [doi:10.1088/1741-2552/acb1d7](https://doi.org/10.1088/1741-2552/acb1d7) · [PMID 36626825](https://pubmed.ncbi.nlm.nih.gov/36626825/) · [PMC9885522](https://pmc.ncbi.nlm.nih.gov/articles/PMC9885522/) | FT | CONTEXT | Human FDI (2296 MUs, triangular ramps to 30 % max): ΔF 4.80 [4.31, 5.30] pps, brace height 34.5 % rTri, acceleration slope 3.12, attenuation slope 0.624 pps/MVT (estimated-metric table). Benchmark for the PIC checks: MyoGen NERLab pool ΔF 0.2 pps. In their model pool ΔF depends on both neuromodulation and the inhibition pattern. | [AlphaMN__Pool](../myogen/simulator/neuron/populations/motor_neurons.py), [alpha_mn_default.yaml](../myogen/config/alpha_mn_default.yaml) |
 | <a id="gorassini2002"></a>**gorassini2002**: Gorassini M, Yang JF, Siu M, Bennett DJ (2002). Intrinsic activation of human motoneurons: possible contribution to motor unit excitation. J Neurophysiol 87:1850–1858. [doi:10.1152/jn.00024.2001](https://doi.org/10.1152/jn.00024.2001) · [PMID 11929906](https://pubmed.ncbi.nlm.nih.gov/11929906/) | LINK | CONTEXT | Origin of the paired motor-unit ΔF method used in the spinal check 2. Not read. | — |
@@ -142,10 +142,10 @@ To add a source: add it to the table in `scripts/build_literature.py`, link it f
 |---|---|---|---|---|
 | <a id="milnerbrown1973a"></a>**milnerbrown1973a**: Milner-Brown HS, Stein RB, Yemm R (1973). The contractile properties of human motor units during voluntary isometric contractions. J Physiol 228:285–306. [doi:10.1113/jphysiol.1973.sp010087](https://doi.org/10.1113/jphysiol.1973.sp010087) | FT | CONTEXT | Methods; twitch shape. | — |
 | <a id="milnerbrown1973b"></a>**milnerbrown1973b**: Milner-Brown HS, Stein RB, Yemm R (1973). The orderly recruitment of human motor units during voluntary isometric contractions. J Physiol 230:359–370. [doi:10.1113/jphysiol.1973.sp010192](https://doi.org/10.1113/jphysiol.1973.sp010192) · [PMC1350367](https://pmc.ncbi.nlm.nih.gov/articles/PMC1350367/) | FT | CONTEXT | 145 FDI units: twitch vs threshold slope ~1, CT 30–100 ms; supports ForceModel defaults; >half of units below ~10 % MVC (threshold check). Twitch range ~100-fold vs model 50-fold (open). | [ForceModel](../myogen/simulator/core/force/force_model.py) |
-| <a id="fuglevand1999"></a>**fuglevand1999**: Fuglevand AJ, Macefield VG, Bigland-Ritchie B (1999). Force-frequency and fatigue properties of motor units in muscles that control digits of the human hand. J Neurophysiol 81(4):1718–1729. [doi:10.1152/jn.1999.81.4.1718](https://doi.org/10.1152/jn.1999.81.4.1718) | FT | ADOPTED | 13 hand units (intraneural stimulation): half-max force 9.1 Hz (CT 65.6 ms) / 15.5 Hz (CT 45.9 ms); 50 Hz gives 91 % of maximum, maximum at 50–100 Hz; twitch:tetanus 0.36 / 0.25; no relation between twitch force and CT. Used in the force-frequency refit (a 0.93, c 1.29). | [ForceModel](../myogen/simulator/core/force/force_model.py), [force_utils](../myogen/simulator/core/force/force_utils.py), [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py)<br>Cited at: [force_utils.py:6](../myogen/simulator/core/force/force_utils.py#L6) |
-| <a id="mcnulty2000"></a>**mcnulty2000**: McNulty PA, Falland KJ, Macefield VG (2000). Comparison of contractile properties of single motor units in human intrinsic and extrinsic finger muscles. J Physiol 526(2):445–456. [doi:10.1111/j.1469-7793.2000.t01-2-00445.x](https://doi.org/10.1111/j.1469-7793.2000.t01-2-00445.x) | FT | ADOPTED | 28 hand units, 14 in FDI: intrinsic CT 70.3 ms (32–111), half-max 11.0 Hz, 34 % at 8 Hz, maximum force at 34.7 Hz (14–66), twitch:tetanus 0.27; CT unrelated to twitch force. Used in the force-frequency refit. | [ForceModel](../myogen/simulator/core/force/force_model.py), [force_utils](../myogen/simulator/core/force/force_utils.py), [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py)<br>Cited at: [force_utils.py:6](../myogen/simulator/core/force/force_utils.py#L6) |
-| <a id="macefield1996"></a>**macefield1996**: Macefield VG, Fuglevand AJ, Bigland-Ritchie B (1996). Contractile properties of single motor units in human toe extensors assessed by intraneural motor axon stimulation. J Neurophysiol 75(6):2509–2519. [doi:10.1152/jn.1996.75.6.2509](https://doi.org/10.1152/jn.1996.75.6.2509) | FT | ADOPTED | 13 toe-extensor units: mean-force curve (Fig. 6C) 48/65/77/87/92 % at 10/15/20/30/50 Hz, CT 74.8 ms; maximum at 15–100 Hz (median 50); twitch:tetanus 0.22 (0.28 potentiated); CT unrelated to twitch force. Fuglevand gain saturates by ~20 Hz at this CT; used in the force-frequency refit (HUMAN_UNIT_FORCE_FREQUENCY, opt-in via ForceModel(force_frequency_shape=...) and --force-frequency human). | [ForceModel](../myogen/simulator/core/force/force_model.py), [force_utils](../myogen/simulator/core/force/force_utils.py), [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py)<br>Cited at: [force_utils.py:6](../myogen/simulator/core/force/force_utils.py#L6) |
-| <a id="fuglevand1993"></a>**fuglevand1993**: Fuglevand AJ, Winter DA, Patla AE (1993). Models of recruitment and rate coding organization in motor-unit pools. J Neurophysiol 70(6):2470–2488. [doi:10.1152/jn.1993.70.6.2470](https://doi.org/10.1152/jn.1993.70.6.2470) | FT | ADOPTED | Twitch/force model: exponential thresholds (p. 2472 eq. 1), twitch force and contraction time (p. 2474 eqs. 13-14), force-frequency saturation S = 1 - exp(-2 (T/ISI)^3) and gain (p. 2475 eqs. 16-17; twitch:tetanus 0.11). Implementation checked against an independent one. | [ForceModel](../myogen/simulator/core/force/force_model.py), [07_fdi_spikes_force](../examples/02_finetune/07_fdi_spikes_force.py), [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py)<br>Cited at: [force_model.py:26](../myogen/simulator/core/force/force_model.py#L26), [force_model.py:28](../myogen/simulator/core/force/force_model.py#L28), [force_model.py:51](../myogen/simulator/core/force/force_model.py#L51), [force_model.py:55](../myogen/simulator/core/force/force_model.py#L55), [force_model.py:62](../myogen/simulator/core/force/force_model.py#L62), [force_model.py:65](../myogen/simulator/core/force/force_model.py#L65), [force_model.py:94](../myogen/simulator/core/force/force_model.py#L94), [force_model.py:153](../myogen/simulator/core/force/force_model.py#L153), [force_model.py:201](../myogen/simulator/core/force/force_model.py#L201), [force_model.py:204](../myogen/simulator/core/force/force_model.py#L204), [force_model.py:266](../myogen/simulator/core/force/force_model.py#L266), [force_model.py:403](../myogen/simulator/core/force/force_model.py#L403), [force_model.py:425](../myogen/simulator/core/force/force_model.py#L425), [force_model_vectorized.py:31](../myogen/simulator/core/force/force_model_vectorized.py#L31), [force_utils.py:4](../myogen/simulator/core/force/force_utils.py#L4), [force_utils.py:7](../myogen/simulator/core/force/force_utils.py#L7), [force_utils.py:7](../myogen/simulator/core/force/force_utils.py#L7), [force_utils.py:13](../myogen/simulator/core/force/force_utils.py#L13), [force_utils.py:21](../myogen/simulator/core/force/force_utils.py#L21), [physiological_distribution.py:35](../myogen/simulator/core/physiological_distribution.py#L35), [physiological_distribution.py:73](../myogen/simulator/core/physiological_distribution.py#L73), [physiological_distribution.py:79](../myogen/simulator/core/physiological_distribution.py#L79), [fr_optimization_targets.py:11](../myogen/utils/fr_optimization_targets.py#L11), [fr_optimization_targets.py:52](../myogen/utils/fr_optimization_targets.py#L52), [fr_optimization_targets.py:327](../myogen/utils/fr_optimization_targets.py#L327), [06_optimize_dd_multilevel.py:6](../examples/02_finetune/06_optimize_dd_multilevel.py#L6), [06_optimize_dd_multilevel.py:96](../examples/02_finetune/06_optimize_dd_multilevel.py#L96), [06_optimize_dd_multilevel.py:310](../examples/02_finetune/06_optimize_dd_multilevel.py#L310), [06_optimize_dd_multilevel.py:1032](../examples/02_finetune/06_optimize_dd_multilevel.py#L1032), [07_fdi_spikes_force.py:8](../examples/02_finetune/07_fdi_spikes_force.py#L8) |
+| <a id="fuglevand1999"></a>**fuglevand1999**: Fuglevand AJ, Macefield VG, Bigland-Ritchie B (1999). Force-frequency and fatigue properties of motor units in muscles that control digits of the human hand. J Neurophysiol 81(4):1718–1729. [doi:10.1152/jn.1999.81.4.1718](https://doi.org/10.1152/jn.1999.81.4.1718) | FT | ADOPTED | 13 hand units (intraneural stimulation): half-max force 9.1 Hz (CT 65.6 ms) / 15.5 Hz (CT 45.9 ms); 50 Hz gives 91 % of maximum, maximum at 50–100 Hz; twitch:tetanus 0.36 / 0.25; no relation between twitch force and CT. Used in the force-frequency refit (a 0.93, c 1.29). | [ForceModel](../myogen/simulator/core/force/force_model.py), [force_utils](../myogen/simulator/core/force/force_utils.py), [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py)<br>Cited at: [fdi_semg.yaml:48](../myogen/profiles/fdi_semg.yaml#L48), [force_utils.py:6](../myogen/simulator/core/force/force_utils.py#L6) |
+| <a id="mcnulty2000"></a>**mcnulty2000**: McNulty PA, Falland KJ, Macefield VG (2000). Comparison of contractile properties of single motor units in human intrinsic and extrinsic finger muscles. J Physiol 526(2):445–456. [doi:10.1111/j.1469-7793.2000.t01-2-00445.x](https://doi.org/10.1111/j.1469-7793.2000.t01-2-00445.x) | FT | ADOPTED | 28 hand units, 14 in FDI: intrinsic CT 70.3 ms (32–111), half-max 11.0 Hz, 34 % at 8 Hz, maximum force at 34.7 Hz (14–66), twitch:tetanus 0.27; CT unrelated to twitch force. Used in the force-frequency refit. | [ForceModel](../myogen/simulator/core/force/force_model.py), [force_utils](../myogen/simulator/core/force/force_utils.py), [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py)<br>Cited at: [fdi_semg.yaml:48](../myogen/profiles/fdi_semg.yaml#L48), [force_utils.py:6](../myogen/simulator/core/force/force_utils.py#L6) |
+| <a id="macefield1996"></a>**macefield1996**: Macefield VG, Fuglevand AJ, Bigland-Ritchie B (1996). Contractile properties of single motor units in human toe extensors assessed by intraneural motor axon stimulation. J Neurophysiol 75(6):2509–2519. [doi:10.1152/jn.1996.75.6.2509](https://doi.org/10.1152/jn.1996.75.6.2509) | FT | ADOPTED | 13 toe-extensor units: mean-force curve (Fig. 6C) 48/65/77/87/92 % at 10/15/20/30/50 Hz, CT 74.8 ms; maximum at 15–100 Hz (median 50); twitch:tetanus 0.22 (0.28 potentiated); CT unrelated to twitch force. Fuglevand gain saturates by ~20 Hz at this CT; used in the force-frequency refit (HUMAN_UNIT_FORCE_FREQUENCY, opt-in via ForceModel(force_frequency_shape=...) and --force-frequency human). | [ForceModel](../myogen/simulator/core/force/force_model.py), [force_utils](../myogen/simulator/core/force/force_utils.py), [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py)<br>Cited at: [fdi_semg.yaml:48](../myogen/profiles/fdi_semg.yaml#L48), [force_utils.py:6](../myogen/simulator/core/force/force_utils.py#L6) |
+| <a id="fuglevand1993"></a>**fuglevand1993**: Fuglevand AJ, Winter DA, Patla AE (1993). Models of recruitment and rate coding organization in motor-unit pools. J Neurophysiol 70(6):2470–2488. [doi:10.1152/jn.1993.70.6.2470](https://doi.org/10.1152/jn.1993.70.6.2470) | FT | ADOPTED | Twitch/force model: exponential thresholds (p. 2472 eq. 1), twitch force and contraction time (p. 2474 eqs. 13-14), force-frequency saturation S = 1 - exp(-2 (T/ISI)^3) and gain (p. 2475 eqs. 16-17; twitch:tetanus 0.11). Implementation checked against an independent one. | [ForceModel](../myogen/simulator/core/force/force_model.py), [07_fdi_spikes_force](../examples/02_finetune/07_fdi_spikes_force.py), [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py)<br>Cited at: [fdi_semg.published.yaml:18](../myogen/profiles/fdi_semg.published.yaml#L18), [fdi_semg.yaml:30](../myogen/profiles/fdi_semg.yaml#L30), [force_model.py:26](../myogen/simulator/core/force/force_model.py#L26), [force_model.py:28](../myogen/simulator/core/force/force_model.py#L28), [force_model.py:51](../myogen/simulator/core/force/force_model.py#L51), [force_model.py:55](../myogen/simulator/core/force/force_model.py#L55), [force_model.py:62](../myogen/simulator/core/force/force_model.py#L62), [force_model.py:65](../myogen/simulator/core/force/force_model.py#L65), [force_model.py:94](../myogen/simulator/core/force/force_model.py#L94), [force_model.py:153](../myogen/simulator/core/force/force_model.py#L153), [force_model.py:201](../myogen/simulator/core/force/force_model.py#L201), [force_model.py:204](../myogen/simulator/core/force/force_model.py#L204), [force_model.py:266](../myogen/simulator/core/force/force_model.py#L266), [force_model.py:403](../myogen/simulator/core/force/force_model.py#L403), [force_model.py:425](../myogen/simulator/core/force/force_model.py#L425), [force_model_vectorized.py:31](../myogen/simulator/core/force/force_model_vectorized.py#L31), [force_utils.py:4](../myogen/simulator/core/force/force_utils.py#L4), [force_utils.py:7](../myogen/simulator/core/force/force_utils.py#L7), [force_utils.py:7](../myogen/simulator/core/force/force_utils.py#L7), [force_utils.py:13](../myogen/simulator/core/force/force_utils.py#L13), [force_utils.py:21](../myogen/simulator/core/force/force_utils.py#L21), [physiological_distribution.py:35](../myogen/simulator/core/physiological_distribution.py#L35), [physiological_distribution.py:73](../myogen/simulator/core/physiological_distribution.py#L73), [physiological_distribution.py:79](../myogen/simulator/core/physiological_distribution.py#L79), [fr_optimization_targets.py:11](../myogen/utils/fr_optimization_targets.py#L11), [fr_optimization_targets.py:52](../myogen/utils/fr_optimization_targets.py#L52), [fr_optimization_targets.py:327](../myogen/utils/fr_optimization_targets.py#L327), [06_optimize_dd_multilevel.py:6](../examples/02_finetune/06_optimize_dd_multilevel.py#L6), [06_optimize_dd_multilevel.py:95](../examples/02_finetune/06_optimize_dd_multilevel.py#L95), [06_optimize_dd_multilevel.py:224](../examples/02_finetune/06_optimize_dd_multilevel.py#L224), [06_optimize_dd_multilevel.py:831](../examples/02_finetune/06_optimize_dd_multilevel.py#L831), [07_fdi_spikes_force.py:8](../examples/02_finetune/07_fdi_spikes_force.py#L8) |
 
 ## G. Experimental EMG benchmarks (FDI amplitude and spectra)
 
diff --git a/examples/01_basic/01_simulate_recruitment_thresholds.py b/examples/01_basic/01_simulate_recruitment_thresholds.py
index cf27138b..ce8e19e9 100644
--- a/examples/01_basic/01_simulate_recruitment_thresholds.py
+++ b/examples/01_basic/01_simulate_recruitment_thresholds.py
@@ -26,7 +26,7 @@
 
 import joblib
 from matplotlib import pyplot as plt
-from myogen import simulator
+from myogen import load_profile, simulator
 
 plt.style.use("fivethirtyeight")
 
@@ -163,6 +163,27 @@
 plt.tight_layout()
 plt.show()
 
+##############################################################################
+# The FDI Profile Card
+# ^^^^^^^^^^^^^^^^^^^^
+#
+# The following examples simulate the **first dorsal interosseous** (FDI) with the ``fdi_semg`` **profile card**,
+# which holds every validated setting with its source (``myogen/profiles/fdi_semg.yaml``). Its thresholds are the
+# ones the FDI drive was calibrated with: 120 units, Fuglevand model, recruitment range 50.
+#
+# !!! note
+#     ``load_profile("fdi_semg.published")`` gives each component as its source paper describes it.
+
+card = load_profile()  # fdi_semg
+recruitment_thresholds = card.thresholds()
+
+plt.plot(recruitment_thresholds * 100, "-o")
+plt.title("FDI Card Recruitment Thresholds")
+plt.xlabel("Motor Unit (#)")
+plt.ylabel("Recruitment Threshold (%)")
+plt.tight_layout()
+plt.show()
+
 ##############################################################################
 # Save Recruitment Thresholds
 # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -170,6 +191,6 @@
 # !!! note
 #     All **MyoGen** objects can be saved to a file using `joblib`. This is useful to **avoid re-running expensive simulations** if you need to use the same parameters.
 
-joblib.dump(combined_results[5], save_path / "thresholds.pkl")
+joblib.dump(recruitment_thresholds, save_path / "thresholds.pkl")
 
 # mkdocs_gallery_thumbnail_path = "gallery_thumbs/01_simulate_recruitment_thresholds.png"
diff --git a/examples/01_basic/03_simulate_spike_trains_descending_drive.py b/examples/01_basic/03_simulate_spike_trains_descending_drive.py
index dfdf398e..7369c3ab 100644
--- a/examples/01_basic/03_simulate_spike_trains_descending_drive.py
+++ b/examples/01_basic/03_simulate_spike_trains_descending_drive.py
@@ -2,23 +2,19 @@
 Spike Train Generation with Descending Drive
 ============================================
 
-This example demonstrates **realistic spike train simulation** using **sinusoidal descending drive (DD)**
-instead of direct current injection. This approach provides more physiologically accurate motor control
-patterns by modeling cortical input through descending drive populations.
+This example simulates **motor neuron spike trains** driven by a **descending drive (DD)**: a population of
+Poisson neurons modelling cortical input, connected to the motor neuron pool through synapses. This is more
+realistic than direct current injection (example 02) because it models the synaptic input from upper motor neurons.
+
+The drive here is the **calibrated FDI drive** of the ``fdi_semg`` profile card: its network, noise and drive rates
+were fitted so the pool recruits and fires like a human first dorsal interosseous, and the card maps force levels in
+**% MVC** to drive. One call runs it.
 
 !!! note
-    This example bridges the gap between simple current injection (example 02) and full spinal network
-    simulation (11_simulate_spinal_network.py). It uses:
-
-    - **DescendingDrive__Pool**: Poisson process neurons modeling cortical input
-    - **AlphaMN__Pool**: Biophysically detailed motor neurons (Powers2017 model)
-    - **Network**: Synaptic connections between DD and motor neuron populations
-    - **Sinusoidal patterns**: Smooth, physiologically relevant input at 0.5-2 Hz
-
-!!! important
-    **Descending Drive (DD)** refers to the cortical and subcortical neural pathways that provide
-    voluntary motor commands to spinal motor neurons. This is more realistic than direct current
-    injection because it models the actual synaptic input patterns from upper motor neurons.
+    The card is ``myogen/profiles/fdi_semg.yaml``: every value with its source, and notes on what else to check when
+    you change one. ``load_profile("fdi_semg.published")`` gives each component as its source paper describes it.
+    To build the pools yourself, see ``card.motor_neuron_pool()`` and
+    ``myogen.simulator.neuron.descending_drive_network``.
 """
 
 # %%
@@ -30,39 +26,18 @@
 # !!! important
 #     In **MyoGen** all **random number generation** is handled by the RNG returned from
 #     [`get_random_generator`][myogen.get_random_generator], a thin wrapper around `numpy.random`.
-#
-#     Always fetch the generator at the call site so the current seed is honored:
-#
-#     ```python
-#     from myogen import simulator, get_random_generator
-#     get_random_generator().normal(0, 1)
-#     ```
-#
-#     To change the seed, use [`set_random_seed`][myogen.set_random_seed]:
-#
-#     ```python
-#     from myogen import set_random_seed
-#     set_random_seed(42)
-#     ```
+#     The profile card reseeds it from its own ``seeds`` before each build, so a card run is reproducible.
 
-# %%
-
-import itertools
 from pathlib import Path
 
 import joblib
 import numpy as np
 import quantities as pq
 from matplotlib import pyplot as plt
-from neo import AnalogSignal, Block, Segment, SpikeTrain
-from neuron import h
-from tqdm import tqdm
+from neo import AnalogSignal
 
-from myogen import get_random_generator
-from myogen.simulator.neuron import Network
-from myogen.simulator.neuron.populations import AlphaMN__Pool, DescendingDrive__Pool
+from myogen import load_profile
 from myogen.utils.nmodl import load_nmodl_mechanisms
-from myogen.utils.types import pps
 
 plt.style.use("fivethirtyeight")
 
@@ -87,269 +62,46 @@ def population_psth(spiketrains, bin_size):
     counts, _ = np.histogram(spikes, bins=edges)
     return counts, edges[:-1]
 
+
 ##############################################################################
-# Create Populations
+# Run the Calibrated Drive
 # ------------------------
 #
-# Like the previous example, we create a **motor neuron pool** using the [`AlphaMN__Pool`][myogen.simulator.neuron.populations.AlphaMN__Pool] class.
-#
-# We also create a [`DescendingDrive__Pool`][myogen.simulator.neuron.populations.DescendingDrive__Pool] to represent the cortical input.
+# A **trapezoidal contraction** to 25 % MVC: 1 s rest, 4 s ramp up, 4 s hold, 4 s ramp down, 1 s rest. Each
+# segment is ``(duration_s, start_pct_mvc, end_pct_mvc)``.
 #
 # !!! note
-#     These neurons are modeled as Poisson point processes to convert the smooth input signal into realistic
-#     spike patterns that represent cortical input to the spinal cord.
-#
+#     This runs 120 biophysical motor neurons for 14 s and takes a few minutes.
 
 load_nmodl_mechanisms()
 
 save_path = Path("./results")
 save_path.mkdir(exist_ok=True)
 
-recruitment_thresholds = joblib.load(save_path / "thresholds.pkl")
+card = load_profile()  # fdi_semg
+protocol = [(1, 0, 0), (4, 0, 25), (4, 25, 25), (4, 25, 0), (1, 0, 0)]
+result = card.simulate(protocol)
 
-motor_neuron_pool = AlphaMN__Pool(
-    recruitment_thresholds__array=recruitment_thresholds,
-    config_file="alpha_mn_default.yaml",
-)
-
-timestep = 0.1 * pq.ms
-h.secondorder = 2  # Crank-Nicolson method (second-order accurate)
-descending_drive_pool = DescendingDrive__Pool(n=100, process_type="poisson", timestep__ms=timestep)
-##############################################################################
-# Generate Trapezoidal Drive Pattern
-# -----------------------------------
-#
-# Create a **trapezoidal ramp contraction pattern** that represents realistic
-# voluntary isometric contractions. This pattern has 4 phases:
-# 1. **Ramp-up**: Linear increase from baseline to peak
-# 2. **Plateau**: Sustained peak drive level
-# 3. **Ramp-down**: Linear decrease from peak to baseline
-# 4. **Rest**: Baseline activity
-#
-# This is a common experimental paradigm used in motor control studies.
-
-simulation_time = 15000 * pq.ms
-time_points = int(simulation_time / timestep)
-
-# Trapezoidal parameters
-dd_baseline__pps = 0.0 * pps  # Baseline drive during rest
-dd_peak__pps = 65 * pps  # Peak drive during plateau
-
-# Phase durations (ms) - Total trapezoid duration: 11000ms
-ramp_up_duration = 500 * pq.ms  # 500ms ramp up
-plateau_duration = 10000 * pq.ms  # 10s hold
-ramp_down_duration = 500 * pq.ms  # 500ms ramp down
-
-# Add rest periods before and after
-rest_before = 1000 * pq.ms  # 1s rest before trapezoid
-rest_after = 1000 * pq.ms  # 1s rest after trapezoid
-
-# Center the trapezoid at 7.5s (middle of 15s simulation)
-# Calculate phase boundaries with rest period before
-trapezoid_start = rest_before  # Start at 1s
-ramp_up_end = trapezoid_start + ramp_up_duration  # 1.5s
-plateau_end = ramp_up_end + plateau_duration  # 11.5s
-ramp_down_end = plateau_end + ramp_down_duration  # 12s
-
-# Create time array
-time_array = np.linspace(0, simulation_time.magnitude, time_points) * pq.ms
-
-# Initialize drive signal (all baseline)
-trapezoid_drive = np.ones(time_points) * dd_baseline__pps
-
-for i, t in enumerate(time_array):
-    if t < trapezoid_start:
-        # Phase 0: Rest before
-        trapezoid_drive[i] = dd_baseline__pps
-    elif t < ramp_up_end:
-        # Phase 1: Ramp up
-        elapsed = t - trapezoid_start
-        trapezoid_drive[i] = dd_baseline__pps + (dd_peak__pps - dd_baseline__pps) * (
-            elapsed / ramp_up_duration
-        )
-    elif t < plateau_end:
-        # Phase 2: Plateau
-        trapezoid_drive[i] = dd_peak__pps
-    elif t < ramp_down_end:
-        # Phase 3: Ramp down
-        elapsed = t - plateau_end
-        trapezoid_drive[i] = dd_peak__pps - (dd_peak__pps - dd_baseline__pps) * (
-            elapsed / ramp_down_duration
-        )
-    else:
-        # Phase 4: Rest after
-        trapezoid_drive[i] = dd_baseline__pps
-
-# Add small noise for realism
-trapezoid_drive = (
-    trapezoid_drive + np.clip(get_random_generator().normal(0, 1.0, size=time_points), 0, None) * pps
-)
+spike_train_block = result.spikes
+mn_segment = spike_train_block.segments[0]
+n_motor_neurons = len(mn_segment.spiketrains)
+simulation_time = sum(d for d, _, _ in protocol) * pq.s
+timestep = float(card.value("drive.timestep_ms")) * pq.ms
+ramp_up_end = 5.0 * pq.s  # start of the hold
+plateau_end = 9.0 * pq.s  # end of the hold
 
-# Create AnalogSignal
-trapezoid_drive_signal = AnalogSignal(
-    signal=trapezoid_drive, sampling_period=timestep.rescale(pq.s)
-)
-
-joblib.dump(trapezoid_drive_signal, save_path / "trapezoid_drive_pattern.pkl")
-print(
-    f"\n Trapezoidal drive pattern (1000ms trapezoid centered in {simulation_time}ms simulation):"
-)
-print(f"\tRest before: 0 - {trapezoid_start} ms ({dd_baseline__pps} pps)")
-print(f"\tRamp up: {trapezoid_start} - {ramp_up_end} ms ({dd_baseline__pps} → {dd_peak__pps} pps)")
-print(
-    f"\tPlateau: {ramp_up_end} - {plateau_end} ms ({dd_peak__pps} pps, center at {(ramp_up_end + plateau_end) / 2:.0f}ms)"
-)
-print(f"\tRamp down: {plateau_end} - {ramp_down_end} ms ({dd_peak__pps} → {dd_baseline__pps} pps)")
-print(f"\tRest after: {ramp_down_end} - {simulation_time} ms ({dd_baseline__pps} pps)")
-
-##############################################################################
-# Create Network and Connections
-# -------------------------------
-#
-# In MyoGen, populations can be connected using the [`Network`][myogen.simulator.Network] class from the
-# `myogen.simulator.neuron` module.
-#
-# The **Network** class provides a high-level interface for creating and managing
-# connections between neuron populations.
-
-# Use the **Network** class to create synaptic connections between the descending drive
-# population and the motor neuron pool. This creates realistic synaptic transmission
-# with appropriate delays and weights.
-#
-
-network = Network({"DD": descending_drive_pool, "aMN": motor_neuron_pool})
-
-# Connect DD neurons to motor neurons with realistic synaptic parameters
-network.connect(source="DD", target="aMN", probability=0.5, weight__uS=0.15 * pq.uS)
-
-# Set up external input to DD population
-network.connect_from_external(source="cortical_input", target="DD", weight__uS=1.0 * pq.uS)
-
-# Get NetCons for manual DD stimulation
-dd_netcons = network.get_netcons("cortical_input", "DD")
-
-##############################################################################
-# Setup Spike Recording
-# ---------------------
-#
-# To record spikes, we need to manually set up spike detection for the motor neurons
-# and track spike times for the DD neurons.
-#
-
-# Manual spike tracking for DD neurons (they use Poisson processes)
-dd_spike_times = [[] for _ in range(len(descending_drive_pool))]
-
-# Record spikes from motor neurons
-mn_spike_recorders = []
-for cell in motor_neuron_pool:
-    spike_recorder = h.Vector()
-    nc = h.NetCon(cell.soma(0.5)._ref_v, None, sec=cell.soma)
-    nc.threshold = 50  # Standard threshold for motor neurons
-    nc.record(spike_recorder)
-    mn_spike_recorders.append(spike_recorder)
-
-##############################################################################
-# Run Simulation
-# ------------------------------------
-#
-# Execute the NEURON simulation with real-time injection of the sinusoidal drive pattern.
-# The DD neurons receive time-varying input that drives their Poisson processes.
-
-h.load_file("stdrun.hoc")  # Load standard run library for NEURON
-h.dt = timestep
-h.tstop = simulation_time
-
-# Initialize voltages for all pools
-for section, voltage in itertools.chain.from_iterable(
-    zip(*pool.get_initialization_data()) for pool in [motor_neuron_pool, descending_drive_pool]
-):
-    section.v = voltage
-
-
-h.finitialize()
-
-# Calculate total simulation steps for progress bar
-total_steps = int(simulation_time / timestep)
-
-step_counter = 0
-with tqdm(
-    total=float(simulation_time),
-    desc="Running simulation",
-    unit="ms",
-    bar_format="{l_bar}{bar}| {n:.2f}/{total:.2f} ms [{elapsed}<{remaining}, {rate_fmt}]",
-) as pbar:
-    while h.t < h.tstop:
-        current_drive = trapezoid_drive_signal[min(step_counter, len(trapezoid_drive_signal) - 1)]
-
-        # Drive DD neurons with current input level
-        for dd_cell in descending_drive_pool:
-            if dd_cell.integrate(current_drive):
-                # Record spike time for DD neuron
-                dd_spike_times[dd_cell.pool__ID].append(h.t)
-                # Generate spike in DD neuron
-                spike_time = h.t + 1
-                if spike_time < h.tstop:  # Avoid scheduling beyond simulation end
-                    dd_netcons[dd_cell.pool__ID].event(spike_time)
-
-        # Progress simulation
-        h.fadvance()
-        step_counter += 1
-        pbar.update(float(timestep))
-
-
-##############################################################################
-# Convert Spike Data to Neo Format
-# ---------------------------------
-#
-
-spike_train_block = Block(name="Trapezoidal DD Spike Trains")
-
-dd_segment = Segment(name="Descending Drive")
-dd_segment.spiketrains = [
-    SpikeTrain(
-        (spike_times * pq.ms).rescale(pq.s),  # type: ignore
-        t_stop=simulation_time.rescale(pq.s),
-        sampling_rate=(1 / (h.dt * pq.ms)).rescale(pq.Hz),
-        sampling_period=h.dt * pq.ms,
-        name=f"DD_{i}",
-    )
-    for i, spike_times in enumerate(dd_spike_times)
-]
-
-mn_segment = Segment(name="Motor Neurons")
-mn_segment.spiketrains = [
-    SpikeTrain(
-        (recorder.as_numpy() * pq.ms).rescale(pq.s),  # type: ignore
-        t_stop=simulation_time.rescale(pq.s),
-        sampling_rate=(1 / (h.dt * pq.ms)).rescale(pq.Hz),
-        sampling_period=h.dt * pq.ms,
-        name=f"MN_{i}",
-    )
-    for i, recorder in enumerate(mn_spike_recorders)
-]
-
-# We only save the motor neuron spikes  segment
-spike_train_block.segments.append(mn_segment)
+# Excitation (0-1) per drive time step, saved for the surface EMG example
+excitation_signal = AnalogSignal(result.excitation[:, None], units=pq.dimensionless,
+                                 sampling_period=timestep.rescale(pq.s))
 
+# The spike trains and the excitation are used by the surface EMG example (06) and the force example (08)
 joblib.dump(spike_train_block, save_path / "trapezoid_dd_spike_trains.pkl")
+joblib.dump(excitation_signal, save_path / "trapezoid_drive_pattern.pkl")
 
 ##############################################################################
 # Calculate Firing Rate Statistics
 # ---------------------------------
-#
-
-print("\nFiring rate analysis:")
-
-# Calculate DD firing rates
-dd_firing_rates = np.array(
-    [
-        mean_firing_rate(st__s.time_slice(st__s.min(), st__s.max()))
-        for st__s in dd_segment.spiketrains
-        if len(st__s) > 1  # need >=2 spikes for a rate over the spike span
-    ]
-)
 
-# Calculate MN firing rates
 mn_firing_rates = np.array(
     [
         mean_firing_rate(st__s.time_slice(st__s.min(), st__s.max()))
@@ -358,97 +110,66 @@ def population_psth(spiketrains, bin_size):
     ]
 )
 
-print("Descending Drive neurons:")
-print(f"\tActive neurons: {len(dd_firing_rates)}/{descending_drive_pool.n}")
-if len(dd_firing_rates) > 0:
-    print(f"\tMean firing rate: {np.mean(dd_firing_rates):.1f} ± {np.std(dd_firing_rates):.1f} pps")
-    print(f"\tRate range: {np.min(dd_firing_rates):.1f} - {np.max(dd_firing_rates):.1f} pps")
-
 print("Motor neurons:")
-print(f"\tActive neurons: {len(mn_firing_rates)}/{motor_neuron_pool.n}")
+print(f"\tActive neurons: {len(mn_firing_rates)}/{n_motor_neurons}")
 if len(mn_firing_rates) > 0:
     print(f"\tMean firing rate: {np.mean(mn_firing_rates):.1f} ± {np.std(mn_firing_rates):.1f} pps")
     print(f"\tRate range: {np.min(mn_firing_rates):.1f} - {np.max(mn_firing_rates):.1f} pps")
 
 ##############################################################################
-# Advanced Visualization
-# -----------------------
+# Visualization
+# -------------
 #
-# Create comprehensive visualizations showing:
-# 1. Sinusoidal drive input pattern
-# 2. DD population raster plot with drive overlay
-# 3. Motor neuron raster plot showing recruitment
-# 4. Population firing rates over time
-
-# Create figure with subplots
-fig, axes = plt.subplots(4, 1, figsize=(15, 12), sharex=True)
-
-# 1. Plot trapezoidal drive pattern
-time_s = trapezoid_drive_signal.times.rescale(pq.s).magnitude
-axes[0].plot(time_s, trapezoid_drive_signal, "b-", linewidth=2, label="DD Input")
-axes[0].axhline(dd_baseline__pps, color="r", linestyle="--", alpha=0.7, label="Baseline")
-axes[0].set_ylabel("Drive (Hz)")
-axes[0].set_title("Trapezoidal Descending Drive Pattern (Ramp Contraction)")
+# 1. Target force and the simulated force (% MVC)
+# 2. Motor neuron raster plot showing recruitment
+# 3. Population firing rate over time
+
+fig, axes = plt.subplots(3, 1, figsize=(15, 10), sharex=True)
+
+t_force = np.arange(result.force_pct_mvc.size) / result.fs_force__Hz
+t_target = np.cumsum([0] + [d for d, _, _ in protocol])
+target = [protocol[0][1]] + [b for _, _, b in protocol]
+axes[0].plot(t_target, target, "k--", linewidth=2, label="Target")
+axes[0].plot(t_force, result.force_pct_mvc, "b-", linewidth=2, label="Force")
+axes[0].set_ylabel("Force (% MVC)")
+axes[0].set_title("Trapezoidal Contraction (calibrated FDI drive)")
 axes[0].legend()
 axes[0].grid(True, alpha=0.3)
 
-# 2. DD population raster plot
-dd_colors = plt.get_cmap("Blues")(np.linspace(0.3, 0.8, len(dd_segment.spiketrains)))
-for i, (spiketrain, color) in enumerate(zip(dd_segment.spiketrains, dd_colors)):
-    if len(spiketrain) > 0:
-        axes[1].scatter(spiketrain.magnitude, [i] * len(spiketrain), c=[color], s=0.8, alpha=0.8)
-
-axes[1].set_ylabel("DD Neuron ID")
-axes[1].set_title(f"Descending Drive Population Activity (n={descending_drive_pool.n})")
-axes[1].set_ylim(-1, descending_drive_pool.n)
-axes[1].grid(True, alpha=0.3)
-
-# 3. Motor neuron raster plot (recruitment ordered)
-mn_colors = plt.get_cmap("Reds")(np.linspace(0.3, 0.9, len(mn_segment.spiketrains)))
+mn_colors = plt.get_cmap("Reds")(np.linspace(0.3, 0.9, n_motor_neurons))
 active_mn_count = 0
 for i, (spiketrain, color) in enumerate(zip(mn_segment.spiketrains, mn_colors)):
     if len(spiketrain) > 0:
         spike_times = spiketrain.rescale(pq.s).magnitude
-        axes[2].scatter(spike_times, [i] * len(spike_times), c=[color], s=1.0, alpha=0.8)
+        axes[1].scatter(spike_times, [i] * len(spike_times), c=[color], s=1.0, alpha=0.8)
         active_mn_count += 1
 
-axes[2].set_ylabel("Motor Neuron ID\n(Recruitment Order)")
-axes[2].set_title(
-    f"Motor Neuron Population Activity (n={active_mn_count}/{motor_neuron_pool.n} active)"
-)
-axes[2].set_ylim(-1, motor_neuron_pool.n)
-axes[2].grid(True, alpha=0.3)
+axes[1].set_ylabel("Motor Neuron ID\n(Recruitment Order)")
+axes[1].set_title(f"Motor Neuron Population Activity (n={active_mn_count}/{n_motor_neurons} active)")
+axes[1].set_ylim(-1, n_motor_neurons)
+axes[1].grid(True, alpha=0.3)
 
-# 4. Population firing rates over time (binned)
 bin_size_ms = 100
+mn_counts, bin_centers_s = population_psth(mn_segment.spiketrains, bin_size_ms * pq.ms)
+mn_rates_binned = mn_counts / (bin_size_ms / 1000.0) / n_motor_neurons
+axes[2].plot(bin_centers_s, mn_rates_binned, "r-", linewidth=2, label="MN Population", alpha=0.8)
+axes[2].set_xlabel("Time (s)")
+axes[2].set_ylabel("Population Rate (Hz)")
+axes[2].set_title("Population Firing Rate Over Time")
+axes[2].legend()
+axes[2].grid(True, alpha=0.3)
 
-dd_counts, bin_centers_s = population_psth(dd_segment.spiketrains, bin_size_ms * pq.ms)
-dd_rates_binned = dd_counts / (bin_size_ms / 1000.0) / descending_drive_pool.n
-
-mn_counts, _ = population_psth(mn_segment.spiketrains, bin_size_ms * pq.ms)
-mn_rates_binned = mn_counts / (bin_size_ms / 1000.0) / motor_neuron_pool.n
-axes[3].plot(bin_centers_s, dd_rates_binned, "b-", linewidth=2, label="DD Population", alpha=0.8)
-axes[3].plot(bin_centers_s, mn_rates_binned, "r-", linewidth=2, label="MN Population", alpha=0.8)
-
-axes[3].set_xlabel("Time (s)")
-axes[3].set_ylabel("Population Rate (Hz)")
-axes[3].set_title("Population Firing Rates Over Time")
-axes[3].legend()
-axes[3].grid(True, alpha=0.3)
-
-# Format all axes
 for ax in axes:
     ax.set_xlim(0, simulation_time.rescale(pq.s).magnitude)
 
 plt.tight_layout()
 plt.show()
-
 ##############################################################################
 # Individual Motor Neuron Discharge Rates
 # ----------------------------------------
 #
 # Compute smoothed instantaneous firing rates for each motor neuron
-# using Elephant's kernel-based rate estimation.
+# with a 400 ms Hanning window.
 
 print("\nComputing smoothed discharge rates per neuron...")
 
@@ -485,7 +206,7 @@ def population_psth(spiketrains, bin_size):
         plateau_spiketrain = spiketrain.time_slice(ramp_up_end, plateau_end)
 
         # Compute mean firing rate (Hz) from plateau spikes only
-        plateau_duration_s = (plateau_end - ramp_up_end) / 1000.0
+        plateau_duration_s = float((plateau_end - ramp_up_end).rescale(pq.s).magnitude)
         mean_rate = (
             len(plateau_spiketrain) / plateau_duration_s if len(plateau_spiketrain) > 0 else 0.0
         )
@@ -567,6 +288,6 @@ def population_psth(spiketrains, bin_size):
 plt.tight_layout()
 plt.show()
 
-print("\n[DONE] Simulation complete with individual neuron noise and discharge rate analysis!")
+print("\n[DONE] Simulation complete with discharge rate analysis!")
 
 # mkdocs_gallery_thumbnail_path = "gallery_thumbs/03_simulate_spike_trains_descending_drive.png"
diff --git a/examples/01_basic/04_simulate_muscle.py b/examples/01_basic/04_simulate_muscle.py
index 920708f0..7a2c8468 100644
--- a/examples/01_basic/04_simulate_muscle.py
+++ b/examples/01_basic/04_simulate_muscle.py
@@ -18,9 +18,8 @@
 import joblib
 import matplotlib.pyplot as plt
 import numpy as np
-import quantities as pq
 
-from myogen import simulator
+from myogen import load_profile
 from myogen.utils.plotting.muscle import plot_innervation_areas_2d, plot_mf_centers
 
 plt.style.use("fivethirtyeight")
@@ -29,32 +28,20 @@
 # Define Parameters
 # -----------------
 #
-# The **muscle model** is created using the [`Muscle`][myogen.simulator.Muscle] object.
-#
-# The **Muscle** object takes the following parameters:
+# The **muscle model** is a [`Muscle`][myogen.simulator.Muscle] object. Its main parameters are:
 #
 # - ``recruitment_thresholds``: Recruitment thresholds of the motor units
-# - ``radius``: Radius of the muscle in mm
-# - ``fiber_density``: Fiber density per mm²
-# - ``max_innervation_area_to_total_muscle_area__ratio``: Maximum innervation area to total muscle area ratio
-# - ``grid_resolution``: Spatial resolution for muscle discretization
+# - ``radius__mm``: Radius of the muscle
+# - ``fiber_density__fibers_per_mm2``: Fiber density
+# - ``skin_thickness__mm``, ``fat_thickness__mm``: Layers between the muscle and the electrodes
 #
-# Since the **recruitment thresholds** are already generated, we can load them from the previous example using ``joblib``.
+# The ``fdi_semg`` **profile card** holds the FDI values with their sources (radius 6.91 mm, 267 fibres/mm²,
+# skin 1.2 mm, fat 1.3 mm, per-unit conduction velocities from recruitment force, ...). ``card.muscle()`` builds
+# the muscle with them and the card's recruitment thresholds (regenerated from its seed: those of example 01). Keywords override a card value, e.g.
+# ``card.muscle(fat_thickness__mm=3 * pq.mm)``; the card then prints a note on what else that affects.
 
-# Load recruitment thresholds
 save_path = Path("./results")
-
-recruitment_thresholds = joblib.load(save_path / "thresholds.pkl")
-
-# Define muscle parameters
-muscle_radius = 4.9 * pq.mm  # Muscle radius in mm
-mean_fiber_length = 32 * pq.mm  # Mean fiber length in mm
-fiber_length_variation = 3 * pq.mm  # Fiber length variation in mm
-fiber_density = 400 * pq.mm**-2  # Fiber density in fibers per mm²
-
-# Define simulation parameters
-max_innervation_ratio = 1 / 4  # Maximum motor unit territory size
-grid_resolution = 256  # Spatial resolution for muscle discretization
+card = load_profile()  # fdi_semg
 
 ##############################################################################
 # Create Muscle Model
@@ -65,14 +52,8 @@
 #
 #     To **avoid running the simulation every time**, we can save the muscle model using `joblib`.
 
-# Create muscle model
-muscle = simulator.Muscle(
-    recruitment_thresholds=recruitment_thresholds,
-    radius_bone__mm=1.0 * pq.mm,  # Non-zero bone radius required for proper MUAP amplitudes
-    fiber_density__fibers_per_mm2=fiber_density,
-    fat_thickness__mm=10 * pq.mm,
-    autorun=True,
-)
+# Create muscle model (innervation areas, fibres, fibre-to-unit assignment, fibre properties)
+muscle = card.muscle()
 
 # Save muscle model for future use
 joblib.dump(muscle, save_path / "muscle_model.pkl")
@@ -82,8 +63,8 @@
 
 print("Muscle model statistics:")
 print(f"\tTotal muscle fibers: {total_fibers}")
-print(f"\tMean fibers per MU: {total_fibers / len(recruitment_thresholds):.1f}")
-print(f"\tMuscle cross-sectional area: {np.pi * muscle_radius**2:.1f}")
+print(f"\tMean fibers per MU: {total_fibers / len(muscle.resulting_number_of_innervated_fibers):.1f}")
+print(f"\tMuscle cross-sectional area: {np.pi * muscle.radius__mm**2:.1f}")
 
 ##############################################################################
 # Visualize Muscle Fiber Centers
diff --git a/examples/01_basic/05_simulate_surface_muaps.py b/examples/01_basic/05_simulate_surface_muaps.py
index d8702291..7a4ac229 100644
--- a/examples/01_basic/05_simulate_surface_muaps.py
+++ b/examples/01_basic/05_simulate_surface_muaps.py
@@ -17,7 +17,7 @@
 import numpy as np
 import quantities as pq
 
-from myogen import simulator
+from myogen import load_profile, simulator
 from myogen.utils.neo import signal_to_grid
 
 plt.style.use("fivethirtyeight")
@@ -71,19 +71,17 @@
     bending_radius__mm=muscle.radius__mm + muscle.skin_thickness__mm + muscle.fat_thickness__mm,
 )
 
-surface_emg = simulator.SurfaceEMG(
-    muscle_model=muscle,
-    electrode_arrays=[electrode_array_monopolar],
+# The ``fdi_semg`` card sets the volume conductor and action-potential settings it was validated with (fibre
+# dispersion, Rosenfalck lambda 0.63 fitted to Botelho 2019, 200 mm IAP kernel); the electrodes, the sampling rate
+# and the units to simulate are yours.
+card = load_profile()  # fdi_semg
+surface_emg = card.surface_emg(
+    muscle,
+    [electrode_array_monopolar],
     sampling_frequency__Hz=sampling_frequency,
     sampling_points_in_t_and_z_domains=256,
-    sampling_points_in_theta_domain=32,  # Restored to default (log-space Bessel implementation prevents overflow)
+    sampling_points_in_theta_domain=32,
     MUs_to_simulate=[0, 1, 2, 3],  # Simulate MUs 0, 1, 2, and 3
-    # iap_kernel_length__mm=None,  # Default: uses scaled fiber length (2.5× for boundary clearance)
-    # The IAP kernel is automatically scaled to 2.5× the mean fiber length to:
-    # 1. Prevent boundary truncation artifacts
-    # 2. Allow full IAP wavefront development and decay
-    # 3. Maintain proportionality to actual fiber lengths
-    # Result: Realistic MUAP durations (8-10 ms) independent of sampling resolution N!
 )
 
 ##############################################################################
diff --git a/examples/01_basic/06_simulate_surface_emg.py b/examples/01_basic/06_simulate_surface_emg.py
index 9507a578..d711f6bb 100644
--- a/examples/01_basic/06_simulate_surface_emg.py
+++ b/examples/01_basic/06_simulate_surface_emg.py
@@ -25,7 +25,7 @@
 
 from myogen import simulator
 from myogen.utils.neo import signal_to_grid
-from myogen.utils.types import CURRENT__AnalogSignal, SPIKE_TRAIN__Block
+from myogen.utils.types import SPIKE_TRAIN__Block
 
 plt.style.use("fivethirtyeight")
 
@@ -37,9 +37,7 @@
 save_path = Path("./results")
 
 spike_train__Block: SPIKE_TRAIN__Block = joblib.load(save_path / "trapezoid_dd_spike_trains.pkl")
-input_current__AnalogSignal: CURRENT__AnalogSignal = joblib.load(
-    save_path / "trapezoid_drive_pattern.pkl"
-)
+excitation__AnalogSignal = joblib.load(save_path / "trapezoid_drive_pattern.pkl")  # drive excitation (0-1)
 surface_emg: simulator.SurfaceEMG = joblib.load(save_path / "surface_emg.pkl")
 
 ##############################################################################
@@ -75,8 +73,6 @@
 
 noisy_surface_emg__Block = surface_emg.add_noise(snr__dB=5.0)
 
-# Load input current as AnalogSignal
-
 plt.rcParams.update({"font.size": 24})
 # Create single plot with normalized signals
 fig, ax = plt.subplots(figsize=(12, 6))
@@ -85,14 +81,12 @@
 noisy_signal = noisy_surface_emg__Block.groups[0].segments[0].analogsignals[0]
 emg_grid = signal_to_grid(noisy_signal)  # Convert to 3D (time, rows, cols)
 emg_signal = emg_grid[:, 2, 2]
-current_signal = input_current__AnalogSignal[:, 0].magnitude
+excitation = excitation__AnalogSignal[:, 0].magnitude
 
-# Normalize current between 0 and 1
-current_normalized = (current_signal - np.min(current_signal)) / (
-    np.max(current_signal) - np.min(current_signal)
-)
+# Normalize the excitation between 0 and 1
+excitation_normalized = (excitation - np.min(excitation)) / (np.max(excitation) - np.min(excitation))
 
-current_normalized *= np.max(emg_signal)  # Scale to match EMG amplitude range
+excitation_normalized *= np.max(emg_signal)  # Scale to match EMG amplitude range
 
 # Plot both normalized signals on same axis
 ax.plot(
@@ -103,10 +97,10 @@
 )
 
 ax.plot(
-    input_current__AnalogSignal.times.rescale(pq.s).magnitude,
-    current_normalized,
+    excitation__AnalogSignal.times.rescale(pq.s).magnitude,
+    excitation_normalized,
     linewidth=2,
-    label="Input Current",
+    label="Excitation",
     alpha=1,
     zorder=4,
 )
@@ -116,7 +110,7 @@
 
 sns.despine(trim=True, left=False, bottom=False, right=True, top=True, offset=5)
 
-plt.title("Surface EMG\nand Input Current")
+plt.title("Surface EMG\nand Drive Excitation")
 
 plt.tight_layout()
 plt.show()
diff --git a/examples/01_basic/08_simulate_force.py b/examples/01_basic/08_simulate_force.py
index e5ea81a6..62123261 100644
--- a/examples/01_basic/08_simulate_force.py
+++ b/examples/01_basic/08_simulate_force.py
@@ -29,86 +29,34 @@
 # ----------------
 
 # %%
+from pathlib import Path
+
+import joblib
 import matplotlib.pyplot as plt
 import numpy as np
-import quantities as pq
-
-from myogen import simulator
-from myogen.simulator.core.force.force_model import ForceModel
-from myogen.simulator.neuron.populations import AlphaMN__Pool
-from myogen.utils.currents import create_trapezoid_current
-from myogen.utils.neuron.inject_currents_into_populations import (
-    inject_currents_and_simulate_spike_trains,
-)
-from myogen.utils.nmodl import load_nmodl_mechanisms
+
+from myogen import load_profile
 from myogen.utils.plotting.force import plot_twitch_parameter_assignment, plot_twitches
 
 plt.style.use("fivethirtyeight")
 
 ##############################################################################
-# Load NMODL Mechanisms
-# ---------------------
-#
-# Load the custom NEURON mechanisms required for the motor neuron models.
-
-load_nmodl_mechanisms()
-
-##############################################################################
-# Define Parameters
-# -----------------
+# Create Force Model
+# ------------------
 #
-# The **force model** requires several key parameters:
+# The [`ForceModel`][myogen.simulator.ForceModel] calculates each motor unit's twitch (peak force and contraction
+# time) from its recruitment threshold. Its main parameters are:
 #
 # - ``recruitment_thresholds``: Motor unit recruitment thresholds
-# - ``recording_frequency__Hz``: Sampling frequency for force simulation
-# - ``longest_duration_rise_time__ms``: Maximum twitch rise time
-# - ``contraction_time_range``: Range of contraction times across motor units
-#
-# For demonstration, we'll also create input currents and simulate the complete
-# pipeline from current input to force output.
-
-# Motor unit pool parameters
-n_motor_units = 50
-recruitment_range = 50  # Recruitment range (max_threshold / min_threshold)
-
-# Force model parameters
-recording_frequency__Hz = 2048 * pq.Hz  # 2048 Hz sampling rate
-longest_duration_rise_time__ms = 90.0 * pq.ms  # Maximum twitch rise time
-contraction_time_range = 3  # Contraction time range factor
-
-# Simulation parameters
-simulation_duration__ms = 10000.0  # 10 seconds
-timestep__ms = 0.05 * pq.ms  # 0.05 ms time step
-t_points = int(simulation_duration__ms / timestep__ms.magnitude)
-
-##############################################################################
-# Generate Recruitment Thresholds
-# --------------------------------
+# - ``recording_frequency__Hz``: Sampling frequency of the force
+# - ``force_frequency_shape``: Shape of the force-frequency (gain) curve
 #
-# First, we need to generate recruitment thresholds for our motor unit pool.
-# These will determine both the activation order and the force characteristics
-# of each motor unit.
-
-recruitment_thresholds, _ = simulator.RecruitmentThresholds(
-    N=n_motor_units,
-    recruitment_range__ratio=recruitment_range,
-    mode="combined",
-    deluca__slope=5,
-)
+# The ``fdi_semg`` profile card builds it with the thresholds and the human force-frequency curve
+# (Macefield 1996, Fuglevand 1999) the FDI drive was calibrated with.
+# ``load_profile("fdi_semg.published").force_model()`` uses the curve of Fuglevand 1993 instead.
 
-##############################################################################
-# Create Force Model
-# ------------------
-#
-# The force model calculates individual motor unit twitch properties based on
-# recruitment thresholds and physiological scaling rules.
-
-force_model = ForceModel(
-    recruitment_thresholds=recruitment_thresholds,
-    recording_frequency__Hz=recording_frequency__Hz,
-    longest_duration_rise_time__ms=longest_duration_rise_time__ms,
-    contraction_time_range_factor=contraction_time_range,
-)
+card = load_profile()  # fdi_semg
+force_model = card.force_model()
 
 # Display force model statistics
 print("Force model statistics:")
@@ -144,48 +92,13 @@
 plt.show()
 
 ##############################################################################
-# Generate Input Current
-# ----------------------
-#
-# To demonstrate force generation, we'll create a trapezoid current that will
-# drive motor unit recruitment and firing patterns.
-
-# Parameters for trapezoid current
-trap_amplitude = 15.0 * pq.nA  # Peak amplitude
-trap_rise_time = 5000.0 * pq.ms  # Rise duration (ms)
-trap_plateau_time = 8000.0 * pq.ms  # Plateau duration (ms)
-trap_fall_time = 3000.0 * pq.ms  # Fall duration (ms)
-trap_offset = 5.0 * pq.nA  # Baseline current
-trap_delay = 0.0 * pq.ms  # Initial delay (ms)
-
-input_current__AnalogSignal = create_trapezoid_current(
-    n_pools=1,
-    t_points=t_points,
-    timestep__ms=timestep__ms,
-    amplitudes__nA=[trap_amplitude],
-    rise_times__ms=[trap_rise_time],
-    plateau_times__ms=[trap_plateau_time],
-    fall_times__ms=[trap_fall_time],
-    offsets__nA=[trap_offset],
-    delays__ms=[trap_delay],
-)
-
-##############################################################################
-# Create Motor Neuron Pool and Generate Spike Trains
-# ---------------------------------------------------
+# Load Spike Trains
+# -----------------
 #
-# Now we'll create a motor neuron pool and generate spike trains in response
-# to the input current.
-
+# The spike trains of the 25 % MVC contraction from example 03 (calibrated FDI drive).
 
-motor_neuron_pool = AlphaMN__Pool(recruitment_thresholds__array=recruitment_thresholds)
-
-# Generate spike trains using the convenient utility function
-spike_train__Block = inject_currents_and_simulate_spike_trains(
-    populations=[motor_neuron_pool],
-    input_current__AnalogSignal=input_current__AnalogSignal,
-    spike_detection_thresholds__mV=50 * pq.mV,
-)
+save_path = Path("./results")
+spike_train__Block = joblib.load(save_path / "trapezoid_dd_spike_trains.pkl")
 
 ##############################################################################
 # Generate Force Output
@@ -195,8 +108,9 @@
 # The model will simulate individual motor unit twitches and their temporal
 # summation to produce the total muscle force.
 
-# Generate force from spike trains
+# Generate force from spike trains, and express it in % MVC with the card's MVC (from its MVC ramp)
 force_output = force_model.generate_force(spike_train__Block=spike_train__Block)
+force_pct_mvc = 100.0 * force_output.magnitude[:, 0] / card.value("drive.mvc_model_units")
 
 # Add some realistic noise to the force signal
 noise_level = 0.015  # 0.15% of mean force
@@ -217,12 +131,12 @@
 
 ax1.plot(
     force_output.times.rescale("s"),
-    force_output[:, 0],
+    force_pct_mvc,
     linewidth=2,
     label="Clean Force",
 )
 
-ax1.set_ylabel("Force (a.u.)")
+ax1.set_ylabel("Force (% MVC)")
 ax1.set_title("Simulated Force Output")
 ax1.grid(True, alpha=0.3)
 
diff --git a/examples/01_basic/README.md b/examples/01_basic/README.md
index 2b269a6f..3e93ee6d 100644
--- a/examples/01_basic/README.md
+++ b/examples/01_basic/README.md
@@ -14,3 +14,44 @@ motor-unit recruitment to EMG signal generation.
 The examples are self-contained and can be run independently, though they may
 depend on each other for context. Following the numbered order provides the
 best learning experience.
+
+## Profiles
+
+Examples 01, 03, 04, 05, 06 and 08 simulate the first dorsal interosseous (FDI)
+with the `fdi_semg` **profile card** (`myogen/profiles/fdi_semg.yaml`): every
+validated setting in one place, each value with its source.
+
+```python
+from myogen import load_profile
+
+card = load_profile()                       # fdi_semg
+run = card.simulate([(1, 0, 0), (4, 0, 25), (4, 25, 25), (4, 25, 0), (1, 0, 0)])  # % MVC
+muscle = card.muscle()
+semg = card.surface_emg(muscle, electrode_arrays)   # your electrodes
+```
+
+- **Override** a value with a keyword, e.g. `card.muscle(fat_thickness__mm=3 * pq.mm)`. The override
+  stays on that card object for every later build (`card.thresholds(n_units=60)` then `card.muscle()`
+  builds 60 units); passing the card value again clears it. Use a fresh `load_profile()` for the card as is.
+- **Other keywords** (electrodes, sampling, units to simulate, ...) pass through to the class and are
+  listed under `keywords` in the built object's `profile_record`.
+- **Notes**: when an overridden value is one the drive calibration or the λ fit
+  depends on, the card prints a note (once per card and field per run) saying what else to check.
+  Silence them with `warnings.simplefilter("ignore", myogen.ProfileNote)`.
+- **Record**: `card.record()` lists the card, its version and every override; built
+  objects carry it as `profile_record`, spike and surface-EMG blocks as the
+  `myogen_profile` annotation.
+- **Seeds**: each builder draws its build from its own seed (the card's `seeds`, or `seed=`; recorded as `seed`),
+  so a build does not depend on what ran before it, and then gives your generator and seed back. This differs from
+  the rest of MyoGen, where library code never seeds and the script owns the seed (one `set_random_seed` call, or
+  the default 180319, for the whole script).
+- **What reproduces**: on any platform, `card.thresholds()` and `card.muscle()` give the validated thresholds and
+  fibres (each fibre's motor unit, diameter and conduction velocity). `card.simulate()` gives the same spikes up
+  to NEURON's floating point, which differs between platforms: spike times can shift by a time step, and a unit can
+  gain or lose a spike near the end of a run (in the library's test run, Linux fires 2001 spikes where macOS and
+  Windows fire 2000). sEMG from a card therefore matches the validation within its limits, not bit for bit.
+- **Sampling**: generate sEMG at 2048 Hz (the `SurfaceEMG` default) and iEMG at 10240 Hz. The card's validation
+  (`benchmarks/profiles/validate_card.py`) records sEMG at 10240 Hz only to score the 20-2000 Hz band of its Delsys
+  comparator.
+- `load_profile("fdi_semg.published")` gives each component as its source paper
+  describes it (for comparison with those papers and older MyoGen results).
diff --git a/examples/02_finetune/06_optimize_dd_multilevel.py b/examples/02_finetune/06_optimize_dd_multilevel.py
index 835330ee..e3e7d8c4 100644
--- a/examples/02_finetune/06_optimize_dd_multilevel.py
+++ b/examples/02_finetune/06_optimize_dd_multilevel.py
@@ -30,7 +30,6 @@
 import numpy as np
 import optuna
 import quantities as pq
-from neo import Block, Segment, SpikeTrain
 from neuron import h
 
 SCRIPT_DIR = Path(__file__).resolve().parent
@@ -39,17 +38,17 @@
 from myogen.simulator import RecruitmentThresholds
 from myogen.simulator.core.force.force_model import ForceModel
 from myogen.simulator.core.force.force_utils import FUGLEVAND_FORCE_FREQUENCY, HUMAN_UNIT_FORCE_FREQUENCY
-from myogen.simulator.neuron import Network
-from myogen.simulator.neuron.populations import AlphaMN__Pool, DescendingDrive__Pool
 from myogen.utils.fr_optimization_targets import compute_fr_targets
-from myogen.utils.helper import get_gamma_shape_for_mvc
 from myogen.utils.nmodl import load_nmodl_mechanisms
 from myogen.utils.ramp_analysis import extract_ramp_recruitment
+from myogen.simulator.neuron.descending_drive_network import (  # noqa: E402  (one authority, shared with the profiles)
+    NETWORK, build_dd_mn_network, excitation_to_drive_hz, run_descending_drive,
+    insert_scaled_gfluctdv as _insert_gfluctdv,
+)
 
 # ─── Constants ────────────────────────────────────────────────
 GAMMA_SHAPE = 3.0
 SYNAPTIC_WEIGHT = 0.05  # µS, DD input weight on the smallest (best-coupled) motor neuron
-NETWORK = {"fan_in": "equal", "weight_scaling": "dendritic_transfer"}  # DD -> MN wiring, see build_dd_mn_network
 WEIGHT_EXPONENT_RANGE = (0.0, 1.0)  # fitted range of the DD weight exponent (--fit-weight-exponent)
 FIT_WEIGHT_EXPONENT = False
 CEILING_WEIGHT = 3.0       # Penalty weight for the recruitment ceiling (force at the last recruitment vs LAST_RECRUIT_PCT_MVC)
@@ -147,81 +146,6 @@ def detect_threshold(thresholds, dd_params, min_fr=5.0):
     return float(threshold_hz)
 
 
-def _insert_gfluctdv(mn_pool, scale):
-    """Insert Gfluctdv independent noise on all MN dendrites at given scale."""
-    if scale <= 0:
-        return
-    for cell in mn_pool:
-        cell.insert_Gfluctdv()
-        for d in cell.dend:
-            d.std_e_Gfluctdv *= scale
-            d.std_i_Gfluctdv *= scale
-            d.g_e0_Gfluctdv *= scale
-            d.g_i0_Gfluctdv *= scale
-
-
-def dendritic_transfer(mn_pool):
-    """Fraction of a steady current injected at each MN's DD synapse site that reaches its clamped soma.
-
-    DC transfer impedance (synapse site -> soma) over soma input impedance, from ``h.Impedance`` linearised at
-    the initial state (active channels included). Call before ``_insert_gfluctdv``.
-    """
-    for s, v in zip(*mn_pool.get_initialization_data()):
-        s.v = v
-    h.finitialize()
-    k = []
-    for cell in mn_pool:
-        syn = next(s for s in cell.synapse__list if s.e >= -40)
-        seg = syn.get_segment()
-        imp = h.Impedance()
-        imp.loc(0.5, sec=cell.soma)
-        imp.compute(0, 1)
-        k.append(imp.transfer(seg.x, sec=seg.sec) / imp.input(0.5, sec=cell.soma))
-    return np.array(k)
-
-
-def build_dd_mn_network(mn_pool, dd_pool, dd_params):
-    """Wire the DD pool to the motor neurons and the cortical input to the DD pool.
-
-    ``dd_params["network"]`` must be ``NETWORK``; parameter files without it (calibrated with random fan-in and
-    one weight for all neurons) are refused rather than reinterpreted.
-
-    - Equal fan-in: every MN receives exactly round(conn_probability * dd_neurons) DD inputs, drawn without
-      replacement. Modelling assumption (descending drive shared across the pool); it removes the binomial
-      +/-13 % input-count scatter that reordered the largest units, whose thresholds are ~3 % apart.
-    - Dendritic-transfer weights: weight_i = synaptic_weight * (k_max / k_i) ** p (``dendritic_transfer``), p from
-      ``network["weight_exponent"]`` (default 1). With p = 1 each input delivers the same current to a clamped soma
-      at low drive whatever the neuron's size, and recruitment order follows rheobase: equal effective synaptic
-      current per input, the equal-share case of Heckman & Binder 1991 (J Neurophysiol 65:952) [heckmanbinder1991]. p < 1 gives larger
-      neurons less effective current per input (p = 0: one weight for all). Modelling assumption either way;
-      measured effective currents are input-system specific (cat data) and none are known for corticospinal
-      input to a hand muscle. Not compensated: the Gfluctdv mean-conductance shunt and the driving-force
-      saturation at high drive.
-    """
-    net_cfg = dict(dd_params.get("network") or {})
-    p = net_cfg.pop("weight_exponent", 1.0)
-    if not isinstance(p, (int, float)) or isinstance(p, bool) or not np.isfinite(p) or p < 0:
-        raise ValueError(f"network weight_exponent must be a finite number >= 0, got {p!r}")
-    if net_cfg != NETWORK:
-        raise ValueError(
-            f"dd_params['network'] must be {NETWORK}, got {net_cfg!r}; parameter files from before 2026-11-01 "
-            "were calibrated with random fan-in and uniform weights and must be recalibrated."
-        )
-    n_dd = len(dd_pool)
-    n_in = int(round(dd_params["conn_probability"] * n_dd))
-    if n_in < 1:
-        raise ValueError(f"conn_probability * dd_neurons rounds to {n_in} inputs per motor neuron")
-    transfer = dendritic_transfer(mn_pool)
-    weights = dd_params["synaptic_weight"] * (transfer.max() / transfer) ** p
-    network = Network({"DD": dd_pool, "aMN": mn_pool})
-    dd_to_mn = network.connect_fan_in("DD", "aMN", n_inputs=n_in, weights__uS=weights * pq.uS,
-                                      delay__ms=1.0 * pq.ms)
-    network.connect_from_external(source="cortical_input", target="DD", weight__uS=1.0 * pq.uS)
-    return SimpleNamespace(network=network, dd_to_mn=dd_to_mn,
-                           cortical_netcons=network.get_netcons("cortical_input", "DD"),
-                           transfer=transfer, weights=weights)
-
-
 def network_block(weight_exponent=None):
     """The ``network`` block of a parameter set: NETWORK, plus the DD weight exponent when it was fitted."""
     return dict(NETWORK) if weight_exponent is None else {**NETWORK, "weight_exponent": float(weight_exponent)}
@@ -296,16 +220,6 @@ def ceiling_penalty(pct_last, target_pct=None):
     return CEILING_WEIGHT * abs(pct_last - target_pct) / target_pct
 
 
-def mvc_drive(segments, threshold_hz, max_hz):
-    """Piecewise-linear excitation [(duration_s, start, end)] and its DD rate (Hz): threshold_hz + e (max - threshold),
-    0.8 threshold_hz at rest, plus 1 % half-normal noise."""
-    parts = [np.linspace(a, b, int(round(d * 1000 / TIMESTEP_MS)), endpoint=False) for d, a, b in segments]
-    e = np.concatenate(parts)
-    hz = np.where(e > 0, threshold_hz + e * (max_hz - threshold_hz), 0.8 * threshold_hz)
-    hz = hz + np.abs(np.random.default_rng(42).normal(0, hz.max() * 0.01, size=hz.size))
-    return e, hz
-
-
 def make_force_model(thresholds):
     """ForceModel (Fuglevand 1993 [fuglevand1993] twitches) at FORCE_FS_HZ with the FORCE_FREQUENCY_SHAPE saturation."""
     return ForceModel(recruitment_thresholds=thresholds, recording_frequency__Hz=FORCE_FS_HZ * pq.Hz,
@@ -341,136 +255,21 @@ def run_mvc_ramp(thresholds, dd_params, threshold_hz):
 
 def run_dd_constant(thresholds, drive_hz, dd_params, duration_ms=3000.0):
     """Run a constant-drive DD simulation (for threshold detection)."""
-    gamma_shape = get_gamma_shape_for_mvc(100, mvc_shape_value=dd_params["gamma_shape"])
-    gfluctdv_scale = dd_params.get("gfluctdv_scale", 0.0)
-    mn_pool = AlphaMN__Pool(
-        recruitment_thresholds__array=thresholds,
-        config_file="alpha_mn_default.yaml",
-    )
-    dd_pool = DescendingDrive__Pool(
-        n=dd_params["dd_neurons"],
-        timestep__ms=TIMESTEP_MS * pq.ms,
-        process_type="gamma",
-        shape=gamma_shape,
-    )
-    wiring = build_dd_mn_network(mn_pool, dd_pool, dd_params)
-    _insert_gfluctdv(mn_pool, gfluctdv_scale)
-    dd_netcons = wiring.cortical_netcons
-
-    mn_recorders = []
-    for cell in mn_pool:
-        rec = h.Vector()
-        nc = h.NetCon(cell.soma(0.5)._ref_v, None, sec=cell.soma)
-        nc.threshold = 50
-        nc.record(rec)
-        mn_recorders.append(rec)
-
-    h.load_file("stdrun.hoc")
-    h.dt = TIMESTEP_MS
-    h.tstop = duration_ms
-
-    for s, v in zip(*mn_pool.get_initialization_data()):
-        s.v = v
-    for s, v in zip(*dd_pool.get_initialization_data()):
-        s.v = v
-    h.finitialize()
-
     n_steps = int(duration_ms / TIMESTEP_MS)
     drive = np.full(n_steps, drive_hz) + np.clip(
         np.random.default_rng(42).normal(0, drive_hz * 0.02, size=n_steps), 0, None,
     )
-
-    step = 0
-    while h.t < h.tstop:
-        d = drive[min(step, n_steps - 1)]
-        for dd_cell in dd_pool:
-            if dd_cell.integrate(d):
-                if h.t < h.tstop:
-                    dd_netcons[dd_cell.pool__ID].event(h.t + 1)
-        h.fadvance()
-        step += 1
-
-    fs = 10240.0
-    t_stop_s = duration_ms / 1000.0
-    block = Block()
-    seg = Segment()
-    block.segments.append(seg)
-    seg.spiketrains = [
-        SpikeTrain(
-            rec.as_numpy() / 1000 * pq.s,
-            t_stop=t_stop_s * pq.s,
-            sampling_rate=fs * pq.Hz,
-            sampling_period=(1 / fs) * pq.s,
-            name=str(i),
-        )
-        for i, rec in enumerate(mn_recorders)
-    ]
-    return block
+    return run_descending_drive(thresholds, drive, dd_params, duration_ms, timestep_ms=TIMESTEP_MS)
 
 
 def run_dd_ramp(thresholds, drive_signal, dd_params, total_ms):
-    """Run a DD simulation with time-varying drive (trapezoidal ramp)."""
-    gamma_shape = get_gamma_shape_for_mvc(100, mvc_shape_value=dd_params["gamma_shape"])
-    gfluctdv_scale = dd_params.get("gfluctdv_scale", 0.0)
-    mn_pool = AlphaMN__Pool(
-        recruitment_thresholds__array=thresholds,
-        config_file="alpha_mn_default.yaml",
-    )
-    dd_pool = DescendingDrive__Pool(
-        n=dd_params["dd_neurons"],
-        timestep__ms=TIMESTEP_MS * pq.ms,
-        process_type="gamma",
-        shape=gamma_shape,
-    )
-    wiring = build_dd_mn_network(mn_pool, dd_pool, dd_params)
-    _insert_gfluctdv(mn_pool, gfluctdv_scale)
-    dd_netcons = wiring.cortical_netcons
-
-    mn_recorders = []
-    for cell in mn_pool:
-        rec = h.Vector()
-        nc = h.NetCon(cell.soma(0.5)._ref_v, None, sec=cell.soma)
-        nc.threshold = 50
-        nc.record(rec)
-        mn_recorders.append(rec)
-
-    h.load_file("stdrun.hoc")
-    h.dt = TIMESTEP_MS
-    h.tstop = total_ms
-
-    for s, v in zip(*mn_pool.get_initialization_data()):
-        s.v = v
-    for s, v in zip(*dd_pool.get_initialization_data()):
-        s.v = v
-    h.finitialize()
-
-    n_steps = len(drive_signal)
-    step = 0
-    while h.t < h.tstop:
-        d = drive_signal[min(step, n_steps - 1)]
-        for dd_cell in dd_pool:
-            if dd_cell.integrate(d):
-                if h.t < h.tstop:
-                    dd_netcons[dd_cell.pool__ID].event(h.t + 1)
-        h.fadvance()
-        step += 1
-
-    fs = 10240.0
-    t_stop_s = total_ms / 1000.0
-    block = Block()
-    seg = Segment()
-    block.segments.append(seg)
-    seg.spiketrains = [
-        SpikeTrain(
-            rec.as_numpy() / 1000 * pq.s,
-            t_stop=t_stop_s * pq.s,
-            sampling_rate=fs * pq.Hz,
-            sampling_period=(1 / fs) * pq.s,
-            name=str(i),
-        )
-        for i, rec in enumerate(mn_recorders)
-    ]
-    return block
+    """Run a DD simulation with time-varying drive (``run_descending_drive`` at TIMESTEP_MS)."""
+    return run_descending_drive(thresholds, drive_signal, dd_params, total_ms, timestep_ms=TIMESTEP_MS)
+
+
+def mvc_drive(segments, threshold_hz, max_hz):
+    """``excitation_to_drive_hz`` at TIMESTEP_MS (see there)."""
+    return excitation_to_drive_hz(segments, threshold_hz, max_hz, timestep_ms=TIMESTEP_MS)
 
 
 def make_objective(thresholds, fug_targets):
diff --git a/myogen/__init__.py b/myogen/__init__.py
index 785ce613..5b26bf53 100644
--- a/myogen/__init__.py
+++ b/myogen/__init__.py
@@ -1,4 +1,6 @@
 import warnings
+from contextlib import contextmanager
+from copy import deepcopy
 
 from myogen._cuda_env import setup as _setup_cuda
 _setup_cuda()
@@ -91,6 +93,36 @@ def set_random_seed(seed: int = _DEFAULT_SEED) -> None:
     print(f"Random seed set to {seed}.")
 
 
+@contextmanager
+def _seed_scope(seed: int, state: dict | None = None, reset_cell_ids: bool = True):
+    """Draw from a generator seeded with ``seed`` (set to ``state`` if given) inside the block, then give the
+    caller's generator and seed back.
+
+    Used by the profile-card builders (``myogen.profiles``), the one place where library code seeds: a card replays
+    the per-stage seeds of the runs it was validated with. Unlike ``set_random_seed`` it prints nothing and leaves
+    the caller's random stream (``get_random_generator``, ``get_random_seed``, ``derive_subseed``) as it was. The
+    cell ID counters are reset at the start, as ``set_random_seed`` does, and not restored; an error while
+    resetting them is raised after the caller's generator is given back.
+    """
+    global _random_generator, _current_seed
+    saved = _random_generator, _current_seed
+    _current_seed = int(seed)
+    _random_generator = np.random.default_rng(_current_seed)
+    if state is not None:
+        _random_generator.bit_generator.state = deepcopy(state)
+    try:
+        if reset_cell_ids:
+            try:
+                from myogen.simulator.neuron.cells import reset_cell_id_counters
+            except ImportError:  # NEURON not installed: there are no cell counters to reset
+                pass
+            else:
+                reset_cell_id_counters()
+        yield _random_generator
+    finally:
+        _random_generator, _current_seed = saved
+
+
 def __getattr__(name: str):
     """Backwards-compatible access for the deprecated ``RANDOM_GENERATOR`` and ``SEED`` module attributes."""
     if name == "RANDOM_GENERATOR":
@@ -305,7 +337,12 @@ def error(msg):
 # Use strict=True in your code to raise exceptions on failure
 _nmodl_loaded = load_nmodl_mechanisms(quiet=True, strict=False)
 
+from myogen.profiles import ProfileError, ProfileNote, load_profile  # noqa: E402
+
 __all__ = [
+    "load_profile",
+    "ProfileError",
+    "ProfileNote",
     "get_random_generator",
     "get_random_seed",
     "derive_subseed",
diff --git a/myogen/profiles/__init__.py b/myogen/profiles/__init__.py
new file mode 100644
index 00000000..3c6835e3
--- /dev/null
+++ b/myogen/profiles/__init__.py
@@ -0,0 +1,11 @@
+"""Named simulation profiles ("cards"): every validated setting of a simulation in one place.
+
+``load_profile()`` gives the ``fdi_semg`` card (FDI surface EMG); ``load_profile("fdi_semg.published")`` its
+variant with each component as its source paper describes it. The card builds the library objects and runs the
+calibrated drive; keywords passed to a builder override the card, and a note is printed when an overridden value
+is one the card's calibration or fit depends on.
+"""
+from myogen.profiles._card import CardField, Profile, ProfileError, available_profiles, load_profile
+from myogen.profiles._notes import ProfileNote
+
+__all__ = ["CardField", "Profile", "ProfileError", "ProfileNote", "available_profiles", "load_profile"]
diff --git a/myogen/profiles/_builders.py b/myogen/profiles/_builders.py
new file mode 100644
index 00000000..4214aef6
--- /dev/null
+++ b/myogen/profiles/_builders.py
@@ -0,0 +1,416 @@
+"""Card builders: the library objects a card describes, with keyword overrides noted and recorded.
+
+Seeds. Each builder that draws takes a ``seed`` (default: the card's, from ``seeds``; recorded as ``seed`` in the
+built object's ``profile_record``) and draws its build from a generator seeded with it, so a build does not depend on
+what ran before it, and reproduces the per-stage seeds of the runs the card was validated with. Afterwards the
+caller's generator and seed are given back (``myogen._seed_scope``): the user's ``set_random_seed`` stream is as it
+was, and nothing is printed. This deviates from base MyoGen, where library code never seeds and the script owns
+the seed (one ``set_random_seed`` call, or the import-time default, for the whole script). The cell ID counters are
+reset at the start of a seeded build, as ``set_random_seed`` does. On any platform the seeded builds give the
+validated thresholds and fibres; ``simulate`` gives the same spikes up to NEURON's floating point, which differs
+between platforms (spike times can shift by a time step, and a unit can gain or lose a spike near the end of a run),
+so sEMG matches the validation within its limits.
+A keyword that names a card field overrides it on this card instance: the override applies to
+every later build from the same instance (``thresholds(n_units=60)`` then ``muscle()`` builds 60 units), is recorded
+(``overrides()``, ``record()``) and, when the field carries a note, prints the note once per card and field per run
+(``ProfileNote``).
+A keyword equal to the card value (unit-aware) is not an override, and clears an earlier one. Constructor keywords
+the card does not set (electrodes, sampling, units to simulate, ...) are passed through and listed under
+``keywords`` in the built object's ``profile_record``.
+"""
+from __future__ import annotations
+
+import copy
+import inspect
+from dataclasses import asdict, dataclass, is_dataclass
+from typing import Any
+
+import numpy as np
+import quantities as pq
+
+from myogen.profiles._notes import emit_note
+
+# Card field -> constructor keyword, where they differ
+_THRESHOLD_KEYWORDS = {"n_units": "N", "recruitment_range": "recruitment_range__ratio", "mode": "mode"}
+_UNIT_CV_KEYWORD = "motor_unit_conduction_velocities__m_per_s"
+_PER_UNIT = "per_unit__m_per_s"  # muscle.unit_cv override holding explicit per-unit CVs
+
+
+def _magnitude(value: Any, unit: str | None) -> Any:
+    if isinstance(value, pq.Quantity):
+        return value.rescale(pq.Quantity(1.0, unit)).magnitude if unit else value.magnitude
+    return value
+
+
+def _same(a: Any, b: Any) -> bool:
+    if a is None or b is None:
+        return a is None and b is None
+    if isinstance(a, dict) and isinstance(b, dict):
+        return {str(k) for k in a} == {str(k) for k in b} and all(
+            _same(a[k], next(b[j] for j in b if str(j) == str(k))) for k in a)
+    if isinstance(a, str) or isinstance(b, str):
+        return a == b
+    try:
+        x, y = np.asarray(a, dtype=float), np.asarray(b, dtype=float)
+    except (TypeError, ValueError):
+        return a == b
+    return x.shape == y.shape and bool(np.allclose(x, y, rtol=1e-12, atol=0.0))
+
+
+def _jsonable(value: Any) -> Any:
+    if isinstance(value, np.ndarray):
+        return value.tolist()
+    if isinstance(value, np.generic):
+        return value.item()
+    if isinstance(value, tuple):
+        return list(value)
+    return value
+
+
+def _describe(value: Any) -> Any:
+    """A JSON-friendly description of a passed-through keyword for the record."""
+    if isinstance(value, pq.Quantity):
+        mag = _jsonable(np.asarray(value.magnitude))
+        return f"{mag} {value.dimensionality.string}"
+    if is_dataclass(value):
+        return {k: _jsonable(v) for k, v in asdict(value).items()}
+    if isinstance(value, (str, int, float, bool, type(None), np.ndarray, np.generic, list, tuple, dict)):
+        return _jsonable(value)
+    return repr(value)
+
+
+def _dataclass_from(kind: type, spec: dict) -> Any:
+    return kind(**{k: tuple(x) if isinstance(x, list) else x for k, x in spec.items()})
+
+
+class CardBuilders:
+    """Builder methods of ``Profile`` (mixed in; uses ``field``, ``value``, ``seeds``, ``name`` and ``version``)."""
+
+    # -- overrides --------------------------------------------------------------------------------------------
+    def overrides(self) -> dict[str, Any]:
+        """Every override in force on this card instance, as ``{"section.field": value}`` (card units)."""
+        return copy.deepcopy(self._overrides)
+
+    def record(self, keywords: dict | None = None, seed: int | None = None) -> dict[str, Any]:
+        """What a result was built with: card name, version, parent, overrides (and passed-through keywords, seed)."""
+        rec = {"profile": self.name, "version": self.version, "parent": self.parent,
+               "parent_version": self.parent_version, "overrides": self.overrides()}
+        if keywords is not None:
+            rec["keywords"] = dict(keywords)
+        if seed is not None:
+            rec["seed"] = seed
+        return rec
+
+    def _seed(self, stage: str, seed: int | None) -> int:
+        """The seed of a build: ``seed`` if given, else the card's seed for ``stage``."""
+        if seed is None:
+            return int(self.seeds[stage])
+        if isinstance(seed, bool) or not isinstance(seed, (int, np.integer)):
+            raise TypeError(f"{self.name}: seed must be an integer, got {seed!r}")
+        return int(seed)
+
+    def _card_form(self, path: str, field, key: str, given: Any) -> Any:
+        """``given`` in the card's format and units (what the card value would be)."""
+        if key == _UNIT_CV_KEYWORD:
+            return {_PER_UNIT: np.asarray(_magnitude(given, "m/s"), dtype=float).tolist()}
+        if is_dataclass(given):
+            return {k: _jsonable(v) for k, v in asdict(given).items()}
+        if isinstance(given, pq.Quantity):
+            if field.unit is None:
+                raise TypeError(f"{self.name}: {path} has no unit; got the quantity {given!r}")
+            return _jsonable(_magnitude(given, field.unit))
+        return _jsonable(given)
+
+    def _card_value_in_form_of(self, path: str, key: str, given: Any) -> Any:
+        """The card value, in the same form as ``_card_form`` gives for ``given``."""
+        value = self.value(path)
+        if key == _UNIT_CV_KEYWORD:
+            cvs = self._unit_cvs(value)
+            return None if cvs is None else {_PER_UNIT: cvs.tolist()}
+        if is_dataclass(given) and isinstance(value, dict):
+            return {k: _jsonable(v) for k, v in asdict(_dataclass_from(type(given), value)).items()}
+        return value
+
+    def _take(self, section: str, overrides: dict, allowed: set[str], aliases: dict[str, str] | None = None,
+              keywords: dict | None = None) -> dict:
+        """Card values of ``section`` with this instance's overrides and ``overrides`` applied (card units).
+
+        Keywords naming a card field (or an alias of one) override it; keywords in ``allowed`` pass through and are
+        described in ``keywords``; anything else raises ``TypeError``.
+        """
+        aliases = aliases or {}
+        fields = self.section(section)
+        out = {name: self.value(f"{section}.{name}") for name in fields}
+        for path, val in self._overrides.items():
+            sec, name = path.split(".", 1)
+            if sec == section:
+                out[name] = copy.deepcopy(val)
+        for key, given in overrides.items():
+            name = aliases.get(key, key)
+            path = f"{section}.{name}"
+            if name in fields:
+                f = self.field(path)
+                plain = self._card_form(path, f, key, given)
+                if _same(plain, self._card_value_in_form_of(path, key, given)):
+                    self._overrides.pop(path, None)  # back to the card value
+                    out[name] = self.value(path)
+                    continue
+                out[name] = copy.deepcopy(plain)
+                self._override(path, f, plain)
+            elif key in allowed:  # a constructor keyword the card does not set: the user's choice
+                out[key] = given
+                if keywords is not None:
+                    keywords[key] = _describe(given)
+            else:
+                raise TypeError(f"{self.name}: unknown keyword {key!r} for {section} "
+                                f"(card fields: {sorted(fields)})")
+        return out
+
+    def _override(self, path: str, field, value: Any) -> None:
+        others = sorted(p for p in self._overrides if p != path)
+        self._overrides[path] = copy.deepcopy(value)
+        note = field.note
+        if not note and field.calibrated:
+            note = (f"The descending-drive calibration (trial {self.calibration.get('best_trial', '?')}) was fitted "
+                    "with this value, so the drive calibration no longer holds.")
+        if note:
+            text = f"{self.name}: {path} overridden ({value!r}). {note}"
+            if others:
+                text += f" You also changed: {', '.join(others)}."
+            emit_note(f"{self.name}:{path}", text)
+
+    def _q(self, section: str, name: str, value: Any) -> Any:
+        """Attach the card unit to a plain value."""
+        unit = self.field(f"{section}.{name}").unit
+        if unit is None or value is None or isinstance(value, pq.Quantity):
+            return value
+        return value * pq.Quantity(1.0, unit)
+
+    def _unit_cvs(self, spec: Any) -> np.ndarray | None:
+        """Per-unit fibre CVs (m/s) for a ``muscle.unit_cv`` value, or None for surface CVs."""
+        from myogen.utils.ramp_analysis import cv_from_recruitment_force
+
+        if spec is None:
+            return None
+        if _PER_UNIT in spec:
+            return np.asarray(spec[_PER_UNIT], dtype=float)
+        if spec.get("from") != "recruitment_force":
+            raise ValueError(f"{self.name}: muscle.unit_cv supports from: recruitment_force, got {spec!r}")
+        rt = self._take("drive", {}, set())["recruitment_pct_mvc"]
+        cv, _ = cv_from_recruitment_force(np.asarray(rt, dtype=float),
+                                          residual_sd__m_s=float(spec["residual_sd__m_s"]),
+                                          rng=np.random.default_rng(int(self.seeds["cv_scatter"])))
+        return np.asarray(cv, dtype=float)
+
+    # -- builders ---------------------------------------------------------------------------------------------
+    def thresholds(self, seed: int | None = None, **overrides) -> np.ndarray:
+        """Recruitment thresholds, normalised to a largest value of 1."""
+        from myogen import _seed_scope
+        from myogen.simulator import RecruitmentThresholds
+
+        seed = self._seed("thresholds", seed)
+        v = self._take("thresholds", overrides, set())
+        with _seed_scope(seed):
+            thr, _ = RecruitmentThresholds(**{_THRESHOLD_KEYWORDS[k]: v[k] for k in _THRESHOLD_KEYWORDS})
+        thr = np.asarray(thr, dtype=float)
+        return thr / thr.max()
+
+    def muscle_kwargs(self, **overrides) -> dict[str, Any]:
+        """The ``Muscle(...)`` keywords of the card (nothing is built)."""
+        return self._muscle_kwargs(overrides, {})
+
+    def _muscle_kwargs(self, overrides: dict, keywords: dict) -> dict[str, Any]:
+        from myogen.simulator import Muscle
+
+        allowed = set(inspect.signature(Muscle.__init__).parameters) - {"self"}
+        v = self._take("muscle", overrides, allowed, aliases={_UNIT_CV_KEYWORD: "unit_cv"}, keywords=keywords)
+        kw = {k: val for k, val in v.items() if k in allowed}
+        for name in ("radius__mm", "fiber_density__fibers_per_mm2", "skin_thickness__mm", "fat_thickness__mm",
+                     "skin_conductivity__S_per_m"):
+            kw[name] = self._q("muscle", name, v[name])
+        log_sd = v["within_unit_fiber_diameter_log_sd"]
+        kw["within_unit_fiber_diameter_log_sd"] = None if log_sd is None else tuple(log_sd)
+        cvs = self._unit_cvs(v["unit_cv"])
+        if cvs is not None:
+            kw[_UNIT_CV_KEYWORD] = cvs
+        kw.setdefault("recruitment_thresholds", self.thresholds())
+        return kw
+
+    def muscle(self, n_jobs: int = -2, seed: int | None = None, **overrides):
+        """The card's ``Muscle``, built (innervation centres, fibres, assignment, fibre properties).
+
+        Fibre properties are generated only when the card's fibre dispersion uses Muscle fibres
+        (``fiber_properties: muscle``), with the global generator's state restored afterwards, as in the benchmark
+        the card was validated with, so every later draw is the same with or without them.
+        """
+        from myogen import _seed_scope
+        from myogen.simulator import Muscle
+
+        seed = self._seed("muscle", seed)
+        keywords: dict = {}
+        kw = self._muscle_kwargs(overrides, keywords)
+        disp = self._take("surface_emg", {}, set())["fiber_dispersion"]
+        with _seed_scope(seed) as rng:
+            m = Muscle(**kw)
+            m.distribute_innervation_centers()
+            m.generate_muscle_fiber_centers(verbose=False)
+            m.assign_mfs2mns(verbose=False, n_jobs=n_jobs)
+            if isinstance(disp, dict) and disp.get("fiber_properties") == "muscle":
+                state = rng.bit_generator.state
+                m._generate_fiber_properties()
+                rng.bit_generator.state = state
+        m.profile_record = self.record(keywords, seed=seed)
+        return m
+
+    def force_model(self, **overrides):
+        """The card's ``ForceModel`` (twitches and force-frequency curve the drive was calibrated with)."""
+        from myogen.simulator.core.force.force_model import ForceModel
+
+        allowed = set(inspect.signature(ForceModel.__init__).parameters) - {"self"}
+        keywords: dict = {}
+        v = self._take("force", overrides, allowed, keywords=keywords)
+        kw = {k: val for k, val in v.items() if k in allowed and k not in ("force_frequency_shape",
+                                                                             "recording_frequency__Hz")}
+        kw.setdefault("recruitment_thresholds", self.thresholds())
+        fm = ForceModel(recording_frequency__Hz=self._q("force", "recording_frequency__Hz", v["recording_frequency__Hz"]),
+                        force_frequency_shape=tuple(float(x) for x in v["force_frequency_shape"]), **kw)
+        fm.profile_record = self.record(keywords)
+        return fm
+
+    def dd_params(self, **overrides) -> dict[str, Any]:
+        """Descending-drive parameters in the calibration's format (``run_descending_drive``).
+
+        Keywords may name fields of the drive or the motor_neurons section.
+        """
+        mn_fields = set(self.section("motor_neurons"))
+        drive = self._take("drive", {k: v for k, v in overrides.items() if k not in mn_fields}, set())
+        mn = self._take("motor_neurons", {k: v for k, v in overrides.items() if k in mn_fields}, set())
+        p = {k: drive[k] for k in ("dd_neurons", "conn_probability", "dd_drive__Hz", "synaptic_weight", "gamma_shape",
+                                   "network")}
+        p["gfluctdv_scale"] = float(mn["gfluctdv_scale"])
+        return p
+
+    def motor_neuron_pool(self, seed: int | None = None, **overrides):
+        """Motor neuron pool, descending-drive pool and their calibrated wiring (``mn_pool``, ``dd_pool``, ``wiring``)."""
+        from myogen import _seed_scope
+        from myogen.simulator.neuron.descending_drive_network import build_pools
+
+        seed = self._seed("drive", seed)
+        params = self.dd_params(**overrides)
+        mn = self._take("motor_neurons", {}, set())
+        drive = self._take("drive", {}, set())
+        thr = self.thresholds()
+        with _seed_scope(seed):
+            pools = build_pools(thr, params, float(drive["timestep_ms"]), mn["config_file"])
+        pools.profile_record = self.record(seed=seed)
+        return pools
+
+    def surface_emg_kwargs(self, **overrides) -> dict[str, Any]:
+        """The card's ``SurfaceEMG`` keywords besides the muscle and electrodes (nothing is built)."""
+        return self._surface_emg_kwargs(overrides, {})
+
+    def _surface_emg_kwargs(self, overrides: dict, keywords: dict) -> dict[str, Any]:
+        from myogen.simulator import SurfaceEMG
+        from myogen.simulator.core.emg.surface.fiber_dispersion import SurfaceFiberDispersion
+
+        allowed = set(inspect.signature(SurfaceEMG.__init__).parameters) - {"self", "muscle_model", "electrode_arrays"}
+        v = self._take("surface_emg", overrides, allowed, keywords=keywords)
+        kw = {k: val for k, val in v.items() if k in allowed}
+        disp = v["fiber_dispersion"]
+        kw["fiber_dispersion"] = _dataclass_from(SurfaceFiberDispersion, disp) if isinstance(disp, dict) else disp
+        kw["iap_lambda__per_mm"] = float(v["iap_lambda__per_mm"])
+        kw["iap_kernel_length__mm"] = float(v["iap_kernel_length__mm"])
+        return kw
+
+    def surface_emg(self, muscle, electrode_arrays, seed: int | None = None, **overrides):
+        """The card's ``SurfaceEMG`` for ``muscle`` (from ``muscle()``) and the user's electrode arrays.
+
+        The seed and the generator state after construction are kept on the object (``profile_seed``,
+        ``profile_rng_state``); ``simulate_muaps`` draws from them and gives the caller's generator back, so the MUAPs
+        do not depend on what ran between building and simulating.
+        """
+        from myogen import _seed_scope
+        from myogen.simulator import SurfaceEMG
+
+        seed = self._seed("surface_emg", seed)
+        keywords: dict = {}
+        kw = self._surface_emg_kwargs(overrides, keywords)
+        with _seed_scope(seed) as rng:
+            semg = SurfaceEMG(muscle, electrode_arrays, **kw)
+            semg.profile_rng_state = copy.deepcopy(rng.bit_generator.state)
+        semg.profile_seed = seed
+        semg.profile_record = self.record(keywords, seed=seed)
+        return semg
+
+    # -- running in % MVC -------------------------------------------------------------------------------------
+    def excitation_for_pct_mvc(self, pct):
+        """Excitation (0-1) for a force level in % MVC, interpolated on the card's MVC-ramp curve.
+
+        The curve is extended to rest (0 % MVC, excitation 0) and MVC (100 %, excitation 1) when it stops short.
+        """
+        from myogen.profiles._card import ProfileError
+
+        curve = self._take("drive", {}, set())["excitation_at_pct_mvc"]
+        if curve is None:
+            raise ProfileError(f"{self.name}: drive.excitation_at_pct_mvc is not filled in, so % MVC levels cannot be "
+                               "mapped to excitation (run benchmarks/profiles/validate_card.py --mvc-mapping)")
+        pct_arr = np.asarray(pct, dtype=float)
+        if not np.all(np.isfinite(pct_arr)) or np.any(pct_arr < 0) or np.any(pct_arr > 100):
+            raise ValueError(f"force levels must be within 0-100 % MVC, got {pct!r}")
+        points = {float(k): float(v) for k, v in curve.items()}
+        points.setdefault(0.0, 0.0)
+        points.setdefault(100.0, 1.0)
+        levels = sorted(points)
+        e = np.interp(pct_arr, levels, [points[k] for k in levels])
+        return float(e) if np.ndim(e) == 0 else e
+
+    def simulate(self, protocol, seed: int | None = None, **overrides) -> "SimulationResult":
+        """Run the calibrated drive for ``protocol``, segments ``(duration_s, start_pct_mvc, end_pct_mvc)``.
+
+        Each segment's ends are mapped to excitation with ``excitation_for_pct_mvc`` and the excitation is linear in
+        between, as in the calibration. Returns the motor-unit spike trains, force in % MVC and the excitation.
+        Keywords override card fields of the thresholds, motor_neurons, drive and force sections; ``seed`` (default:
+        the card's drive seed) seeds the drive noise and the network.
+        """
+        from myogen import _seed_scope
+        from myogen.simulator.neuron import descending_drive_network as ddn
+
+        seed = self._seed("drive", seed)
+
+        by_section = {s: {} for s in ("thresholds", "motor_neurons", "drive", "force")}
+        for key, val in overrides.items():
+            owner = next((s for s in by_section if key in self.section(s)), None)
+            if owner is None:
+                raise TypeError(f"{self.name}: unknown keyword {key!r} for simulate "
+                                f"(fields of {', '.join(by_section)})")
+            by_section[owner][key] = val
+        thr = self.thresholds(**by_section["thresholds"])
+        params = self.dd_params(**by_section["drive"], **by_section["motor_neurons"])
+        self._take("force", by_section["force"], set())
+        drive = self._take("drive", {}, set())
+        mn = self._take("motor_neurons", {}, set())
+        segments = [(float(d), self.excitation_for_pct_mvc(a), self.excitation_for_pct_mvc(b)) for d, a, b in protocol]
+        dt = float(drive["timestep_ms"])
+        e, hz = ddn.excitation_to_drive_hz(segments, float(drive["threshold_hz"]), params["dd_drive__Hz"],
+                                           timestep_ms=dt, noise_seed=seed)
+        with _seed_scope(seed):
+            block = ddn.run_descending_drive(thr, hz, params, len(e) * dt, timestep_ms=dt,
+                                             config_file=mn["config_file"])
+        fm = self.force_model()
+        force = np.asarray(fm.generate_force(block, verbose=False).magnitude, dtype=float)[:, 0]
+        record = self.record(seed=seed)
+        block.annotations["myogen_profile"] = record
+        fs = float(self._take("force", {}, set())["recording_frequency__Hz"])
+        return SimulationResult(spikes=block, force_pct_mvc=100.0 * force / float(drive["mvc_model_units"]),
+                                excitation=e, fs_force__Hz=fs, profile=record)
+
+
+@dataclass
+class SimulationResult:
+    """Spike trains (Neo block, annotated with the card record), force in % MVC at ``fs_force__Hz``, excitation per
+    drive time step, and the card record."""
+    spikes: Any
+    force_pct_mvc: np.ndarray
+    excitation: np.ndarray
+    fs_force__Hz: float
+    profile: dict
diff --git a/myogen/profiles/_card.py b/myogen/profiles/_card.py
new file mode 100644
index 00000000..71d1cb74
--- /dev/null
+++ b/myogen/profiles/_card.py
@@ -0,0 +1,170 @@
+"""Profile cards: YAML files holding every validated setting of one simulation, with sources and notes.
+
+A card is a mapping of sections (``thresholds``, ``motor_neurons``, ``drive``, ``force``, ``muscle``,
+``surface_emg``) to fields. Each field has a ``value``, a ``source`` (a literature-register key with a locator, the
+word ``assumption``, or ``derived: ...``), and optionally a ``unit``, a ``note`` (printed when the value is
+overridden) and ``calibrated: true`` (the descending-drive calibration was fitted with this value). A variant names
+its ``parent`` and lists only the fields it changes.
+"""
+from __future__ import annotations
+
+import copy
+from dataclasses import dataclass
+from pathlib import Path
+from types import MappingProxyType
+from typing import Any
+
+import yaml
+
+from myogen.profiles._builders import CardBuilders
+
+CARD_DIR = Path(__file__).resolve().parent
+DEFAULT_PROFILE = "fdi_semg"
+SECTIONS = ("thresholds", "motor_neurons", "drive", "force", "muscle", "surface_emg")
+FIELD_KEYS = {"value", "source", "unit", "note", "calibrated"}
+# Fields and seeds the builders read: a card (with its parent) must have all of them
+REQUIRED_FIELDS = {
+    "thresholds": ("n_units", "mode", "recruitment_range"),
+    "motor_neurons": ("config_file", "gfluctdv_scale"),
+    "drive": ("dd_neurons", "conn_probability", "dd_drive__Hz", "synaptic_weight", "gamma_shape", "network",
+              "threshold_hz", "timestep_ms", "mvc_model_units", "excitation_at_pct_mvc", "recruitment_pct_mvc"),
+    "force": ("force_frequency_shape", "recording_frequency__Hz"),
+    "muscle": ("radius__mm", "fiber_density__fibers_per_mm2", "skin_thickness__mm", "fat_thickness__mm",
+               "skin_conductivity__S_per_m", "within_unit_fiber_diameter_log_sd", "unit_cv"),
+    "surface_emg": ("fiber_dispersion", "iap_lambda__per_mm", "iap_kernel_length__mm"),
+}
+REQUIRED_SEEDS = ("thresholds", "muscle", "cv_scatter", "drive", "surface_emg")
+
+
+class ProfileError(ValueError):
+    """A card is missing, malformed, or used in a way it does not support."""
+
+
+@dataclass(frozen=True)
+class CardField:
+    value: Any
+    source: str
+    unit: str | None = None
+    note: str | None = None
+    calibrated: bool = False
+
+
+def available_profiles() -> list[str]:
+    """Names of the cards shipped with MyoGen."""
+    return sorted(p.stem for p in CARD_DIR.glob("*.yaml"))
+
+
+def _read(name_or_path: str | Path) -> tuple[dict, Path]:
+    path = Path(name_or_path)
+    if path.suffix not in (".yaml", ".yml"):
+        path = CARD_DIR / f"{name_or_path}.yaml"
+    if not path.is_file():
+        raise ProfileError(f"no profile {str(name_or_path)!r}; available: {', '.join(available_profiles())}")
+    with open(path, encoding="utf-8") as fh:
+        raw = yaml.safe_load(fh)
+    if not isinstance(raw, dict):
+        raise ProfileError(f"{path.name}: a card must be a YAML mapping")
+    return raw, path
+
+
+def _field(path: str, spec: Any) -> CardField:
+    if not isinstance(spec, dict):
+        raise ProfileError(f"{path}: a field must be a mapping with 'value' and 'source'")
+    unknown = set(spec) - FIELD_KEYS
+    if unknown:
+        raise ProfileError(f"{path}: unknown keys {sorted(unknown)} (allowed: {sorted(FIELD_KEYS)})")
+    for key in ("value", "source"):
+        if key not in spec:
+            raise ProfileError(f"{path}: missing '{key}'")
+    if not isinstance(spec["source"], str) or not spec["source"].strip():
+        raise ProfileError(f"{path}: 'source' must be a register key with a locator, 'assumption' or 'derived: ...'")
+    return CardField(value=spec["value"], source=spec["source"], unit=spec.get("unit"), note=spec.get("note"),
+                     calibrated=bool(spec.get("calibrated", False)))
+
+
+def _fields(where: str, sections: Any) -> dict[str, CardField]:
+    if not isinstance(sections, dict):
+        raise ProfileError(f"{where}: 'sections' must be a mapping of {', '.join(SECTIONS)}")
+    out = {}
+    for section, fields in sections.items():
+        if section not in SECTIONS:
+            raise ProfileError(f"{where}: unknown section {section!r} (allowed: {', '.join(SECTIONS)})")
+        if not isinstance(fields, dict):
+            raise ProfileError(f"{where}: section {section!r} must be a mapping of fields")
+        for name, spec in fields.items():
+            out[f"{section}.{name}"] = _field(f"{section}.{name}", spec)
+    return out
+
+
+class Profile(CardBuilders):
+    """A loaded card: its fields (read-only), the builders, and the overrides applied on this instance."""
+
+    def __init__(self, name: str, version: int, description: str, parent: str | None, seeds: dict,
+                 fields: dict[str, CardField], calibration: dict, validation: dict | None, path: Path,
+                 parent_version: int | None = None):
+        self.name = name
+        self.version = version
+        self.description = description
+        self.parent = parent
+        self.parent_version = parent_version
+        self.seeds = MappingProxyType(dict(seeds))
+        self._fields = MappingProxyType(dict(fields))
+        self.calibration = MappingProxyType(dict(calibration))
+        self.validation = None if validation is None else MappingProxyType(dict(validation))
+        self.path = path
+        self._overrides: dict[str, Any] = {}
+
+    def __repr__(self) -> str:
+        return f"Profile({self.name!r}, version={self.version})"
+
+    def fields(self) -> dict[str, CardField]:
+        return dict(self._fields)
+
+    def field(self, path: str) -> CardField:
+        try:
+            return self._fields[path]
+        except KeyError:
+            raise ProfileError(f"{self.name} has no field {path!r}") from None
+
+    def value(self, path: str) -> Any:
+        """The card value of ``section.field`` (a copy, so callers cannot change the card)."""
+        return copy.deepcopy(self.field(path).value)
+
+    def section(self, section: str) -> dict[str, CardField]:
+        prefix = f"{section}."
+        return {k[len(prefix):]: v for k, v in self._fields.items() if k.startswith(prefix)}
+
+
+def load_profile(name: str | Path = DEFAULT_PROFILE) -> Profile:
+    """Load a profile card by name (``"fdi_semg"``, ``"fdi_semg.published"``) or from a YAML path."""
+    raw, path = _read(name)
+    where = path.name
+    for key in ("name", "version"):
+        if key not in raw:
+            raise ProfileError(f"{where}: missing top-level '{key}'")
+    if "sections" not in raw:
+        raise ProfileError(f"{where}: missing top-level 'sections'")
+    own = _fields(where, raw["sections"])
+    parent_name = raw.get("parent")
+    parent_version = None
+    if parent_name is None:
+        if "seeds" not in raw:
+            raise ProfileError(f"{where}: missing top-level 'seeds'")
+        fields, seeds, calibration = own, dict(raw["seeds"] or {}), dict(raw.get("calibration") or {})
+    else:
+        parent = load_profile(parent_name)
+        parent_version = parent.version
+        fields = parent.fields()
+        for key, f in own.items():
+            if key not in fields:
+                raise ProfileError(f"{where}: {key} is not a field of the parent card {parent_name!r}")
+            fields[key] = f
+        seeds = {**parent.seeds, **(raw.get("seeds") or {})}
+        calibration = {**parent.calibration, **(raw.get("calibration") or {})}
+    missing = [f"{s}.{n}" for s, names in REQUIRED_FIELDS.items() for n in names if f"{s}.{n}" not in fields]
+    missing += [f"seeds.{k}" for k in REQUIRED_SEEDS if k not in seeds]
+    if missing:
+        raise ProfileError(f"{where}: missing required {', '.join(missing)}")
+    return Profile(name=str(raw["name"]), version=int(raw["version"]), description=str(raw.get("description", "")),
+                   parent=parent_name, seeds=seeds, fields=fields, calibration=calibration,
+                   validation=raw.get("validation"), path=path, parent_version=parent_version)
diff --git a/myogen/profiles/_notes.py b/myogen/profiles/_notes.py
new file mode 100644
index 00000000..03731c42
--- /dev/null
+++ b/myogen/profiles/_notes.py
@@ -0,0 +1,37 @@
+"""Notes on the back of a profile card.
+
+They are ``ProfileNote`` warnings, printed at most once per key per Python process, so a run's console shows each
+note once. ``warnings.simplefilter("ignore", ProfileNote)`` silences them; ``"error"`` turns them into errors.
+"""
+from __future__ import annotations
+
+import sys
+import warnings
+
+_shown: set[str] = set()
+
+
+class ProfileNote(UserWarning):
+    """A note from a profile card."""
+
+
+def _caller_stacklevel() -> int:
+    """The ``stacklevel`` (for a ``warnings.warn`` in ``emit_note``) of the first frame outside MyoGen."""
+    level, frame = 2, sys._getframe(2)
+    while frame is not None and (str(frame.f_globals.get("__name__")).split(".")[0] == "myogen"
+                                 or frame.f_code.co_filename.startswith("<@beartype")):  # type-check wrappers
+        level, frame = level + 1, frame.f_back
+    return level
+
+
+def emit_note(key: str, text: str) -> None:
+    """Warn ``text`` as a ``ProfileNote`` the first time ``key`` is seen in this process, at the user's line."""
+    if key in _shown:
+        return
+    _shown.add(key)
+    warnings.warn(text, ProfileNote, stacklevel=_caller_stacklevel())
+
+
+def reset_notes() -> None:
+    """Forget which notes were shown (tests)."""
+    _shown.clear()
diff --git a/myogen/profiles/fdi_semg.published.yaml b/myogen/profiles/fdi_semg.published.yaml
new file mode 100644
index 00000000..ce51b16f
--- /dev/null
+++ b/myogen/profiles/fdi_semg.published.yaml
@@ -0,0 +1,27 @@
+# Variant of fdi_semg: each component as its source paper describes it. Fields not listed are inherited.
+name: fdi_semg.published
+parent: fdi_semg
+version: 1
+description: "fdi_semg with the published component models: Fuglevand 1993 force-frequency, Rosenfalck lambda 1, no fibre dispersion, surface (single) fibre CVs. For comparison with those papers and older MyoGen results; the drive is the calibrated one."
+validation:
+  date: "2026-11-07"
+  script: benchmarks/profiles/validate_card.py
+  result: benchmarks/profiles/results/fdi_semg.published.json
+  git_commit: a6fabd87d741a1bbfc7735d18703c0e7e683e07b
+  setup: "as fdi_semg (25 % MVC trapezoid, 8 Delsys 5-pin sensors at 10240 Hz); run and recorded, no limits (published components, not fitted)"
+  cc: {rms_mV: 0.2470, mdf_Hz: 227.5}
+  adj: {rms_mV: 0.3058, mdf_Hz: 220.0}
+  hold_force_pct_mvc: 25.49
+  checked: true
+sections:
+  force:
+    force_frequency_shape: {value: [2.0, 3.0], source: "[fuglevand1993 p. 2475 eq. 16] (S = 1 - exp(-2 r^3))", calibrated: true, note: "This variant's MVC mapping (drive.mvc_model_units, drive.excitation_at_pct_mvc; benchmarks/profiles/results/fdi_semg.published_mvc.json) was measured with this value. The drive itself (trial 85) was fitted with [0.93, 1.29]. Changing it leaves that mapping stale; rerun benchmarks/profiles/validate_card.py --mvc-mapping."}
+  muscle:
+    within_unit_fiber_diameter_log_sd: {value: null, source: "assumption (no within-unit spread; surface CVs)"}
+    unit_cv: {value: null, source: "assumption (surface CVs)"}
+  surface_emg:
+    fiber_dispersion: {value: null, source: "derived: the published model of [farina2004] has no fibre dispersion"}
+    iap_lambda__per_mm: {value: 1.0, source: "assumption (the library default λ of [rosenfalck1969]; secondary)", unit: "1/mm"}
+  drive:
+    mvc_model_units: {value: 13470.808768515415, source: "derived: MVC ramp with the Fuglevand curve (benchmarks/profiles/validate_card.py --mvc-mapping, 2026-11-04; results/fdi_semg.published_mvc.json)"}
+    excitation_at_pct_mvc: {value: {0: 0.0, 5: 0.05927734375, 10: 0.09799804687500001, 20: 0.153955078125, 25: 0.184521484375, 30: 0.212890625, 50: 0.31396484375, 100: 1.0}, source: "derived: MVC ramp with the Fuglevand curve (benchmarks/profiles/validate_card.py --mvc-mapping, 2026-11-04; results/fdi_semg.published_mvc.json); 0 and 100 are the ends of the ramp"}
diff --git a/myogen/profiles/fdi_semg.yaml b/myogen/profiles/fdi_semg.yaml
new file mode 100644
index 00000000..28b73bd8
--- /dev/null
+++ b/myogen/profiles/fdi_semg.yaml
@@ -0,0 +1,62 @@
+# FDI surface EMG profile card. Every value as validated, with its source; see myogen/profiles/_card.py for the format.
+# Sources are docs/literature.md keys with the locator read, 'assumption', or 'derived: ...'.
+name: fdi_semg
+version: 1
+description: "First dorsal interosseous surface EMG with every validated mechanism on: calibrated descending drive with Gfluctdv noise, human force-frequency curve, Muscle fibres typed by unit with within-unit spread and per-unit CVs from recruitment force, surface fibre dispersion, Rosenfalck lambda fitted to Botelho 2019."
+seeds: {thresholds: 42, muscle: 42, surface_emg: 7, drive: 42, cv_scatter: 20261103}
+calibration:
+  file: results/dd_fdi_120mu_fuglevand50_rec41_deluca82_human59_wexp/dd_optimized_params.json (written by examples/02_finetune/06_optimize_dd_multilevel.py; not in the repo, its values are on this card)
+  best_trial: 85
+  objective_value: 1.2470154737228665
+  isi_cv_mean: 0.20307921948055463
+  max_recruitment: 0.41
+  last_recruit_target_pct_mvc: 59.0
+  mvc_ramp: examples/02_finetune/07_fdi_spikes_force.py (2026-11-02)
+  unit_cv_inputs: benchmarks/semg_dispersion/unit_cv/rt_force_cv_scatter.json
+validation:
+  date: "2026-11-07"
+  script: benchmarks/profiles/validate_card.py
+  result: benchmarks/profiles/results/fdi_semg.json
+  git_commit: a6fabd87d741a1bbfc7735d18703c0e7e683e07b
+  setup: "25 % MVC trapezoid, 8 Delsys 5-pin sensors at 10240 Hz; limits frozen 2026-11-03 (cc / adj RMS +-10 % of 0.196 / 0.233 mV, MDF +-5 % of 162.5 / 157.5 Hz)"
+  cc: {rms_mV: 0.2004, mdf_Hz: 161.2}
+  adj: {rms_mV: 0.2340, mdf_Hz: 156.2}
+  hold_force_pct_mvc: 25.62
+  checked: true
+  passed: true
+sections:
+  thresholds:
+    n_units: {value: 120, source: "derived: [feinstein1955] ~119 FDI motor units (secondary; no locator read)", calibrated: true, note: "The FDI drive calibration (trial 85) was fitted with this value. Changing it leaves drive.threshold_hz, drive.excitation_at_pct_mvc and drive.recruitment_pct_mvc stale, so % MVC levels and the per-unit fibre CVs no longer hold; rerun examples/02_finetune/06_optimize_dd_multilevel.py (then 07_fdi_spikes_force.py) to recalibrate."}
+    mode: {value: "fuglevand", source: "assumption (exponential thresholds, [fuglevand1993 p. 2472 eq. 1])", calibrated: true, note: "The FDI drive calibration (trial 85) was fitted with this value. Changing it leaves drive.threshold_hz, drive.excitation_at_pct_mvc and drive.recruitment_pct_mvc stale, so % MVC levels and the per-unit fibre CVs no longer hold; rerun examples/02_finetune/06_optimize_dd_multilevel.py (then 07_fdi_spikes_force.py) to recalibrate."}
+    recruitment_range: {value: 50.0, source: "assumption (50-fold, as calibrated)", calibrated: true, note: "The FDI drive calibration (trial 85) was fitted with this value. Changing it leaves drive.threshold_hz, drive.excitation_at_pct_mvc and drive.recruitment_pct_mvc stale, so % MVC levels and the per-unit fibre CVs no longer hold; rerun examples/02_finetune/06_optimize_dd_multilevel.py (then 07_fdi_spikes_force.py) to recalibrate."}
+  motor_neurons:
+    config_file: {value: "alpha_mn_default.yaml", source: "assumption (library NERLab model, gamma 0.2, as calibrated)", calibrated: true, note: "The FDI drive calibration (trial 85) was fitted with this value. Changing it leaves drive.threshold_hz, drive.excitation_at_pct_mvc and drive.recruitment_pct_mvc stale, so % MVC levels and the per-unit fibre CVs no longer hold; rerun examples/02_finetune/06_optimize_dd_multilevel.py (then 07_fdi_spikes_force.py) to recalibrate."}
+    gfluctdv_scale: {value: 0.358228116077664, source: "derived: fitted by the drive calibration (ISI CV 0.20 at 30 % MVC)", calibrated: true, note: "The FDI drive calibration (trial 85) was fitted with this value. Changing it leaves drive.threshold_hz, drive.excitation_at_pct_mvc and drive.recruitment_pct_mvc stale, so % MVC levels and the per-unit fibre CVs no longer hold; rerun examples/02_finetune/06_optimize_dd_multilevel.py (then 07_fdi_spikes_force.py) to recalibrate."}
+  drive:
+    dd_neurons: {value: 147, source: "derived: fitted by the drive calibration", calibrated: true, note: "The FDI drive calibration (trial 85) was fitted with this value. Changing it leaves drive.threshold_hz, drive.excitation_at_pct_mvc and drive.recruitment_pct_mvc stale, so % MVC levels and the per-unit fibre CVs no longer hold; rerun examples/02_finetune/06_optimize_dd_multilevel.py (then 07_fdi_spikes_force.py) to recalibrate."}
+    conn_probability: {value: 0.3596767670514289, source: "derived: fitted by the drive calibration", calibrated: true, note: "The FDI drive calibration (trial 85) was fitted with this value. Changing it leaves drive.threshold_hz, drive.excitation_at_pct_mvc and drive.recruitment_pct_mvc stale, so % MVC levels and the per-unit fibre CVs no longer hold; rerun examples/02_finetune/06_optimize_dd_multilevel.py (then 07_fdi_spikes_force.py) to recalibrate."}
+    dd_drive__Hz: {value: 759.6041875607232, source: "derived: fitted by the drive calibration", unit: "Hz", calibrated: true, note: "The FDI drive calibration (trial 85) was fitted with this value. Changing it leaves drive.threshold_hz, drive.excitation_at_pct_mvc and drive.recruitment_pct_mvc stale, so % MVC levels and the per-unit fibre CVs no longer hold; rerun examples/02_finetune/06_optimize_dd_multilevel.py (then 07_fdi_spikes_force.py) to recalibrate."}
+    synaptic_weight: {value: 0.05, source: "assumption (fixed in the calibration)", unit: "uS", calibrated: true, note: "The FDI drive calibration (trial 85) was fitted with this value. Changing it leaves drive.threshold_hz, drive.excitation_at_pct_mvc and drive.recruitment_pct_mvc stale, so % MVC levels and the per-unit fibre CVs no longer hold; rerun examples/02_finetune/06_optimize_dd_multilevel.py (then 07_fdi_spikes_force.py) to recalibrate."}
+    gamma_shape: {value: 3.0, source: "assumption (fixed in the calibration)", calibrated: true, note: "The FDI drive calibration (trial 85) was fitted with this value. Changing it leaves drive.threshold_hz, drive.excitation_at_pct_mvc and drive.recruitment_pct_mvc stale, so % MVC levels and the per-unit fibre CVs no longer hold; rerun examples/02_finetune/06_optimize_dd_multilevel.py (then 07_fdi_spikes_force.py) to recalibrate."}
+    network: {value: {"fan_in": "equal", "weight_scaling": "dendritic_transfer", "weight_exponent": 0.6660301395551126}, source: "derived: equal-share wiring [heckmanbinder1991]; exponent fitted by the drive calibration", calibrated: true, note: "The FDI drive calibration (trial 85) was fitted with this value. Changing it leaves drive.threshold_hz, drive.excitation_at_pct_mvc and drive.recruitment_pct_mvc stale, so % MVC levels and the per-unit fibre CVs no longer hold; rerun examples/02_finetune/06_optimize_dd_multilevel.py (then 07_fdi_spikes_force.py) to recalibrate."}
+    threshold_hz: {value: 103.2586942465358, source: "derived: binary search of the drive calibration", unit: "Hz"}
+    timestep_ms: {value: 0.1, source: "assumption (as calibrated)", unit: "ms", calibrated: true, note: "The FDI drive calibration (trial 85) was fitted with this value. Changing it leaves drive.threshold_hz, drive.excitation_at_pct_mvc and drive.recruitment_pct_mvc stale, so % MVC levels and the per-unit fibre CVs no longer hold; rerun examples/02_finetune/06_optimize_dd_multilevel.py (then 07_fdi_spikes_force.py) to recalibrate."}
+    mvc_model_units: {value: 5741.448370053264, source: "derived: MVC ramp (07_fdi_spikes_force.py; mean force over the last 2 s of a 3 s hold at excitation 1)"}
+    excitation_at_pct_mvc: {value: {0: 0.0, 5: 0.048095703125, 10: 0.08164062500000001, 20: 0.15146484375000002, 25: 0.18818359375000002, 30: 0.216552734375, 50: 0.330078125, 100: 1.0}, source: "derived: MVC ramp (fdi_spikes_force.py); 0 and 100 are the ends of the ramp"}
+    recruitment_pct_mvc: {value: [0.46986020323903693, 0.20365455549991282, 0.03166309929804988, 0.40699515423736293, 0.011177170864003425, 1.6758365143738203, 0.15598253843591156, 1.047930805845789, 0.26288766847810846, 1.345581745074569, 1.3110052937003769, 0.004888734189571558, 0.013057205829349263, 0.33019226017154574, 1.233993167425587, 0.04158729017219931, 0.12995102716269244, 0.2984513695425307, 0.8102340294119393, 0.029215060564397922, 0.764533757333056, 0.12995102716269244, 1.6842446807453215, 0.08771088115748524, 0.197991854411853, 0.6508543586566129, 0.8528533378877912, 1.1555454297776637, 0.02222488605192154, 0.13453719555174798, 0.6166207753350973, 3.339957112631774, 1.1888868196062683, 1.3645789157827526, 1.4064255858255854, 0.8941700833142736, 0.058295212664633465, 2.175204188478618, 0.2733587184909703, 0.41804007169088575, 1.1864711701539916, 2.187645980983252, 2.57061952200972, 0.46986020323903693, 2.7479796713329745, 0.7020379878697227, 1.2218919646232764, 1.5227575303173753, 2.0454903197003556, 2.414698833221715, 3.588153768667672, 4.079761278129645, 2.361993518218798, 1.3332428253581197, 1.546268474067718, 1.0240429247315532, 3.2074242245399875, 3.439552094696828, 5.022680375369172, 4.2163736749481675, 3.4748885087089567, 2.5339736138511406, 5.138412274793791, 4.5349004767403915, 1.6933035223080033, 3.7496557892020173, 5.322370906050619, 9.014480309440712, 6.094457923225501, 5.234339352112088, 5.148353677258317, 5.748741605182785, 4.5349004767403915, 5.748741605182785, 7.7442244368661, 7.809433345615285, 8.579071045030739, 8.991553411035962, 11.942718818747275, 10.218601575385467, 8.887154037276314, 9.827794176983302, 11.383546604028536, 9.71543993783791, 9.710584865497216, 11.196236831879483, 10.245796911877068, 12.232767657775055, 17.174502492220903, 14.529522555498156, 14.56262124623736, 17.30600476352565, 16.384794170918614, 21.878415541336654, 19.804187614252506, 17.78648308116625, 19.341434771464694, 24.032869851765838, 28.754174034459318, 23.751445778748074, 25.00601951764433, 27.36797525027772, 29.153911002524218, 32.979968125639445, 31.798633621143995, 31.743242730470456, 32.60342198341561, 32.77561427853402, 43.86945014864069, 40.311785082610015, 42.015527332860096, 45.10120072012499, 50.49296850480229, 49.15915991924408, 44.534854028051306, 51.66439068799985, 52.095461304932044, 57.48199613445284, 62.12949871797975, 65.88336085391157], source: "derived: MVC ramp; first of >= 3 spikes within 500 ms"}
+  force:
+    force_frequency_shape: {value: [0.93, 1.29], source: "[macefield1996 Fig. 6C], [fuglevand1999 Table 1, Fig. 6A], [mcnulty2000 Table 2]", calibrated: true, note: "The FDI drive calibration (trial 85) was fitted with this value. Changing it leaves drive.threshold_hz, drive.excitation_at_pct_mvc and drive.recruitment_pct_mvc stale, so % MVC levels and the per-unit fibre CVs no longer hold; rerun examples/02_finetune/06_optimize_dd_multilevel.py (then 07_fdi_spikes_force.py) to recalibrate."}
+    recording_frequency__Hz: {value: 2048.0, source: "assumption (as calibrated)", unit: "Hz"}
+  muscle:
+    radius__mm: {value: 6.91, source: "assumption (class default, attributed to [jacobson1992]; Table 1 unverified)", unit: "mm"}
+    fiber_density__fibers_per_mm2: {value: 267.0, source: "derived: [feinstein1955] ~40,000 FDI fibres (secondary; no locator read) over pi * 6.91^2 = 150 mm2", unit: "mm**-2", note: "Amplitude was checked against Botelho 2019 (RMS 0.19-0.26 mV at 25 % MVC) with this fibre density; amplitude scales with the number of fibres per unit."}
+    skin_thickness__mm: {value: 1.2, source: "assumption (median hand-dorsum skin of the healthy controls of [li2018]; locator not recorded)", unit: "mm", note: "Amplitude was checked against Botelho 2019 (RMS 0.19-0.26 mV at 25 % MVC) with these layers (skin 1.2 mm, fat 1.3 mm, skin 1 S/m); thicker layers or a resistive skin change it several-fold."}
+    fat_thickness__mm: {value: 1.3, source: "derived: [nishikawa2024] subcutaneous depth minus the skin", unit: "mm", note: "Amplitude was checked against Botelho 2019 (RMS 0.19-0.26 mV at 25 % MVC) with these layers (skin 1.2 mm, fat 1.3 mm, skin 1 S/m); thicker layers or a resistive skin change it several-fold."}
+    skin_conductivity__S_per_m: {value: 1.0, source: "[roeleveld1997 Fig. 4] model III effective value (also the value of Farina et al. 2004)", unit: "S/m", note: "Amplitude was checked against Botelho 2019 (RMS 0.19-0.26 mV at 25 % MVC) with these layers (skin 1.2 mm, fat 1.3 mm, skin 1 S/m); thicker layers or a resistive skin change it several-fold."}
+    fiber_model: {value: "typed", source: "assumption (one fibre type per motor unit in recruitment order, [polgar1973 p. 312 Table 4])", note: "lambda 0.63 (surface_emg.iap_lambda__per_mm) was fitted with this setting; changing it moves the spectrum and amplitude, so the fit to Botelho 2019 no longer holds."}
+    within_unit_fiber_diameter_log_sd: {value: [0.08, 0.13], source: "[martin1988 p. C46 Table 1]", note: "lambda 0.63 (surface_emg.iap_lambda__per_mm) was fitted with this setting; changing it moves the spectrum and amplitude, so the fit to Botelho 2019 no longer holds."}
+    unit_cv: {value: {"from": "recruitment_force", "residual_sd__m_s": 0.39}, source: "[delvecchio2018 p. 2]; residual derived from R2 0.70", note: "lambda 0.63 (surface_emg.iap_lambda__per_mm) was fitted with this setting; changing it moves the spectrum and amplitude, so the fit to Botelho 2019 no longer holds."}
+  surface_emg:
+    fiber_dispersion: {value: {"nmj_model": "two_layer_intra", "nmj_scatter_scale": 1.0, "nerve_delays": true, "nerve_cvs__m_per_s": [5.0, 2.0], "fiber_cv_threshold_slope": 0.0, "fiber_cv_within_mu_cov": 0.0, "nmj_coefficients": "published", "tendon_end_band__mm": 5.0, "tendon_end_correlation__mm": 5.0, "tendon_end_jitter__mm": 0.3, "diameter_scales_source": true, "fiber_properties": "muscle"}, source: "[akhmadeev2020 eq. 5-6] end-plates; other fields as documented in SurfaceFiberDispersion", note: "lambda 0.63 (surface_emg.iap_lambda__per_mm) was fitted with this setting; changing it moves the spectrum and amplitude, so the fit to Botelho 2019 no longer holds."}
+    iap_lambda__per_mm: {value: 0.63, source: "derived: fit to [botelho2019] (PR doc refit 2026-11-03)", unit: "1/mm", note: "lambda 0.63 was fitted to Botelho 2019 (FDI, 25 % MVC) with Muscle fibres, per-unit CVs from recruitment force and the fibre dispersion on this card. With surface CVs the fit is 0.7; at 1.0 the median frequency is ~35 Hz higher. Below 1 the IAP is longer, so keep iap_kernel_length__mm long enough (the window check warns)."}
+    iap_kernel_length__mm: {value: 200.0, source: "assumption (long enough for lambda 0.63)", unit: "mm"}
diff --git a/myogen/simulator/core/emg/surface/surface_emg.py b/myogen/simulator/core/emg/surface/surface_emg.py
index 266251f1..ec4fd6e5 100644
--- a/myogen/simulator/core/emg/surface/surface_emg.py
+++ b/myogen/simulator/core/emg/surface/surface_emg.py
@@ -319,6 +319,14 @@ def simulate_muaps(
         The motor units are processed in parallel using joblib, with each motor unit's
         fibers processed sequentially to maintain optimization efficiency.
         """
+        if getattr(self, "profile_rng_state", None) is not None:  # built from a profile card: draw as at build time,
+            from myogen import _seed_scope                          # then give the caller's generator back
+
+            with _seed_scope(self.profile_seed, state=self.profile_rng_state, reset_cell_ids=False):
+                return self._simulate_muaps(n_jobs, verbose, use_gpu, execution)
+        return self._simulate_muaps(n_jobs, verbose, use_gpu, execution)
+
+    def _simulate_muaps(self, n_jobs, verbose, use_gpu, execution) -> SURFACE_MUAP__Block:
         if execution is not None:
             from ..runtime.semg import simulate_muaps
             return simulate_muaps(self, execution, n_jobs, verbose, use_gpu)
@@ -981,6 +989,8 @@ def simulate_surface_emg(self, spike_train__Block: SPIKE_TRAIN__Block, verbose:
                 )
 
         # Store results privately
+        if getattr(self, "profile_record", None) is not None:  # built from a profile card
+            block.annotations["myogen_profile"] = self.profile_record
         self._surface_emg__Block = block
         return block
 
diff --git a/myogen/simulator/neuron/descending_drive_network.py b/myogen/simulator/neuron/descending_drive_network.py
new file mode 100644
index 00000000..7ce3e8d5
--- /dev/null
+++ b/myogen/simulator/neuron/descending_drive_network.py
@@ -0,0 +1,163 @@
+"""Calibrated descending-drive network: the DD pool, its wiring to the motor neurons, and the drive loop.
+
+This is the drive the FDI calibration (``examples/02_finetune/06_optimize_dd_multilevel.py``) was fitted with and the
+``fdi_semg`` profile runs: a gamma-process descending-drive pool, equal fan-in onto every motor neuron with
+dendritic-transfer weights, Gfluctdv conductance noise at a fitted scale, and a drive rate in Hz stepped with the
+NEURON integration. The optimiser imports these functions, so the calibration and the profile share one authority.
+"""
+from types import SimpleNamespace
+
+import numpy as np
+import quantities as pq
+from neo import Block, Segment, SpikeTrain
+from neuron import h
+
+from myogen.simulator.neuron.network import Network
+from myogen.simulator.neuron.populations import AlphaMN__Pool, DescendingDrive__Pool
+from myogen.utils.helper import get_gamma_shape_for_mvc
+
+NETWORK = {"fan_in": "equal", "weight_scaling": "dendritic_transfer"}  # DD -> MN wiring, see build_dd_mn_network
+SPIKE_SAMPLING_RATE__HZ = 10240.0  # sampling rate attached to the recorded spike trains
+
+
+def insert_scaled_gfluctdv(mn_pool, scale):
+    """Insert Gfluctdv independent noise on all MN dendrites at given scale."""
+    if scale <= 0:
+        return
+    for cell in mn_pool:
+        cell.insert_Gfluctdv()
+        for d in cell.dend:
+            d.std_e_Gfluctdv *= scale
+            d.std_i_Gfluctdv *= scale
+            d.g_e0_Gfluctdv *= scale
+            d.g_i0_Gfluctdv *= scale
+
+
+def dendritic_transfer(mn_pool):
+    """Fraction of a steady current injected at each MN's DD synapse site that reaches its clamped soma.
+
+    DC transfer impedance (synapse site -> soma) over soma input impedance, from ``h.Impedance`` linearised at
+    the initial state (active channels included). Call before ``insert_scaled_gfluctdv``.
+    """
+    for s, v in zip(*mn_pool.get_initialization_data()):
+        s.v = v
+    h.finitialize()
+    k = []
+    for cell in mn_pool:
+        syn = next(s for s in cell.synapse__list if s.e >= -40)
+        seg = syn.get_segment()
+        imp = h.Impedance()
+        imp.loc(0.5, sec=cell.soma)
+        imp.compute(0, 1)
+        k.append(imp.transfer(seg.x, sec=seg.sec) / imp.input(0.5, sec=cell.soma))
+    return np.array(k)
+
+
+def build_dd_mn_network(mn_pool, dd_pool, dd_params):
+    """Wire the DD pool to the motor neurons and the cortical input to the DD pool.
+
+    ``dd_params["network"]`` must be ``NETWORK``; parameter files without it (calibrated with random fan-in and
+    one weight for all neurons) are refused rather than reinterpreted.
+
+    - Equal fan-in: every MN receives exactly round(conn_probability * dd_neurons) DD inputs, drawn without
+      replacement. Modelling assumption (descending drive shared across the pool); it removes the binomial
+      +/-13 % input-count scatter that reordered the largest units, whose thresholds are ~3 % apart.
+    - Dendritic-transfer weights: weight_i = synaptic_weight * (k_max / k_i) ** p (``dendritic_transfer``), p from
+      ``network["weight_exponent"]`` (default 1). With p = 1 each input delivers the same current to a clamped soma
+      at low drive whatever the neuron's size, and recruitment order follows rheobase: equal effective synaptic
+      current per input, the equal-share case of Heckman & Binder 1991 (J Neurophysiol 65:952) [heckmanbinder1991]. p < 1 gives larger
+      neurons less effective current per input (p = 0: one weight for all). Modelling assumption either way;
+      measured effective currents are input-system specific (cat data) and none are known for corticospinal
+      input to a hand muscle. Not compensated: the Gfluctdv mean-conductance shunt and the driving-force
+      saturation at high drive.
+    """
+    net_cfg = dict(dd_params.get("network") or {})
+    p = net_cfg.pop("weight_exponent", 1.0)
+    if not isinstance(p, (int, float)) or isinstance(p, bool) or not np.isfinite(p) or p < 0:
+        raise ValueError(f"network weight_exponent must be a finite number >= 0, got {p!r}")
+    if net_cfg != NETWORK:
+        raise ValueError(
+            f"dd_params['network'] must be {NETWORK}, got {net_cfg!r}; parameter files from MyoGen 0.13.0 or "
+            "earlier were calibrated with random fan-in and uniform weights and must be recalibrated."
+        )
+    n_dd = len(dd_pool)
+    n_in = int(round(dd_params["conn_probability"] * n_dd))
+    if n_in < 1:
+        raise ValueError(f"conn_probability * dd_neurons rounds to {n_in} inputs per motor neuron")
+    transfer = dendritic_transfer(mn_pool)
+    weights = dd_params["synaptic_weight"] * (transfer.max() / transfer) ** p
+    network = Network({"DD": dd_pool, "aMN": mn_pool})
+    dd_to_mn = network.connect_fan_in("DD", "aMN", n_inputs=n_in, weights__uS=weights * pq.uS,
+                                      delay__ms=1.0 * pq.ms)
+    network.connect_from_external(source="cortical_input", target="DD", weight__uS=1.0 * pq.uS)
+    return SimpleNamespace(network=network, dd_to_mn=dd_to_mn,
+                           cortical_netcons=network.get_netcons("cortical_input", "DD"),
+                           transfer=transfer, weights=weights)
+
+
+def excitation_to_drive_hz(segments, threshold_hz, max_hz, timestep_ms=0.1, noise_seed=42):
+    """Piecewise-linear excitation [(duration_s, start, end)] and its DD rate (Hz): threshold_hz + e (max - threshold),
+    0.8 threshold_hz at rest, plus 1 % half-normal noise."""
+    parts = [np.linspace(a, b, int(round(d * 1000 / timestep_ms)), endpoint=False) for d, a, b in segments]
+    e = np.concatenate(parts)
+    hz = np.where(e > 0, threshold_hz + e * (max_hz - threshold_hz), 0.8 * threshold_hz)
+    hz = hz + np.abs(np.random.default_rng(noise_seed).normal(0, hz.max() * 0.01, size=hz.size))
+    return e, hz
+
+
+def build_pools(thresholds, dd_params, timestep_ms=0.1, config_file="alpha_mn_default.yaml"):
+    """Motor neuron pool, DD pool and their wiring (``build_dd_mn_network``), with Gfluctdv at the fitted scale."""
+    gamma_shape = get_gamma_shape_for_mvc(100, mvc_shape_value=dd_params["gamma_shape"])
+    mn_pool = AlphaMN__Pool(recruitment_thresholds__array=thresholds, config_file=config_file)
+    dd_pool = DescendingDrive__Pool(n=dd_params["dd_neurons"], timestep__ms=timestep_ms * pq.ms,
+                                    process_type="gamma", shape=gamma_shape)
+    wiring = build_dd_mn_network(mn_pool, dd_pool, dd_params)
+    insert_scaled_gfluctdv(mn_pool, dd_params.get("gfluctdv_scale", 0.0))
+    return SimpleNamespace(mn_pool=mn_pool, dd_pool=dd_pool, wiring=wiring)
+
+
+def run_descending_drive(thresholds, drive_hz, dd_params, total_ms, timestep_ms=0.1,
+                         config_file="alpha_mn_default.yaml"):
+    """Run the DD network with a time-varying drive rate (Hz, one value per time step); return the MN spike trains."""
+    pools = build_pools(thresholds, dd_params, timestep_ms, config_file)
+    mn_pool, dd_pool, dd_netcons = pools.mn_pool, pools.dd_pool, pools.wiring.cortical_netcons
+
+    mn_recorders = []
+    for cell in mn_pool:
+        rec = h.Vector()
+        nc = h.NetCon(cell.soma(0.5)._ref_v, None, sec=cell.soma)
+        nc.threshold = 50
+        nc.record(rec)
+        mn_recorders.append(rec)
+
+    h.load_file("stdrun.hoc")
+    h.dt = timestep_ms
+    h.tstop = total_ms
+
+    for s, v in zip(*mn_pool.get_initialization_data()):
+        s.v = v
+    for s, v in zip(*dd_pool.get_initialization_data()):
+        s.v = v
+    h.finitialize()
+
+    n_steps = len(drive_hz)
+    step = 0
+    while h.t < h.tstop:
+        d = drive_hz[min(step, n_steps - 1)]
+        for dd_cell in dd_pool:
+            if dd_cell.integrate(d):
+                if h.t < h.tstop:
+                    dd_netcons[dd_cell.pool__ID].event(h.t + 1)
+        h.fadvance()
+        step += 1
+
+    fs = SPIKE_SAMPLING_RATE__HZ
+    block = Block()
+    seg = Segment()
+    block.segments.append(seg)
+    seg.spiketrains = [
+        SpikeTrain(rec.as_numpy() / 1000 * pq.s, t_stop=total_ms / 1000.0 * pq.s, sampling_rate=fs * pq.Hz,
+                   sampling_period=(1 / fs) * pq.s, name=str(i))
+        for i, rec in enumerate(mn_recorders)
+    ]
+    return block
diff --git a/pyproject.toml b/pyproject.toml
index f72a4dbb..15214157 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -42,6 +42,7 @@ dependencies = [
     "numpy>=1.26",
     "optuna>=4.5.0",
     "pandas>=2.0",
+    "pyyaml>=6.0",
     "quantities>=0.16.2",
     "scienceplots>=2.1.1",
     "scikit-fmm>=2025.1.29",
diff --git a/tests/test_descending_drive_network.py b/tests/test_descending_drive_network.py
new file mode 100644
index 00000000..ca57a413
--- /dev/null
+++ b/tests/test_descending_drive_network.py
@@ -0,0 +1,86 @@
+"""The calibrated descending-drive network (moved from the optimiser, examples/02_finetune/06_optimize_dd_multilevel.py).
+
+The stored spikes were computed with the optimiser's own functions before the move (2026-11-03), so the library
+module reproduces the calibrated drive exactly. Gfluctdv is off in the pinned run: its ``normrand`` stream is
+NEURON's and not seeded from Python.
+"""
+import numpy as np
+import pytest
+
+from myogen import load_nmodl_mechanisms, set_random_seed
+from myogen.simulator import RecruitmentThresholds
+
+# FDI calibration _human59_wexp, best trial 85 (copied into the fdi_semg card)
+DD_PARAMS = {
+    "dd_neurons": 147,
+    "conn_probability": 0.3596767670514289,
+    "dd_drive__Hz": 759.6041875607232,
+    "synaptic_weight": 0.05,
+    "gamma_shape": 3.0,
+    "network": {"fan_in": "equal", "weight_scaling": "dendritic_transfer", "weight_exponent": 0.6660301395551126},
+    "gfluctdv_scale": 0.0,
+}
+THRESHOLD_HZ = 103.2586942465358
+# Stored spikes: count per motor unit, the first spike time (s) of six units, and the sum of all spike times (s).
+SPIKES_PER_UNIT = [
+    22, 23, 21, 23, 22, 23, 21, 22, 21, 23, 22, 22, 21, 21, 23, 22, 22, 22, 21, 22, 21, 21, 21, 21, 22,
+    22, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 20, 21, 21, 21, 21, 20, 20,
+    20, 20, 20, 20, 21, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 19, 20, 19, 19, 19, 18, 18, 19, 18, 18,
+    18, 18, 18, 18, 18, 17, 17, 17, 16, 16, 16, 16, 15, 15, 15, 14, 14, 13, 13, 12, 12, 12, 11, 11, 10,
+    10, 9, 8, 8, 7, 6, 5, 4, 5, 3, 3, 2, 1, 1, 0, 0, 0, 0, 0, 0
+]
+FIRST_SPIKE__S = {0: 0.18880000000010866, 19: 0.2636000000001256, 38: 0.5264000000001691, 57: 0.546200000000151, 76: 0.6731000000000356, 95: 0.9732999999997626}
+SPIKE_TIME_SUM__S = 2193.9812999992955
+
+
+def _thresholds():
+    set_random_seed(42)
+    thr, _ = RecruitmentThresholds(N=120, recruitment_range__ratio=50, mode="fuglevand")
+    thr = np.asarray(thr, dtype=float)
+    return thr / thr.max()
+
+
+def test_run_matches_stored_spikes():
+    from myogen.simulator.neuron.descending_drive_network import excitation_to_drive_hz, run_descending_drive
+
+    load_nmodl_mechanisms()
+    thr = _thresholds()
+    e, hz = excitation_to_drive_hz([(0.5, 0, 0), (1.0, 0, 0.3)], THRESHOLD_HZ, DD_PARAMS["dd_drive__Hz"])
+    block = run_descending_drive(thr, hz, DD_PARAMS, len(e) * 0.1)
+    trains = block.segments[0].spiketrains
+    # NEURON's floating point differs between platforms. Once the network diverges, spike times shift by whole time
+    # steps, so a unit can gain or lose a spike near the end of the run. Against the stored run (2000 spikes): macOS
+    # arm64 matches every count and its sum is 1.1e-4 higher (2194.218 s); Linux x86_64 fires once more in unit 8
+    # (at 1.496 s of 1.5 s) and its sum is 5.9e-4 higher (2195.276 s). First spikes match on both to 1e-9.
+    counts = np.array([len(st) for st in trains])
+    assert np.abs(counts - SPIKES_PER_UNIT).max() <= 1
+    for unit, t0 in FIRST_SPIKE__S.items():
+        assert float(trains[unit].rescale("s").magnitude[0]) == pytest.approx(t0, rel=1e-9)
+    total = sum(float(np.sum(st.rescale("s").magnitude)) for st in trains)
+    assert total == pytest.approx(SPIKE_TIME_SUM__S, rel=1e-3)
+
+
+def test_gfluctdv_scale_changes_noise():
+    from myogen.simulator.neuron.descending_drive_network import insert_scaled_gfluctdv
+    from myogen.simulator.neuron.populations import AlphaMN__Pool
+
+    load_nmodl_mechanisms()
+    full = AlphaMN__Pool(recruitment_thresholds__array=np.array([0.5]))
+    half = AlphaMN__Pool(recruitment_thresholds__array=np.array([0.5]))
+    insert_scaled_gfluctdv(full, 1.0)
+    insert_scaled_gfluctdv(half, 0.5)
+    assert half[0].dend[0].std_e_Gfluctdv == pytest.approx(0.5 * full[0].dend[0].std_e_Gfluctdv)
+    assert half[0].dend[0].g_i0_Gfluctdv == pytest.approx(0.5 * full[0].dend[0].g_i0_Gfluctdv)
+
+
+def test_network_block_is_required():
+    from myogen.simulator.neuron.descending_drive_network import build_dd_mn_network
+    from myogen.simulator.neuron.populations import AlphaMN__Pool, DescendingDrive__Pool
+    import quantities as pq
+
+    load_nmodl_mechanisms()
+    mn = AlphaMN__Pool(recruitment_thresholds__array=np.array([0.5, 1.0]))
+    dd = DescendingDrive__Pool(n=4, timestep__ms=0.1 * pq.ms, process_type="gamma", shape=3.0)
+    params = {k: v for k, v in DD_PARAMS.items() if k != "network"}
+    with pytest.raises(ValueError, match="network"):
+        build_dd_mn_network(mn, dd, params)
diff --git a/tests/test_profiles_builders.py b/tests/test_profiles_builders.py
new file mode 100644
index 00000000..18ba7531
--- /dev/null
+++ b/tests/test_profiles_builders.py
@@ -0,0 +1,173 @@
+"""Card builders: card values reach the library objects, overrides are noted and recorded, builds are reseeded."""
+import json
+import warnings
+from pathlib import Path
+
+import numpy as np
+import pytest
+import quantities as pq
+
+from myogen import load_profile
+from myogen.profiles import ProfileNote
+from myogen.profiles._notes import reset_notes
+
+ROOT = Path(__file__).resolve().parents[1]
+# The validated arm-F muscle (benchmarks/semg_dispersion arm delsys5r_wexp_musclertsc_lam0.63): fibres per unit,
+# and fibre diameter (mm) and CV (mm/s) at evenly spaced fibre indices, with their sums.
+ARM_F_FIBRES_PER_UNIT = [
+    22, 21, 27, 29, 30, 26, 27, 28, 28, 29, 30, 32, 28, 33, 31, 43, 49, 43, 46, 42, 47, 45, 57, 52, 43,
+    54, 69, 64, 67, 60, 54, 72, 71, 86, 79, 79, 83, 93, 73, 78, 91, 90, 101, 100, 122, 108, 112, 100,
+    121, 105, 112, 119, 125, 140, 143, 170, 186, 167, 174, 216, 186, 227, 203, 203, 230, 243, 219, 265,
+    262, 265, 302, 284, 269, 282, 317, 331, 307, 330, 347, 340, 351, 364, 393, 370, 422, 399, 422, 418,
+    447, 442, 483, 539, 441, 537, 646, 641, 650, 605, 686, 713, 716, 712, 777, 739, 834, 850, 868, 822,
+    867, 1058, 1074, 1028, 1087, 1142, 1162, 1264, 1171, 1334, 1319, 1274
+]
+ARM_F_INDICES = [0, 5721, 11442, 17164, 22885, 28607, 34328, 40050]
+ARM_F_DIAMETERS__MM = [
+    0.06728408646276396, 0.06839268111349746, 0.047277283638866574, 0.06324227369016999,
+    0.06061115896061368, 0.054665281872642656, 0.0429479657676721, 0.04763986418537935
+]
+ARM_F_CVS__MM_PER_S = [
+    4809.070973006686, 4864.50070554336, 3806.142092441618, 4758.021558504406, 4626.46582202659,
+    4213.794291232723, 3634.5395187905533, 3946.1445527315545
+]
+ARM_F_SUMS = (2294.4213299653193, 175695618.33762425)  # diameters (mm), CVs (mm/s)
+
+
+def _notes(rec):
+    return [str(w.message) for w in rec if issubclass(w.category, ProfileNote)]
+
+
+def test_muscle_kwargs_match_card():
+    card = load_profile()
+    kw = card.muscle_kwargs()
+    assert kw["radius__mm"].rescale(pq.mm).magnitude == pytest.approx(6.91)
+    assert kw["fiber_density__fibers_per_mm2"].rescale(pq.mm**-2).magnitude == pytest.approx(267.0)
+    assert kw["skin_thickness__mm"].rescale(pq.mm).magnitude == pytest.approx(1.2)
+    assert kw["fat_thickness__mm"].rescale(pq.mm).magnitude == pytest.approx(1.3)
+    assert kw["skin_conductivity__S_per_m"].rescale(pq.S / pq.m).magnitude == pytest.approx(1.0)
+    assert tuple(kw["within_unit_fiber_diameter_log_sd"]) == (0.08, 0.13)
+    cv = np.asarray(kw["motor_unit_conduction_velocities__m_per_s"])
+    assert cv.shape == (120,)
+    np.testing.assert_array_equal(cv, np.asarray(card.muscle_kwargs()["motor_unit_conduction_velocities__m_per_s"]))
+    assert len(kw["recruitment_thresholds"]) == 120
+
+
+def test_unit_cvs_are_the_validated_ones():
+    validated = json.load(open(ROOT / "benchmarks" / "semg_dispersion" / "unit_cv" / "rt_force_cv_scatter.json"))["cv__m_s"]
+    cv = np.asarray(load_profile().muscle_kwargs()["motor_unit_conduction_velocities__m_per_s"])
+    np.testing.assert_array_equal(cv, validated)
+
+
+def test_override_equal_to_card_is_not_an_override():
+    reset_notes()
+    card = load_profile()
+    with warnings.catch_warnings(record=True) as rec:
+        warnings.simplefilter("always")
+        card.muscle_kwargs(skin_thickness__mm=0.0012 * pq.m)
+    assert _notes(rec) == []
+    assert card.overrides() == {}
+
+
+def test_override_prints_note_and_is_recorded():
+    reset_notes()
+    card = load_profile()
+    with warnings.catch_warnings(record=True) as rec:
+        warnings.simplefilter("always")
+        kw = card.surface_emg_kwargs(iap_lambda__per_mm=1.0)
+        card.surface_emg_kwargs(iap_lambda__per_mm=1.0)
+    notes = _notes(rec)
+    assert len(notes) == 1 and "0.63" in notes[0]
+    assert kw["iap_lambda__per_mm"] == 1.0
+    assert card.overrides() == {"surface_emg.iap_lambda__per_mm": 1.0}
+    assert card.record() == {"profile": "fdi_semg", "version": 1, "parent": None, "parent_version": None,
+                             "overrides": {"surface_emg.iap_lambda__per_mm": 1.0}}
+
+
+def test_note_names_other_overrides():
+    reset_notes()
+    card = load_profile()
+    with warnings.catch_warnings(record=True) as rec:
+        warnings.simplefilter("always")
+        card.muscle_kwargs(skin_thickness__mm=2.0 * pq.mm)
+        card.surface_emg_kwargs(iap_lambda__per_mm=1.0)
+    assert "muscle.skin_thickness__mm" in _notes(rec)[-1]
+
+
+def test_calibrated_override_warns_calibration_stale():
+    reset_notes()
+    card = load_profile()
+    with warnings.catch_warnings(record=True) as rec:
+        warnings.simplefilter("always")
+        card.thresholds(recruitment_range=30.0)
+    assert any("06_optimize_dd_multilevel.py" in n for n in _notes(rec))
+    assert card.overrides() == {"thresholds.recruitment_range": 30.0}
+
+
+def test_unknown_override_raises():
+    with pytest.raises(TypeError, match="skin_thicknes__mm"):
+        load_profile().muscle_kwargs(skin_thicknes__mm=1.0)
+
+
+def test_force_model_uses_card_shape():
+    card = load_profile()
+    fm = card.force_model()
+    assert fm.force_frequency_shape == (0.93, 1.29)
+    assert fm.profile_record == card.record({})  # no pass-through keywords
+    assert load_profile("fdi_semg.published").force_model().force_frequency_shape == (2.0, 3.0)
+
+
+def test_builds_are_reseeded_and_match_the_validated_fibres():
+    """Slow (full 120-unit FDI muscle, twice)."""
+    card = load_profile()
+
+    def fibres():
+        m = card.muscle(n_jobs=4)
+        d = np.asarray(m.muscle_fiber_diameters__mm.magnitude, dtype=float)
+        cv = np.asarray(m.muscle_fiber_conduction_velocities__mm_per_s.magnitude, dtype=float)
+        counts = [len(m.resulting_fiber_assignment(i)) for i in range(120)]
+        return counts, d, cv
+
+    counts, d, cv = fibres()
+    np.random.default_rng(0).normal(size=10)
+    from myogen import set_random_seed
+    set_random_seed(999)  # another builder or user code moved the global generator
+    counts2, d2, cv2 = fibres()
+    assert counts2 == counts
+    np.testing.assert_array_equal(d2, d)
+    np.testing.assert_array_equal(cv2, cv)
+    assert counts == ARM_F_FIBRES_PER_UNIT
+    np.testing.assert_allclose(d[ARM_F_INDICES], ARM_F_DIAMETERS__MM, rtol=1e-9)
+    np.testing.assert_allclose(cv[ARM_F_INDICES], ARM_F_CVS__MM_PER_S, rtol=1e-9)
+    np.testing.assert_allclose((d.sum(), cv.sum()), ARM_F_SUMS, rtol=1e-9)
+
+
+def test_fibre_properties_only_with_muscle_fibres(monkeypatch):
+    from myogen.simulator import Muscle
+    calls = []
+    for step in ("distribute_innervation_centers", "generate_muscle_fiber_centers", "assign_mfs2mns"):
+        monkeypatch.setattr(Muscle, step, lambda self, *a, **k: None)
+    monkeypatch.setattr(Muscle, "_generate_fiber_properties", lambda self: calls.append(1))
+    load_profile().muscle()
+    assert calls == [1]
+    load_profile("fdi_semg.published").muscle()
+    assert calls == [1]  # no dispersion on the variant: surface CVs, no fibre properties
+
+
+def test_sampling_frequency_is_the_users_choice():
+    card = load_profile()
+    kw = card.surface_emg_kwargs()
+    assert "sampling_frequency__Hz" not in kw and "internal_sampling_frequency__Hz" not in kw
+    kw = card.surface_emg_kwargs(sampling_frequency__Hz=10240.0 * pq.Hz)
+    assert kw["sampling_frequency__Hz"] == 10240.0 * pq.Hz and "internal_sampling_frequency__Hz" not in kw
+    assert card.overrides() == {}  # not a card field: passed through, not recorded
+
+
+@pytest.mark.parametrize("name", ["fdi_semg", "fdi_semg.published"])
+def test_cards_use_the_typed_fibre_model(name):
+    """The cards were validated with fibre types by motor unit, which Muscle makes opt-in."""
+    from myogen.simulator import Muscle
+
+    kw = load_profile(name).muscle_kwargs()
+    assert kw["fiber_model"] == "typed"
+    Muscle(**kw, autorun=False)  # typed-only arguments are accepted
diff --git a/tests/test_profiles_card.py b/tests/test_profiles_card.py
new file mode 100644
index 00000000..a8976498
--- /dev/null
+++ b/tests/test_profiles_card.py
@@ -0,0 +1,93 @@
+"""Profile cards: loading, variants and sources."""
+import re
+from pathlib import Path
+
+import pytest
+
+from myogen import load_profile
+from myogen.profiles import ProfileError
+
+REPO = Path(__file__).resolve().parents[1]
+REGISTER = (REPO / "docs" / "literature.md").read_text(encoding="utf-8")
+SHIPPED = ("fdi_semg", "fdi_semg.published")
+
+
+def test_default_is_fdi_semg():
+    card = load_profile()
+    assert card.name == "fdi_semg"
+    assert card.version == 1
+    assert card.parent is None
+
+
+def test_card_values():
+    card = load_profile()
+    assert card.value("muscle.fiber_density__fibers_per_mm2") == 267
+    assert card.value("surface_emg.iap_lambda__per_mm") == 0.63
+    assert card.value("thresholds.n_units") == 120
+    assert len(card.value("drive.recruitment_pct_mvc")) == 120
+    assert card.value("drive.excitation_at_pct_mvc")[25] == pytest.approx(0.18818359375000002)
+    assert card.field("motor_neurons.gfluctdv_scale").calibrated is True
+    assert card.field("muscle.skin_thickness__mm").unit == "mm"
+
+
+@pytest.mark.parametrize("name", SHIPPED)
+def test_every_field_has_a_registered_source(name):
+    for path, field in load_profile(name).fields().items():
+        src = field.source.strip()
+        assert src, f"{path}: empty source"
+        if src.startswith("assumption"):
+            continue
+        keys = re.findall(r"\[([a-z]+[0-9]{4}[a-z]?)\b", src)
+        assert keys or src.startswith("derived"), f"{path}: source {src!r} cites no register key"
+        for key in keys:
+            assert f'id="{key}"' in REGISTER, f"{path}: {key} is not in docs/literature.md"
+
+
+def test_variant_inherits_and_changes():
+    card = load_profile("fdi_semg.published")
+    assert card.parent == "fdi_semg"
+    assert card.value("surface_emg.iap_lambda__per_mm") == 1.0
+    assert card.value("force.force_frequency_shape") == [2.0, 3.0]
+    assert card.value("surface_emg.fiber_dispersion") is None
+    assert card.value("muscle.fiber_density__fibers_per_mm2") == 267
+
+
+def _write(tmp_path, name, text):
+    p = tmp_path / f"{name}.yaml"
+    p.write_text(text, encoding="utf-8")
+    return p
+
+
+def test_variant_unknown_field_raises(tmp_path):
+    p = _write(tmp_path, "bad", "name: bad\nparent: fdi_semg\nversion: 1\nsections:\n"
+                                "  muscle:\n    skin_thicknes__mm: {value: 2.0, source: assumption}\n")
+    with pytest.raises(ProfileError, match="muscle.skin_thicknes__mm"):
+        load_profile(p)
+
+
+def test_missing_required_field_raises(tmp_path):
+    p = _write(tmp_path, "bad", "name: bad\nversion: 1\nseeds: {}\nsections:\n"
+                                "  muscle:\n    radius__mm: {value: 6.91}\n")
+    with pytest.raises(ProfileError, match="muscle.radius__mm"):
+        load_profile(p)
+    p = _write(tmp_path, "bad2", "name: bad2\nversion: 1\nseeds: {}\n")
+    with pytest.raises(ProfileError, match="sections"):
+        load_profile(p)
+
+
+def test_unknown_card_name_raises():
+    with pytest.raises(ProfileError, match="fdi_semg"):
+        load_profile("vl_semg")
+
+
+@pytest.mark.parametrize("name", SHIPPED)
+def test_validated_cards_record_their_validation(name):
+    card = load_profile(name)
+    assert card.validation is not None, f"{name} has no validation block"
+    v = card.validation
+    for key in ("date", "script", "result", "cc", "adj", "checked"):
+        assert key in v, f"{name}: validation.{key} missing"
+    assert v["checked"] is True
+    assert (REPO / v["result"]).is_file()
+    if name == "fdi_semg":
+        assert v["passed"] is True
diff --git a/tests/test_profiles_delsys_scorer.py b/tests/test_profiles_delsys_scorer.py
new file mode 100644
index 00000000..4ff5f6d1
--- /dev/null
+++ b/tests/test_profiles_delsys_scorer.py
@@ -0,0 +1,32 @@
+"""The Delsys cc/adj scorer reproduces the archived arm-F numbers the fdi_semg limits were frozen from."""
+import importlib.util
+import os
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+ROOT = Path(__file__).resolve().parents[1]
+# Arm delsys5r_wexp_musclertsc_lam0.63 of benchmarks/semg_dispersion (33 MB, not in the repo); set MYOGEN_ARM_F_NPZ to its path.
+ARM_F = Path(os.environ.get("MYOGEN_ARM_F_NPZ", "arm_F_not_set.npz"))
+
+
+def _delsys():
+    spec = importlib.util.spec_from_file_location("delsys", ROOT / "benchmarks" / "profiles" / "delsys.py")
+    mod = importlib.util.module_from_spec(spec)
+    spec.loader.exec_module(mod)
+    return mod
+
+
+@pytest.mark.skipif(not ARM_F.exists(), reason="set MYOGEN_ARM_F_NPZ to the arm-F result")
+def test_scorer_reproduces_arm_f():
+    d = np.load(ARM_F, allow_pickle=True)
+    res = _delsys().score(d["mono_emg"], float(d["fs"]), float(d["ramp"]), float(d["hold"]))
+    assert set(res) == {"cc", "adj"}
+    assert round(res["cc"][0], 3) == 0.196 and res["cc"][1] == pytest.approx(162.5, abs=0.5)
+    assert round(res["adj"][0], 3) == 0.233 and res["adj"][1] == pytest.approx(157.5, abs=0.5)
+
+
+def test_scorer_rejects_wrong_shape():
+    with pytest.raises(ValueError, match="time, sensors, 3, 3"):
+        _delsys().score(np.zeros((1000, 3, 3)), 10240.0, 2.0, 4.0)
diff --git a/tests/test_profiles_minors.py b/tests/test_profiles_minors.py
new file mode 100644
index 00000000..56a9da2b
--- /dev/null
+++ b/tests/test_profiles_minors.py
@@ -0,0 +1,75 @@
+"""Deferred minors of the named-profiles review: note text, warning location, note keys, parent version."""
+import warnings
+
+import numpy as np
+import pytest
+import quantities as pq
+import yaml
+
+from myogen import load_profile
+from myogen.profiles import ProfileNote
+from myogen.profiles._card import CARD_DIR
+from myogen.profiles._notes import reset_notes
+
+
+def _notes(rec):
+    return [w for w in rec if issubclass(w.category, ProfileNote)]
+
+
+def _card_without_note(tmp_path, field="thresholds.recruitment_range"):
+    raw = yaml.safe_load((CARD_DIR / "fdi_semg.yaml").read_text(encoding="utf-8"))
+    section, name = field.split(".")
+    del raw["sections"][section][name]["note"]
+    assert raw["sections"][section][name]["calibrated"] is True
+    path = tmp_path / "no_note.yaml"
+    path.write_text(yaml.safe_dump(raw), encoding="utf-8")
+    return load_profile(path)
+
+
+def test_calibrated_field_without_note_still_warns(tmp_path):
+    reset_notes()
+    card = _card_without_note(tmp_path)
+    with warnings.catch_warnings(record=True) as rec:
+        warnings.simplefilter("always")
+        card.thresholds(recruitment_range=30.0)
+    texts = [str(w.message) for w in _notes(rec)]
+    assert len(texts) == 1
+    assert "thresholds.recruitment_range" in texts[0] and "calibration" in texts[0]
+    assert "85" in texts[0]  # the card's calibration.best_trial
+
+
+def test_note_points_at_the_callers_line():
+    reset_notes()
+    card = load_profile()
+    with warnings.catch_warnings(record=True) as rec:
+        warnings.simplefilter("always")
+        card.muscle_kwargs(skin_thickness__mm=2.0 * pq.mm)  # muscle_kwargs -> _muscle_kwargs -> _take -> _override
+        card.thresholds(recruitment_range=30.0)  # thresholds -> _take -> _override
+    notes = _notes(rec)
+    assert len(notes) == 2
+    assert all(w.filename == __file__ for w in notes), [w.filename for w in notes]
+
+
+def test_same_field_on_two_cards_notes_twice():
+    reset_notes()
+    with warnings.catch_warnings(record=True) as rec:
+        warnings.simplefilter("always")
+        load_profile("fdi_semg").thresholds(recruitment_range=30.0)
+        load_profile("fdi_semg.published").thresholds(recruitment_range=30.0)
+    texts = [str(w.message) for w in _notes(rec)]
+    assert len(texts) == 2
+    assert texts[0].startswith("fdi_semg:") and texts[1].startswith("fdi_semg.published:")
+
+
+def test_record_carries_the_parent_version():
+    parent = load_profile("fdi_semg")
+    assert load_profile("fdi_semg.published").record()["parent_version"] == parent.version
+    assert parent.record()["parent_version"] is None
+
+
+@pytest.mark.parametrize("phrase", ["thresholds of example 01"])
+def test_example_04_does_not_claim_example_01_thresholds(phrase):
+    from pathlib import Path
+
+    text = (Path(__file__).resolve().parents[1] / "examples/01_basic/04_simulate_muscle.py").read_text(encoding="utf-8")
+    assert phrase not in text
diff --git a/tests/test_profiles_notes.py b/tests/test_profiles_notes.py
new file mode 100644
index 00000000..d3ae80c2
--- /dev/null
+++ b/tests/test_profiles_notes.py
@@ -0,0 +1,56 @@
+"""ProfileNote: card notes once per run; classes built by hand print nothing."""
+import warnings
+
+import numpy as np
+import pytest
+import quantities as pq
+
+from myogen.profiles import ProfileNote
+from myogen.profiles._notes import emit_note, reset_notes
+from myogen.simulator import Muscle
+
+THRESHOLDS = np.linspace(0.02, 1.0, 10)
+
+
+def _notes(record):
+    return [w for w in record if issubclass(w.category, ProfileNote)]
+
+
+def test_note_prints_once_per_key():
+    reset_notes()
+    with warnings.catch_warnings(record=True) as rec:
+        warnings.simplefilter("always")
+        emit_note("surface_emg.iap_lambda__per_mm", "lambda note")
+        emit_note("surface_emg.iap_lambda__per_mm", "lambda note")
+        emit_note("muscle.skin_thickness__mm", "skin note")
+    assert [str(w.message) for w in _notes(rec)] == ["lambda note", "skin note"]
+
+
+def test_hand_built_classes_do_not_hint():
+    """Classes built by hand say nothing about profile cards (a card is FDI-specific)."""
+    from myogen.simulator import ForceModel
+
+    reset_notes()
+    with warnings.catch_warnings(record=True) as rec:
+        warnings.simplefilter("always")
+        Muscle(recruitment_thresholds=THRESHOLDS)
+        ForceModel(recruitment_thresholds=THRESHOLDS, recording_frequency__Hz=2048 * pq.Hz)
+    assert _notes(rec) == []
+
+
+def test_notes_can_be_silenced():
+    reset_notes()
+    with warnings.catch_warnings(record=True) as rec:
+        warnings.simplefilter("always")
+        warnings.simplefilter("ignore", ProfileNote)
+        emit_note("a", "text")
+        Muscle(recruitment_thresholds=THRESHOLDS)
+    assert _notes(rec) == []
+
+
+def test_notes_can_be_errors():
+    reset_notes()
+    with warnings.catch_warnings():
+        warnings.simplefilter("error", ProfileNote)
+        with pytest.raises(ProfileNote):
+            emit_note("b", "text")
diff --git a/tests/test_profiles_review_fixes.py b/tests/test_profiles_review_fixes.py
new file mode 100644
index 00000000..927060c7
--- /dev/null
+++ b/tests/test_profiles_review_fixes.py
@@ -0,0 +1,209 @@
+"""Final-review fixes: overrides are applied where they are recorded, carry across builders, and records are exact."""
+import copy
+import warnings
+
+import numpy as np
+import pytest
+import quantities as pq
+import yaml
+
+from myogen import load_profile
+from myogen.profiles import ProfileError, ProfileNote
+from myogen.profiles._card import CARD_DIR
+from myogen.profiles._notes import reset_notes
+
+CV = "motor_unit_conduction_velocities__m_per_s"
+
+
+class _StubForce:
+    def __init__(self, n):
+        self.n = n
+
+    def generate_force(self, block, verbose=False):
+        return np.full((50, 1), 100.0) * pq.dimensionless
+
+
+@pytest.fixture
+def stub_drive(monkeypatch):
+    """Capture what simulate() hands to the drive, without running NEURON."""
+    import neo
+
+    from myogen.profiles import _builders
+    from myogen.simulator.neuron import descending_drive_network as ddn
+
+    seen = {}
+
+    def fake_excitation(segments, threshold_hz, max_hz, timestep_ms=0.1, noise_seed=42):
+        seen.update(threshold_hz=threshold_hz, timestep_ms=timestep_ms)
+        return np.zeros(10), np.zeros(10)
+
+    def fake_run(thresholds, drive_hz, dd_params, total_ms, timestep_ms=0.1, config_file="alpha_mn_default.yaml"):
+        seen.update(n_units=len(thresholds), run_timestep_ms=timestep_ms, config_file=config_file)
+        block = neo.Block()
+        block.segments.append(neo.Segment())
+        return block
+
+    monkeypatch.setattr(ddn, "excitation_to_drive_hz", fake_excitation)
+    monkeypatch.setattr(ddn, "run_descending_drive", fake_run)
+    monkeypatch.setattr(_builders.CardBuilders, "force_model",
+                        lambda self, **kw: (seen.__setitem__("force_n", len(self.thresholds())), _StubForce(0))[1])
+    return seen
+
+
+def test_simulate_applies_drive_and_force_overrides(stub_drive):
+    card = load_profile()
+    res = card.simulate([(1, 0, 25)], timestep_ms=0.05, threshold_hz=90.0, mvc_model_units=1000.0,
+                        recording_frequency__Hz=1024.0, config_file="other.yaml")
+    assert stub_drive["timestep_ms"] == stub_drive["run_timestep_ms"] == 0.05
+    assert stub_drive["threshold_hz"] == 90.0
+    assert stub_drive["config_file"] == "other.yaml"
+    assert res.fs_force__Hz == 1024.0
+    assert res.force_pct_mvc[0] == pytest.approx(10.0)  # 100 / 1000 * 100
+    assert set(res.profile["overrides"]) == {"drive.timestep_ms", "drive.threshold_hz", "drive.mvc_model_units",
+                                             "force.recording_frequency__Hz", "motor_neurons.config_file"}
+
+
+def test_simulate_uses_overridden_excitation_curve(stub_drive):
+    card = load_profile()
+    card.simulate([(1, 0, 25)], excitation_at_pct_mvc={0: 0.0, 25: 0.5, 100: 1.0})
+    assert card.excitation_for_pct_mvc(25) == pytest.approx(0.5)
+
+
+def test_overrides_carry_to_later_builds(stub_drive):
+    card = load_profile()
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", ProfileNote)
+        card.thresholds(n_units=60)
+        assert len(card.muscle_kwargs()["recruitment_thresholds"]) == 60
+        card.simulate([(1, 0, 25)])
+    assert stub_drive["n_units"] == stub_drive["force_n"] == 60
+    card.thresholds(n_units=120)  # back to the card value: no longer an override
+    assert card.overrides() == {}
+
+
+def test_override_in_another_unit_is_converted():
+    reset_notes()
+    card = load_profile()
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", ProfileNote)
+        kw = card.muscle_kwargs(radius__mm=0.008 * pq.m)
+    assert kw["radius__mm"].rescale(pq.mm).magnitude == pytest.approx(8.0)
+    assert card.overrides() == {"muscle.radius__mm": pytest.approx(8.0)}
+
+
+def test_quantity_on_a_unitless_field_raises():
+    with pytest.raises(TypeError, match="thresholds.recruitment_range"):
+        load_profile().thresholds(recruitment_range=50 * pq.mm)
+
+
+def test_dispersion_equal_to_card_is_not_an_override():
+    from myogen.simulator.core.emg.surface.fiber_dispersion import SurfaceFiberDispersion
+    card = load_profile()
+    same = card.surface_emg_kwargs()["fiber_dispersion"]
+    card.surface_emg_kwargs(fiber_dispersion=copy.deepcopy(same))
+    assert card.overrides() == {}
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", ProfileNote)
+        kw = card.surface_emg_kwargs(fiber_dispersion=SurfaceFiberDispersion())
+    assert isinstance(kw["fiber_dispersion"], SurfaceFiberDispersion)
+    assert list(card.overrides()) == ["surface_emg.fiber_dispersion"]
+    assert card.surface_emg_kwargs()["fiber_dispersion"] == SurfaceFiberDispersion()  # carried, rebuilt
+
+
+def test_unit_cvs_equal_to_card_are_not_an_override():
+    card = load_profile()
+    cvs = card.muscle_kwargs()[CV]
+    card.muscle_kwargs(**{CV: np.array(cvs)})
+    assert card.overrides() == {}
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", ProfileNote)
+        kw = card.muscle_kwargs(**{CV: np.full(120, 4.0)})
+    np.testing.assert_array_equal(kw[CV], np.full(120, 4.0))
+    assert list(card.overrides()) == ["muscle.unit_cv"]
+    np.testing.assert_array_equal(card.muscle_kwargs()[CV], np.full(120, 4.0))  # carried
+
+
+def test_pass_through_keywords_are_recorded_on_the_object():
+    card = load_profile()
+    fm = card.force_model(longest_duration_rise_time__ms=60.0 * pq.ms)
+    assert card.overrides() == {}
+    assert fm.profile_record["keywords"] == {"longest_duration_rise_time__ms": "60.0 ms"}
+
+
+def test_internal_sampling_rate_stays_the_library_default():
+    kw = load_profile().surface_emg_kwargs(sampling_frequency__Hz=2048.0 * pq.Hz)
+    assert "internal_sampling_frequency__Hz" not in kw
+
+
+def test_curve_without_ends_extends_to_rest_and_mvc(tmp_path):
+    path = tmp_path / "short_curve.yaml"
+    path.write_text("name: short_curve\nparent: fdi_semg\nversion: 1\nsections:\n  drive:\n"
+                    "    excitation_at_pct_mvc: {value: {10: 0.1, 50: 0.4}, source: assumption}\n", encoding="utf-8")
+    card = load_profile(path)
+    assert card.excitation_for_pct_mvc(80) == pytest.approx(0.4 + 0.6 * 30 / 50)
+    assert card.excitation_for_pct_mvc(5) == pytest.approx(0.05)
+
+
+def test_card_missing_a_required_field_fails_at_load(tmp_path):
+    raw = yaml.safe_load((CARD_DIR / "fdi_semg.yaml").read_text(encoding="utf-8"))
+    del raw["sections"]["muscle"]["radius__mm"]
+    path = tmp_path / "broken.yaml"
+    path.write_text(yaml.safe_dump(raw), encoding="utf-8")
+    with pytest.raises(ProfileError, match="muscle.radius__mm"):
+        load_profile(path)
+    del raw["sections"]["muscle"]
+    raw["sections"]["thresholds"] = raw["sections"]["thresholds"]
+    del raw["seeds"]["cv_scatter"]
+    path.write_text(yaml.safe_dump(raw), encoding="utf-8")
+    with pytest.raises(ProfileError, match="seeds.cv_scatter"):
+        load_profile(path)
+
+
+def _small_muscle():
+    from myogen import set_random_seed
+    from myogen.simulator import Muscle, RecruitmentThresholds
+    set_random_seed(42)
+    thresholds, _ = RecruitmentThresholds(N=12, recruitment_range__ratio=30, mode="fuglevand")
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", ProfileNote)
+        muscle = Muscle(thresholds, radius__mm=6 * pq.mm, length__mm=60 * pq.mm,
+                        fiber_density__fibers_per_mm2=6 * pq.mm**-2, mean_fiber_length__mm=50 * pq.mm,
+                        var_fiber_length__mm=2 * pq.mm, radius_bone__mm=0 * pq.mm, fat_thickness__mm=1 * pq.mm)
+    muscle.distribute_innervation_centers()
+    muscle.generate_muscle_fiber_centers(verbose=False)
+    muscle.assign_mfs2mns(verbose=False, n_jobs=1)
+    return muscle
+
+
+def test_card_surface_emg_muaps_do_not_depend_on_what_ran_in_between():
+    from myogen import set_random_seed
+    from myogen.simulator import SurfaceElectrodeArray
+    from myogen.utils.neo import signal_to_grid
+
+    muscle = _small_muscle()
+    card = load_profile()
+
+    def muaps(between):
+        array = SurfaceElectrodeArray(1, 1, 10 * pq.mm, 2.5 * pq.mm, center_point__mm_deg=(12 * pq.mm, 0 * pq.deg),
+                                      bending_radius__mm=8.3 * pq.mm, differentiation_mode="bipolar_longitudinal")
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", ProfileNote)
+            semg = card.surface_emg(muscle, [array], sampling_frequency__Hz=2048 * pq.Hz,
+                                    internal_sampling_frequency__Hz=2048 * pq.Hz, MUs_to_simulate=[4, 11],
+                                    fiber_dispersion=None, iap_kernel_length__mm=150.0)
+        between()
+        block = semg.simulate_muaps(n_jobs=1, verbose=False, use_gpu=False)
+        return np.array([np.asarray(signal_to_grid(s.analogsignals[0].rescale(pq.mV))) for s in block.groups[0].segments])
+
+    np.testing.assert_array_equal(muaps(lambda: None), muaps(lambda: set_random_seed(999)))
+
+
+@pytest.mark.parametrize("name", ("fdi_semg", "fdi_semg.published"))
+def test_cited_keys_carry_a_locator(name):
+    import re
+    for path, field in load_profile(name).fields().items():
+        src = field.source.strip()
+        if src.startswith(("assumption", "derived")):
+            continue
+        for key, locator in re.findall(r"\[([a-z]+[0-9]{4}[a-z]?)([^\]]*)\]", src):
+            assert re.search(r"\b(p\.|Table|Fig\.|eq\.)", locator), f"{path}: [{key}] has no locator"
diff --git a/tests/test_profiles_seed_scope.py b/tests/test_profiles_seed_scope.py
new file mode 100644
index 00000000..e8fe1d1d
--- /dev/null
+++ b/tests/test_profiles_seed_scope.py
@@ -0,0 +1,124 @@
+"""Card builders seed only their own build: the caller's generator and seed are given back, nothing is printed."""
+import warnings
+
+import numpy as np
+import pytest
+import quantities as pq
+
+from myogen import get_random_generator, get_random_seed, load_profile, set_random_seed
+from myogen.profiles import ProfileNote
+from tests.test_profiles_review_fixes import _small_muscle
+
+
+def _expected_draw(seed):
+    return np.random.default_rng(seed).integers(0, 10**9, 5)
+
+
+def _assert_users_stream(seed):
+    assert get_random_seed() == seed
+    np.testing.assert_array_equal(get_random_generator().integers(0, 10**9, 5), _expected_draw(seed))
+
+
+def _semg(card, muscle, **kw):
+    from myogen.simulator import SurfaceElectrodeArray
+
+    array = SurfaceElectrodeArray(1, 1, 10 * pq.mm, 2.5 * pq.mm, center_point__mm_deg=(12 * pq.mm, 0 * pq.deg),
+                                  bending_radius__mm=8.3 * pq.mm, differentiation_mode="bipolar_longitudinal")
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", ProfileNote)
+        return card.surface_emg(muscle, [array], sampling_frequency__Hz=2048 * pq.Hz,
+                                internal_sampling_frequency__Hz=2048 * pq.Hz, MUs_to_simulate=[4, 11],
+                                fiber_dispersion=None, iap_kernel_length__mm=150.0, **kw)
+
+
+def _muaps(semg):
+    from myogen.utils.neo import signal_to_grid
+
+    block = semg.simulate_muaps(n_jobs=1, verbose=False, use_gpu=False)
+    return np.array([np.asarray(signal_to_grid(s.analogsignals[0].rescale(pq.mV))) for s in block.groups[0].segments])
+
+
+def test_thresholds_leave_the_users_generator_and_seed():
+    set_random_seed(123)
+    load_profile().thresholds()
+    _assert_users_stream(123)
+
+
+def test_builders_do_not_print_the_seed(capsys):
+    set_random_seed(123)
+    capsys.readouterr()
+    load_profile().thresholds()
+    assert "Random seed set" not in capsys.readouterr().out
+
+
+def test_surface_emg_build_and_muaps_leave_the_users_generator_and_seed():
+    muscle = _small_muscle()
+    card = load_profile()
+    set_random_seed(123)
+    semg = _semg(card, muscle)
+    _assert_users_stream(123)
+    set_random_seed(123)
+    _muaps(semg)
+    _assert_users_stream(123)
+
+
+def test_pool_build_leaves_the_users_generator_and_seed():
+    from myogen.utils.nmodl import load_nmodl_mechanisms
+
+    load_nmodl_mechanisms()
+    set_random_seed(123)
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", ProfileNote)
+        pools = load_profile().motor_neuron_pool()
+    _assert_users_stream(123)
+    assert pools.profile_record["seed"] == 42
+
+
+def test_surface_emg_seed_keyword_gives_another_realisation_and_is_recorded():
+    muscle = _small_muscle()
+    card = load_profile()
+    default, other = _semg(card, muscle), _semg(card, muscle, seed=8)
+    assert default.profile_record["seed"] == card.seeds["surface_emg"]
+    assert other.profile_record["seed"] == 8
+    assert not any("seed" in path for path in card.overrides())  # a seed is not a card field
+    assert not np.array_equal(_muaps(default), _muaps(other))
+    np.testing.assert_array_equal(_muaps(default), _muaps(_semg(card, muscle, seed=card.seeds["surface_emg"])))
+
+
+def test_simulate_seed_keyword_reaches_the_drive_and_is_recorded(monkeypatch):
+    import neo
+
+    from myogen.profiles import _builders
+    from myogen.simulator.neuron import descending_drive_network as ddn
+
+    seen = {}
+
+    def fake_excitation(segments, threshold_hz, max_hz, timestep_ms=0.1, noise_seed=42):
+        seen["noise_seed"] = noise_seed
+        return np.zeros(10), np.zeros(10)
+
+    def fake_run(thresholds, drive_hz, dd_params, total_ms, timestep_ms=0.1, config_file="alpha_mn_default.yaml"):
+        seen["seed_in_run"] = get_random_seed()
+        block = neo.Block()
+        block.segments.append(neo.Segment())
+        return block
+
+    class _Force:
+        def generate_force(self, block, verbose=False):
+            return np.full((50, 1), 100.0) * pq.dimensionless
+
+    monkeypatch.setattr(ddn, "excitation_to_drive_hz", fake_excitation)
+    monkeypatch.setattr(ddn, "run_descending_drive", fake_run)
+    monkeypatch.setattr(_builders.CardBuilders, "force_model", lambda self, **kw: _Force())
+    set_random_seed(123)
+    res = load_profile().simulate([(1, 0, 25)], seed=5)
+    assert seen == {"noise_seed": 5, "seed_in_run": 5}
+    assert res.profile["seed"] == 5 and res.spikes.annotations["myogen_profile"]["seed"] == 5
+    _assert_users_stream(123)
+
+
+def test_thresholds_accept_a_seed():
+    card = load_profile()
+    np.testing.assert_array_equal(card.thresholds(seed=5), card.thresholds())  # fuglevand mode draws nothing
+    with pytest.raises(TypeError):
+        card.thresholds(seed=1.5)
diff --git a/tests/test_profiles_simulate.py b/tests/test_profiles_simulate.py
new file mode 100644
index 00000000..3109f3bd
--- /dev/null
+++ b/tests/test_profiles_simulate.py
@@ -0,0 +1,46 @@
+"""Running a card in % MVC: excitation mapping, refusal without a mapping, and a short calibrated-drive run."""
+import numpy as np
+import pytest
+
+from myogen import load_profile
+from myogen.profiles import ProfileError
+
+
+def test_excitation_for_pct_mvc():
+    card = load_profile()
+    assert card.excitation_for_pct_mvc(25) == pytest.approx(0.18818359375000002)
+    assert card.excitation_for_pct_mvc(40) == pytest.approx((0.216552734375 + 0.330078125) / 2)
+    assert card.excitation_for_pct_mvc(100) == pytest.approx(1.0)
+    assert card.excitation_for_pct_mvc(0) == 0.0
+    np.testing.assert_allclose(card.excitation_for_pct_mvc(np.array([5, 10])), [0.048095703125, 0.08164062500000001])
+    for bad in (101, -1, float("nan")):
+        with pytest.raises(ValueError):
+            card.excitation_for_pct_mvc(bad)
+
+
+def test_card_without_mvc_mapping_refuses_pct_mvc(tmp_path):
+    path = tmp_path / "no_mapping.yaml"
+    path.write_text("name: no_mapping\nparent: fdi_semg\nversion: 1\nsections:\n  drive:\n"
+                    "    excitation_at_pct_mvc: {value: null, source: assumption}\n", encoding="utf-8")
+    with pytest.raises(ProfileError, match="drive.excitation_at_pct_mvc"):
+        load_profile(path).excitation_for_pct_mvc(25)
+
+
+def test_published_variant_has_its_own_mvc_mapping():
+    card = load_profile("fdi_semg.published")
+    assert card.value("drive.mvc_model_units") != load_profile().value("drive.mvc_model_units")
+    assert 0.0 < card.excitation_for_pct_mvc(25) < 1.0
+
+
+def test_simulate_short_run():
+    """Slow-ish: 1.5 s of the calibrated 120-unit drive."""
+    card = load_profile()
+    res = card.simulate([(0.5, 0, 0), (1.0, 0, 25)])
+    trains = res.spikes.segments[0].spiketrains
+    assert len(trains) == 120
+    assert sum(len(st) for st in trains) > 0
+    assert res.fs_force__Hz == 2048.0
+    assert res.force_pct_mvc.ndim == 1 and res.force_pct_mvc.size == int(1.5 * 2048)
+    assert res.excitation.max() == pytest.approx(card.excitation_for_pct_mvc(25), rel=1e-3)
+    assert res.spikes.annotations["myogen_profile"] == card.record(seed=card.seeds["drive"])
+    assert 0.0 <= res.force_pct_mvc.max() < 60.0
diff --git a/tests/test_seed_scope_and_dependencies.py b/tests/test_seed_scope_and_dependencies.py
new file mode 100644
index 00000000..cb170d12
--- /dev/null
+++ b/tests/test_seed_scope_and_dependencies.py
@@ -0,0 +1,31 @@
+"""``myogen._seed_scope`` error handling, and the dependencies ``import myogen`` needs."""
+import tomllib
+from pathlib import Path
+
+import pytest
+
+import myogen
+import myogen.simulator.neuron.cells as cells
+
+REPO = Path(__file__).resolve().parents[1]
+
+
+def test_seed_scope_raises_reset_errors_and_gives_the_generator_back(monkeypatch):
+    myogen.set_random_seed(5)
+    before = myogen.get_random_generator()
+
+    def broken():
+        raise RuntimeError("reset failed")
+
+    monkeypatch.setattr(cells, "reset_cell_id_counters", broken)
+    with pytest.raises(RuntimeError, match="reset failed"):
+        with myogen._seed_scope(1):
+            pass
+    assert myogen.get_random_generator() is before
+    assert myogen.get_random_seed() == 5
+
+
+def test_pyyaml_is_declared():
+    """``import myogen`` loads the profile cards, which need PyYAML."""
+    deps = tomllib.loads((REPO / "pyproject.toml").read_text(encoding="utf-8"))["project"]["dependencies"]
+    assert any(d.lower().startswith("pyyaml") for d in deps)
diff --git a/uv.lock b/uv.lock
index fdac1fc1..ea3daddd 100644
--- a/uv.lock
+++ b/uv.lock
@@ -984,6 +984,7 @@ dependencies = [
     { name = "optuna" },
     { name = "pandas" },
     { name = "psutil" },
+    { name = "pyyaml" },
     { name = "quantities" },
     { name = "scienceplots" },
     { name = "scikit-fmm" },
@@ -1036,6 +1037,7 @@ requires-dist = [
     { name = "pandas", specifier = ">=2.0" },
     { name = "psutil", specifier = ">=5.9" },
     { name = "pynwb", marker = "extra == 'nwb'", specifier = ">=2.8.0" },
+    { name = "pyyaml", specifier = ">=6.0" },
     { name = "quantities", specifier = ">=0.16.2" },
     { name = "scienceplots", specifier = ">=2.1.1" },
     { name = "scikit-fmm", specifier = ">=2025.1.29" },