Surface EMG: corrections against Farina 2004, FDI layer defaults, published end-plates

CHANGELOG.md

@@ -7,6 +7,55 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+This is a breaking release. It corrects the surface-EMG pipeline against Farina, Mesin, Martina & Merletti (2004). The problems were found through feedback from @joaohbbittar; thank you. Both the default route and `execution=` are affected.
+
+### Breaking changes
+
+This release changes default outputs; regenerate anything calibrated against earlier versions.
+
+- sEMG outputs change in amplitude, shape, duration and spectrum. Bipolar outputs from 0.5.0–0.13.0 are invalid. Durations, mean/median frequencies and conduction-velocity estimates from those versions are off by about 2× and must be regenerated.
+- Absolute sEMG amplitudes change by one to two orders of magnitude (amplitude no longer depends on conduction velocity). Anything calibrated in µV or mV against earlier outputs must be regenerated. To keep the earlier resistive (measured) skin layer, pass `skin_conductivity__S_per_m=4.88e-4 * pq.S / pq.m` to `Muscle`. To keep the earlier layers, pass `skin_thickness__mm=1.29 * pq.mm, fat_thickness__mm=0.3 * pq.mm`.
+- Pickled `SurfaceEMG` objects from earlier versions lack the new fibre properties; create them again.
+- The `execution=` MUAP disk-cache key changes, so templates cached by earlier versions are recomputed.
+- iEMG outputs change through the narrower end-plate spreads. Pass `IntramuscularEMG(..., nmj_coefficients="reference_code")` to reproduce earlier end-plates.
+
+### Fixed
+
+- **Longitudinal bipolar channels.** The angular integral B(k_z) is now complex. Keeping only its real part made `bipolar_longitudinal` channels about 10⁻⁷ of their true value and distorted rotated arrays. A bipolar channel now equals the difference of the two monopolar electrodes, on the default, `execution=` and `use_unified=True` routes.
+- **MUAP time axis.** MUAP templates are labelled with the rate of the grid they were computed on (window `K / v`). The previous label was 2× too high, so MUAPs looked half as long, spectra sat an octave high, and conduction velocity estimated from the output was 2× the value set.
+- **Rosenfalck IAP length scale.** The IAP is evaluated on the simulation's own z axis (λ = 1 mm, as published), not on `z / 2`. The stretched IAP was 2× too long; it now agrees with the independent time-domain path (`simulate_fiber_hybrid`).
+- **Fibre semi-lengths.** `L1` (towards +z) and `L2` were swapped, which moved both tendons by 2·z_i with the end-plate. Fibres now span `±fiber_length / 2` with the end-plate at z_i (`fiber_semi_lengths`).
+- **Fibre timing within a motor unit.** Each fibre's potential was circularly shifted onto its own peak before summation, misaligning fibres (up to ~60 % MUAP error on bipolar channels). All fibres now share the firing instant at the window centre.
+- **Mixed electrode-array shapes.** On the default route each array's MUAPs have that array's own grid shape; previously every array was broadcast into the first array's shape.
+- **Absolute sEMG amplitude.** Three scaling errors are corrected; together they made MUAPs labelled mV roughly 50-100× too large. A single fibre now matches an exact solution for a homogeneous insulated cylinder within 2 %.
+  - The Fourier transforms dropped their `dz` and `Fs/N` factors, so every potential was divided by the conduction velocity.
+  - The source was the bare Rosenfalck slope with an ad-hoc `A = 96/10`. It is now the core-conductor line current `σᵢ·π·(d/2)²·∂²Vm/∂z²` with the published `A = 96 mV/mm³`.
+  - The default skin conductivity was 4.88·10⁻⁴ S/m (measured composite skin, Botelho et al. 2019), which raised amplitudes 3.5-6×. It is now 1 S/m, as in Farina et al. (2004). That value comes from Roeleveld et al. (1997), who chose it so their model matched how biceps MUAPs fall off with depth and around the arm; it is an effective value, not a measured skin conductivity. MyoGen keeps it as a calibration: with the low-frequency skin conductivity of the IT'IS tissue database (0.15 S/m, the mean of 7 studies ranging 0.06–0.25 S/m; Baumgartner et al. 2025) the FDI plateau RMS is about 2.5-4× above published ranges (`benchmarks/semg_dispersion/`). A check against Roeleveld's fall-off data (`benchmarks/semg_roeleveld1997/`) was inconclusive and does not support 1 S/m as a skin property.
+- **Default skin and fat thickness (FDI).** Skin 1.29 → 1.2 mm (hand dorsum median, Li et al. 2018) and fat 0.3 → 1.3 mm, so the total depth to muscle is 2.5 mm, the low end of the FDI depth measured by Nishikawa et al. (2024); the fat value is derived from that total, not measured. Surface amplitudes from the default `Muscle` fall by about 40 %.
+- **Failed motor units raise.** On the default `SurfaceEMG` route an exception while simulating one motor unit was logged and replaced by an all-zero MUAP, which passes for a valid, silent unit. The exception is now re-raised with the unit and array index; the `ExecutionPolicy` route already raised.
+- **Angular sampling overflow raises.** With `sampling_points_in_theta_domain` ≥ 256 the scaled modified Bessel functions leave the float64 range at the highest angular orders (I_n underflows to 0, K_n overflows). The volume conductor replaced the overflowing entries with zero, which gave a singular system on the CPU and NaN MUAPs on the GPU. It now raises a `ValueError` that names the setting.
+
+### Added
+
+- `SurfaceEMG(fiber_diameter__um=50, 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:
+  - per-fibre end-plates from the two-layer nerve model of Akhmadeev et al. (2020), shared with intramuscular EMG;
+  - axon-branch conduction delays;
+  - motor-unit conduction velocity linear in recruitment threshold (Del Vecchio et al. 2017), and a spread of fibre velocities within each unit, each fibre propagating at its own velocity;
+  - a source scaled by fibre diameter (from the Nandedkar & Stålberg 1983 velocity–diameter line);
+  - tendon-end stagger: each fibre end in a band, independently or following a smooth random field across the cross-section so that neighbouring fibres end close together.
+
+  The realised spreads are reported per motor unit in `SurfaceEMG.fiber_dispersion_report`. The default is off and leaves outputs bit-identical. `SurfaceFiberDispersion.literature_candidate()` gives a preset whose values come from the literature or are marked as assumptions. The `execution=` route and `use_unified=True` refuse an active configuration.
+- `IntramuscularEMG(nmj_coefficients=...)` selects the end-plate spreads (see Changed).
+- `benchmarks/semg_dispersion/`: paired sweeps of these options on a limb-like cylinder and on the default (FDI) muscle.
+
+### Changed
+
+- **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.
+- `Muscle.var_fiber_length__mm` is documented as the half-width of a uniform fibre-length variation (the FDI value of Jacobson et al. 1992), not a variance. Behaviour is unchanged.
+
 ## [0.13.0] - 2026-10-28
 
 Thanks to @joaohbbittar for the detailed report in #26, this release corrects the Fuglevand force model and the Hill muscle activation.

benchmarks/semg_dispersion/README.md

@@ -0,0 +1,66 @@
+# Surface-EMG fibre-dispersion sweep
+
+Paired simulations behind the `SurfaceFiberDispersion` options and the
+`SurfaceFiberDispersion.literature_candidate()` preset. Each arm reuses the same
+muscle, fibres, spike trains and electrode as its scenario's baseline and
+changes only the dispersion configuration, so MUAPs are compared unit by unit.
+Fibre and spike hashes are stored in every output and checked when scoring.
+
+## Scenarios
+
+| Scenario | Muscle | Electrode |
+|---|---|---|
+| `limb` | cylinder r 15 mm, bone 3 mm, fat 4 mm, 70 ± 5 mm fibres, 350 fibres/mm², 100 MUs | bipolar, 5 mm contacts, IED 10 mm, 20 mm from the end-plate band |
+| `fdi` | default `Muscle` (tuned for the first dorsal interosseous), 120 MUs | bipolar, 2 mm contacts, IED 5 mm, 8 mm from the end-plate band |
+
+Both: 30 % excitation trapezoid, Fuglevand thresholds (50-fold range) recruited up
+to 70 % excitation, 8–35 Hz rate coding, ISI CV 0.2, 2048 Hz. Scoring applies a
+20–450 Hz band-pass and no noise.
+
+## Arms
+
+`arms.json` holds every arm per scenario (`baseline` is the undispersed model).
+Arms run before the published end-plate coefficients became the default pin
+`"nmj_coefficients": "reference_code"`, so they reproduce. `_skin_thickness__mm`
+is a harness key for the skin layer.
+
+Arms with `_spike_file` use NEURON spike trains of the calibrated FDI pool
+(spike times per unit and the pool thresholds), stored under `spikes/`; the
+path is relative to this directory.
+
+## Running
+
+```bash
+# from the repository root; PY selects the interpreter. Check that the
+# recorded myogen path (key "myogen" in each output) is this checkout; set
+# PYTHONPATH to the checkout if another one is installed.
+PY=python bash benchmarks/semg_dispersion/run_arms.sh fdi out/fdi           # all FDI arms
+PY=python bash benchmarks/semg_dispersion/run_arms.sh limb out/limb baseline candidate_published
+python benchmarks/semg_dispersion/score.py fdi out/fdi          # scorecard.tsv, semg_dispersion_sweep.pdf
+python benchmarks/semg_dispersion/plot_signals.py fdi out/fdi   # semg_dispersion_signals.pdf
+python benchmarks/semg_dispersion/plot_amplitude.py out/limb out/fdi out/amplitude.pdf
+```
+
+A limb arm takes about 20–25 min on CPU and an FDI arm about 5 min; the runner
+does three at a time and skips arms already present.
+
+## Findings (2026-10-31)
+
+- Within-unit fibre CV spread is the strongest lever on MUAP peak-to-peak; a
+  CV difference between units only rescales time.
+- In the limb scenario, dispersion mostly blunts the end-of-fibre phase of small
+  and mid units; plateau RMS falls by about 6 %. In the FDI scenario the short,
+  sharp MUAPs of large superficial units lose more, and plateau RMS falls by
+  about 17 %.
+- Tendon ends: a 5 mm band of independent ends lowers median MUAP
+  peak-to-peak by about 17 % (limb); ends correlated over 5 mm across the
+  cross-section, so that neighbours end close together, by about 11 %.
+- FDI skin thickness 1.29 → 3 mm lowers RMS by about 60 %; it is the strongest
+  amplitude lever found so far.
+- Skin conductivity matters more than thickness: at skin 1.2 mm and fat 1.3 mm,
+  1 → 0.15 S/m raises FDI plateau RMS about 3.3×. With 1 S/m and these
+  thicknesses the all-effects model falls inside published FDI ranges; with
+  0.15 S/m, the low-frequency skin conductivity of the IT'IS tissue database
+  (mean of 7 studies, 0.06–0.25 S/m), it is about 2.5–4× too high.
+
+`benchmarks/semg_versions` compares outputs before and after the unreleased surface-EMG fixes.