Repository navigation
Surface EMG: corrections against Farina 2004, FDI layer defaults, published end-plates - #40
Merged
Merged
Conversation
- B(k_z) keeps its imaginary part in the Numba kernel, the NumPy fallback and the phi kernel; longitudinal bipolar outputs were previously lost. - Remove the ad-hoc `z /= 2` IAP stretch: the Rosenfalck source now uses the same mm grid as k_z (matches the time-domain hybrid path). - MUAP sampling rate follows the internal grid (muap_sampling_rate__Hz), shared by the default and execution routes. - Drop the per-fibre peak re-centring, which erased propagation delays and fibre-to-fibre timing within a motor unit. - L1/L2 are the +z/-z semi-lengths from the end-plate (fiber_semi_lengths). - Each electrode array keeps its own output shape. - Bump the execution-route cache key so stale templates are recomputed.
… with 1/v The source FFT was normalised by 1/N and the inverse by N, dropping the dz (forward) and Fs/N (inverse) factors of the continuous transforms. Their product is dz * Fs = v, so every potential was divided by the conduction velocity. Use dz on the forward transform and Fs on the inverse. New test compares a single fibre against an exact, independently written Fourier-Bessel solution for a homogeneous insulated cylinder (both routes, two velocities, two geometries): now within 2 %, before it was off by 1/v.
The fibre source was the bare Rosenfalck slope with an ad-hoc A = 96/10 and no fibre properties, so outputs labelled mV were not potentials. Use the published A = 96 mV/mm^3 and the core-conductor line current density sigma_i * pi * (d/2)^2 * d2Vm/dz2; with S/m and mm the result is in mV. SurfaceEMG gains fiber_diameter__um (default 50) and intracellular_conductivity__S_m (default 1.01), passed to every fibre on both routes and included in the execution cache key. The experimental hybrid path (use_unified=True) keeps its own D1 source and is unchanged. Tests: exact-solution comparison now uses the physical source; amplitude scales with d^2 and sigma_i; both routes forward the parameters.
…eld 1997) The previous default, 4.88e-4 S/m (measured composite skin, Botelho et al. 2019, tissue data at 150 Hz), raised surface amplitudes 3.5-6x. 1 S/m is the value used by Farina et al. (2004), taken from Roeleveld et al. (1997), who chose it so their layered cylinder matched how biceps MUAPs fall off with depth and around the arm. It is an effective parameter of the cylinder model, not a measured skin conductivity. Pass 4.88e-4 explicitly to reproduce the old default. The execution cache key already includes the skin conductivity. CHANGELOG: amplitude fixes, new SurfaceEMG fibre parameters, upgrade notes.
…mogeneous solution
nmj_branch_schedule, two_layer_nmj and nerve_delays__ms hold the Konstantin et al. (2020) branch layout, end-plate positions and axon path lengths that MotorUnitSim used inline. Draw order is unchanged, so intramuscular EMG reproduces its previous end-plates and paths exactly.
…e grid v_grid fixes the (k_z, k_t) grid and window; the fibre velocity enters only the source term, whose IAP spectrum is read at k_t / v. Fibres at different velocities then sum on one time axis. delay_fibre_signal applies a fractional-sample delay. v_grid defaults to v, which leaves every existing call unchanged.
…ocities SurfaceFiberDispersion switches each effect on independently: two-layer end-plate scatter (full or within each motor unit), axon-branch delays, size-ordered and within-unit fibre velocity spread, and diameter-scaled source strength. Draws use their own seeded generators; the default (off) is bit-identical to the undispersed model. The realised spread is reported per motor unit. The execution route and use_unified refuse an active configuration.
…iameter line Add fiber_cv_threshold_slope (Del Vecchio et al. 2017: motor-unit CV rises linearly with recruitment threshold). Read fibre diameter from the Nandedkar & Stalberg (1983) line cv = 0.95 m/s + 0.05 m/s/um * d instead of an unsourced slope and reference diameter, and report the mean diameter. Drop an unverified citation for the within-unit CV spread.
Konstantin et al. (2020, eq. 5-6) print a_mu = 1, b_mu = 2.5, a_sigma = 0.25, b_sigma = 1 mm; the reference simulator code (and MyoGen's default) uses 1.5, 4, 0.5, 1.5. Offer both in nmj_branch_schedule and SurfaceFiberDispersion, default unchanged.
Intramuscular and surface EMG now use the z spreads printed in Konstantin et al. (2020, eq. 5-6). Within-unit end-plate SD falls from a median of about 1.3 mm to 0.8 mm (largest units about 1.9 mm instead of 2.9 mm), in line with simulator conventions of about 1 mm and the 2-3.5 mm upper bounds from single-fibre and surface-array studies. The reference-code values remain available as coefficients="reference_code".
tendon_end_band__mm places each fibre end independently, uniform over a band about +-mean_fiber_length/2, instead of the length variation that moves both ends together. Off by default; the dispersion realisation now carries each fibre's ends and the surface model reads L1/L2 from them.
tendon_end_correlation__mm makes each fibre end follow a smooth Gaussian random field over the cross-section (one per tendon, squared-exponential covariance via random Fourier features) plus a per-fibre jitter, keeping the total SD of the uniform band. Neighbouring fibres then end close to each other. 0 keeps the independent band. Off by default.
SurfaceFiberDispersion.literature_candidate() collects the values chosen from the literature: published two-layer end-plates, 5/2 m/s nerve branches, threshold-linear unit CV (Del Vecchio 2017), 5 % within-unit CV and correlated tendon ends (the last two are assumptions). The diameter-scaled source stays off. Opt-in.
IntramuscularEMG(nmj_coefficients="reference_code") restores the end-plate spreads used up to 0.13.0; the default is the published Konstantin et al. (2020) values.
benchmarks/semg_dispersion: paired sweeps of the SurfaceFiberDispersion options on a limb-like cylinder and the default (FDI) muscle; arms run before the published end-plate default pin reference_code coefficients and reproduce the stored outputs bit for bit. benchmarks/semg_versions: the before/after comparison of the surface-EMG fixes. CHANGELOG: fibre dispersion, the literature preset, the published end-plate spreads as iEMG default (nmj_coefficients restores the old ones) and the fibre-length docstring. Coefficient wording corrected to 1.5-2x; preset docstring drops an unverified delay comparison.
gen_arm.py takes _fat_thickness__mm and _skin_conductivity__S_per_m; existing arms still reproduce bit for bit. Eight FDI arms with skin 1.2 mm (Li 2018), fat 0.3/1.3/2.0 mm and skin conductivity 0.15 or 1 S/m, undispersed and with all effects.
The default comes from Farina et al. 2004, who took it from the three-layer model of Roeleveld et al. 1997 (also Blok et al. 2002). It was chosen so the model matched how biceps MUAPs fall off with depth and around the arm; the authors note it is not known whether skin conducts that well. Measured composite skin is ~4.6e-4 to 4.9e-4 S/m (Gabriel 1996 via Lowery 2002 and Botelho 2019). Replaces the earlier, unsourced 'gel bypasses the stratum corneum' rationale. Value unchanged.
Skin 1.29 -> 1.2 mm, the hand-dorsum median of healthy controls in Li et al. 2018 (range 0.79-1.50 mm). Fat 0.3 -> 1.3 mm so that the total depth to muscle is 2.5 mm, the low end of the FDI subcutaneous depth measured by Nishikawa et al. 2024 (group means 2.54-3.30 mm); the fat value is derived from that total, not measured. Default-Muscle surface amplitudes fall by about 40 %. Skin conductivity stays 1 S/m pending further work.
On the default SurfaceEMG route any exception while simulating one motor unit was logged and replaced by an all-zero MUAP, which passes for a valid, silent unit. Seen with sampling_points_in_theta_domain >= 256, where Bessel under/overflow leaves one conductor system singular. The exception is now re-raised with the unit and array index noted; the ExecutionPolicy route already raised.
_fiber_diameter__um sets SurfaceEMG(fiber_diameter__um); the new arms repeat the literature-layer arms at the Hwang et al. 2013 FDI mean.
…2 rates _rate_model selects the De Luca & Contessa (2012) exponential, with parameters anchored to De Luca et al. (1982) Table 1 FDI rates (8.9 Hz recruitment, 25.3 Hz peak at 40 % MVC, 41.4 Hz at 80 %). The fit is a modelling assumption: Table 1 gives per-unit peak rates of decomposed units, not the population mean.
…dent solutions An independent Fourier-Bessel solver for an insulated cylinder with anisotropic muscle, fat and skin (checked against the homogeneous closed form) is compared with MyoGen at FDI layers, 0.5 and 2 mm depth, skin 1 and 0.15 S/m (1 %). A sealed-end cable source (end-plate generation, tendon extinction) convolved with that solver's point-source potential matches the bipolar finite-fibre MUAP at 5 and 10 mm IED (1 %).
_muscle_radius__mm and _fiber_density__per_mm2 set the FDI cylinder; _spike_file replaces the generated trains with NEURON spike trains for the same threshold array (checked on load). Arms: r 5.3 mm, 267 fibres/mm² (~40 k fibres), fat 2.3 mm, NEURON trains.
…rations Spike files from fdi_spikes_force.py (calib/fdi-dd-deluca1982): VL-optimised drive on the FDI pool, and FDI calibrations with recruitment ending at 70 % and 52 % excitation, at 25 % MVC (force model) and 25 % excitation; 267 and 350 fibres/mm², skin 1 and 0.15 S/m. Spike-file paths point to the session scratchpad.
gen_arm.py: _delsys5 records a monopolar 3x3 grid (2.5 mm pitch, 0.25 mm radius pins) per (z, rotation, theta) so the Delsys channels (corner minus centre, adjacent corners) are formed afterwards; _fs__Hz sets the output rate (10 240 Hz for the 20 Hz-2 kHz band of Botelho et al. 2019).
SurfaceEMG(iap_lambda__per_mm=1.0) sets Vm = 96 (lz)^3 e^{-lz} - 90 mV.
The default is the published 1 / mm and leaves outputs unchanged. Values
below 1 stretch the IAP and lower the MUAP's spectral content.
On the FDI defaults, 0.7 matched the median frequency and adjacent-pin RMS
of the Botelho et al. (2019) Delsys recordings at 25 % MVC (one subject);
1 / mm gives MDF ~200 vs ~165 Hz measured. Any value other than 1 is a
fitted modelling assumption; the default is not changed.
Passed through the default and execution= routes (and the disk-cache key);
use_unified=True refuses lambda != 1 because its source would ignore it.
Tests: single fibre against the exact homogeneous cylinder at lambda 0.7
and 0.5, parameter passed to every fibre on both routes, invalid values
rejected. Full suite: 389 passed, 2 skipped.
Repeats Roeleveld et al. (1997) in MyoGen's cylinder (r 35 mm, muscle 32 mm, fat 1, skin 2) with point and distributed motor-unit territories, monopolar contacts 0/12/24/36 mm from the end-plate and around the arm. Decision rule fixed in the README before the runs. Convergence: M 128 is within 0.2 % of 224; M 64 is not (up to 53 % at position D). M >= 256 returns all-zero MUAPs (side finding, not fixed). Replication of their point-unit depth ratios fails the +-15 % limits (e.g. model III 15/5 mm 0.071 vs 0.09), so the fall-off arms are not interpreted. Uses the ExecutionPolicy route (identical output, ~0.12 s per fibre).
…catter replicate.py smooths a single-fibre profile (1 mm rows along z, gen.py z_rows_mm) with a Gaussian of SD 3 mm, which in the z-unbounded cylinder is the expected MUAP of a point unit whose end-plates scatter by that SD. All five replication ratios pass the frozen +-15 % limits (worst -10 %); with SD 0 the first run is reproduced exactly. Production code unchanged. The fall-off arms lack the scatter and still need a rerun.
…sive gen.py places one 1 mm-row column per arc position (arcs_mm); run_scatter_arms.sh runs skin 1 / 0.15 / 4.88e-4 S/m x territory radius 0-12.5 mm; falloff.py smooths with the 3 mm end-plate Gaussian and applies the frozen rule. No arm meets depth 1.8 +- 0.3 and lateral <= 1.3: territories steepen the lateral fall-off; 1 S/m is closest (lateral 1.7-2.0, depth 2.5-2.7) but still too steep; measured skin values are further off. Production code unchanged.
Only ~100 S/m reproduces the measured lateral and depth spread; 1 S/m is the least-bad physical-range value and measured skin values are 3-6x further off.
MyoGen cited Akhmadeev, Yu, Le Carpentier, Aoustin & Farina (2020, IEEE TBME 67:2005) as "Konstantin et al."; Konstantin is the first author's given name. Docstrings, comments, the recruitment-threshold example and CHANGELOG corrected. The konstantin threshold mode keeps its name.
…nductor leaves float64 simulate_fiber_hybrid allocated B(k_z) as real, dropping its imaginary part, so use_unified=True bipolar_longitudinal channels were ~1e-7 of their true value (the #29 fix covered the other routes only). The bipolar-equals-difference tests now also run the hybrid route. At sampling_points_in_theta_domain >= 256 the scaled Bessel functions leave the float64 range (I_n underflows, K_n overflows). The volume conductor zeroed those entries, giving a singular system on the CPU and NaN MUAPs on the GPU. It now raises a ValueError naming the setting; NaN in the bone rows that the no-bone solve drops is still ignored.
The Roeleveld 1997 benchmark was inconclusive and is not evidence for 1 S/m. State the actual reason (a measured 0.15 S/m puts FDI plateau RMS 3-4x above published ranges) in the CHANGELOG and Muscle docstring. Move the upgrade notes to a Breaking changes section and record the use_unified bipolar and angular-sampling fixes.
arms.json pointed at spike files in a local temporary folder. The five files these arms use (120 units each, ~2 MB) now live under benchmarks/semg_dispersion/spikes/ and _spike_file paths are relative to the benchmark directory.
The calibration argument for 1 S/m compared against a 'measured 0.15 S/m' without a source, next to the 4.88e-4 S/m composite skin also called measured. 0.15 S/m is the IT'IS tissue database's low-frequency skin value (0.148 +/- 0.042 S/m, 7 studies, 0.057-0.25).
This was referenced Nov 7, 2026
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Corrects the surface-EMG pipeline against Farina, Mesin, Martina & Merletti (2004) and sets literature-based FDI layer defaults. First of four stacked PRs: this → #43 (ported fixes) → #44 (tooling) → #45 (fibre model); together they replace #41.
Closes #29, closes #32. Part of #30 (skin conductivity and skin/fat thickness; the fibre-diameter part is in #45).
Thanks to @joaohbbittar, whose feedback to us led to the problems fixed in #29.
Fixes (#29)
bipolar_longitudinalchannels were ~10⁻⁷ of their value; rotated arrays distorted.z/2(2× too long).dzandFs/Nfactors (amplitude ∝ 1/CV) and the ad-hocA = 96/10source, replaced by the core-conductor line currentσᵢ·π·(d/2)²·∂²Vm/∂z². A single fibre matches an exact homogeneous-cylinder solution within 2 %.use_unified=Trueroute: B(k_z) is complex too, so its bipolar channels equal the monopolar difference.sampling_points_in_theta_domain≥ 256 raises aValueErrornaming the setting instead of giving a singular system (CPU) or NaN MUAPs (GPU).Defaults (#30, breaking)
benchmarks/semg_dispersion/). The fall-off check inbenchmarks/semg_roeleveld1997/was inconclusive and does not support 1 S/m as a skin value.End-plates (#32, breaking for iEMG)
nmj_coefficients="reference_code"reproduces the old ones); shared nerve-model module; citation corrected from "Konstantin et al. 2020".Also
SurfaceFiberDispersion; off by default, bit-identical when off) andSurfaceEMG(iap_lambda__per_mm=1.0)(default unchanged). These are tracked in Surface fibre dispersion andMusclefibres in surface EMG #36 and Rosenfalck IAP length scale option and the FDI fit #37, which close with Opt-in fibre types by motor unit, surface EMG from Muscle fibres, FDI fibre diameter #45.Tests. Full suite at
5619927a0: 394 passed, 2 skipped. New since the review:use_unifiedin the bipolar = difference-of-monopolars test, and the float64-range errors on the CPU, hybrid and GPU solves. Earlier: layered conductor and finite fibre against independent solutions; equal-conductivity layers reduce to the homogeneous solution.Upgrade notes. sEMG outputs from 0.5.0–0.13.0 change in amplitude (one to two orders of magnitude), duration and spectrum and must be regenerated; see
CHANGELOG.md(Unreleased).