diff --git a/CHANGELOG.md b/CHANGELOG.md
index ea86361d..98ffa7d0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,7 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
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.
-It also ports three fixes that never reached `main` (Gfluctdv noise, fibre conduction velocities, large fibre layouts) and corrects the units of the intramuscular SFAP.
+It also ports three fixes that never reached `main` (Gfluctdv noise, fibre conduction velocities, large fibre layouts), corrects the units of the intramuscular SFAP and adds the multilevel descending-drive optimiser used to calibrate an FDI pool.
### Breaking changes
@@ -62,6 +62,17 @@ This release changes default outputs; regenerate anything calibrated against ear
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).
+- `ForceModel(force_frequency_shape=(a, c))`: shape of the force-frequency saturation S(r) = 1 − exp(−a r^c) in the Fuglevand gain. The default (2, 3) is Fuglevand et al. 1993 and leaves outputs unchanged. `force_utils.HUMAN_UNIT_FORCE_FREQUENCY` (0.93, 1.29) is fitted to human single motor units stimulated through their axons (Macefield et al. 1996, Fuglevand et al. 1999, McNulty et al. 2000): units keep gaining force up to ~30–50 Hz, where Fuglevand's curve saturates by ~20 Hz for slow units. It also gives the measured twitch:tetanus ratio (0.23; Fuglevand 0.11).
+- **Multilevel descending-drive optimiser.** `examples/02_finetune/06_optimize_dd_multilevel.py` fits the descending-drive network (and the Gfluctdv noise scale) to rate-coding targets over a trapezoidal ramp, sequentially or with parallel Optuna workers. Supporting modules: `myogen.utils.fr_optimization_targets` (Fuglevand, De Luca and experimental targets), `myogen.utils.ramp_analysis` and `myogen.utils.parallel_budget` (worker count from free RAM and CPUs).
+ - Every motor neuron receives the same number of descending-drive inputs, and each input's weight is scaled by the neuron's dendrite-to-soma current transfer, so one input delivers the same current to a clamped soma whatever the neuron's size and recruitment follows rheobase. With random input counts and one weight, large neurons received ~2× less current per input and the largest units were recruited late and out of order. Both are modelling assumptions (equal effective synaptic current per input, the equal-share case of Heckman & Binder 1991).
+ - `--last-recruit-mvc PCT`: each trial also runs an MVC ramp (0.5 s rest, excitation 0 → 1 at 0.1/s, 3 s hold; force from `ForceModel`, MVC = mean of the last 2 s) and penalises the distance between the force (% MVC) at which the last unit is recruited and PCT. The calibration ramp stops at 30 % MVC, so this ceiling was previously unconstrained. Off by default. `07_fdi_spikes_force.py` uses the same MVC ramp.
+ - `--ceiling-rank K` scores the unit with K units recruited after it instead of the last unit (less sensitive to one unit's noise). `--force-frequency human` uses the fitted force-frequency curve below for the MVC ramp. `--high-force-rates` also scores the MVC-ramp firing rates at 40 and 80 % MVC against De Luca et al. 1982 (25.3 and 41.4 Hz); the calibration ramp never reaches them. All off by default.
+ - `--fit-weight-exponent` also fits p in [0, 1] for the DD input weights `w·(k_max/k_i)^p`; it is recorded as `network.weight_exponent`. The default p = 1 (equal current per input at rest) is unchanged, and files without the key use it. With p = 1, large units fire as fast as small ones at high force (no onion skin). Lower p gives large units less current per input.
+ - Parameter files record the wiring in a `network` block; files without it are refused.
+- `examples/02_finetune/07_fdi_spikes_force.py` maps excitation to % MVC with `ForceModel` on an MVC ramp and writes spike trains for force-matched trapezoids.
+- `myogen.utils.ramp_analysis`: recruitment and firing rates at % MVC checkpoints of a ramp (`extract_ramp_recruitment`), the force (% MVC) at which each unit is recruited on an MVC ramp (`recruitment_force_pct_mvc`), and unit conduction velocities from recruitment force (`cv_from_recruitment_force`, Del Vecchio et al. 2018; with the residual scatter of their fit by default).
+- `Network.connect_fan_in(source, target, n_inputs, weights__uS)` gives every target neuron exactly `n_inputs` excitatory inputs, drawn without replacement, with one weight per target.
+- `docs/literature.md`: register of the sources behind these defaults and calibrations, each linked to the code that uses it, with verification status (full text, abstract, secondary) and corrections. It is generated by `scripts/build_literature.py`; sourced constants cite it in the code as `[key locator]`, e.g. `[deluca1982 p. 118 Table 1]`, and `tests/test_literature_register.py` fails when the register is stale or the code cites an unknown key.
- `benchmarks/semg_dispersion/`: paired sweeps of these options on a limb-like cylinder and on the default (FDI) muscle.
### Changed
diff --git a/docs/gallery_conf.py b/docs/gallery_conf.py
index e613e74f..b7a915d1 100644
--- a/docs/gallery_conf.py
+++ b/docs/gallery_conf.py
@@ -109,10 +109,12 @@ def reset_neuron(gallery_conf, fname):
# ``__file__`` (undefined under mkdocs-gallery's exec()) to locate a results
# directory and load each other's ``.pkl`` outputs, so they cannot execute
# cleanly in the gallery — they render source-only. The adaptive CPU CLI
- # also requires arguments and user-provided input files, so render its
- # source without executing it. Keep scripts/build_docs.py in sync.
+ # also requires arguments and user-provided input files, and the multilevel
+ # descending-drive optimiser and its spike/force script use ``__file__`` and
+ # long runs (07 takes arguments), so render their source without executing
+ # them. Keep SOURCE_ONLY in scripts/build_docs.py in sync.
# When MKDOCS_GALLERY_PLOT is false, everything is source-only regardless.
- "filename_pattern": r"/(01_basic|02_finetune)/(?!16_adaptive_cpu_execution\.py$)",
+ "filename_pattern": r"/(01_basic|02_finetune)/(?!16_adaptive_cpu_execution\.py$|06_optimize_dd_multilevel\.py$|07_fdi_spikes_force\.py$)",
"ignore_pattern": r"(14_calibrate_noise_from_real|_oscillating_dc_helpers|_optimize_dc_worker|_pic_protocols)\.py",
"within_subsection_order": FileNameSortKey,
"image_scrapers": ("matplotlib",),
diff --git a/docs/literature.md b/docs/literature.md
new file mode 100644
index 00000000..5893829f
--- /dev/null
+++ b/docs/literature.md
@@ -0,0 +1,170 @@
+# Literature register
+
+Sources behind MyoGen's surface-EMG defaults, the fibre-dispersion options and the FDI descending-drive calibration (work of 2026-10-30 to 2026-11-02). Each parameter links to its sources; each source links to the code, tests and notes that use it. Every parameter needs a literature basis or must be labelled a modelling assumption.
+
+**Citing in code:** a sourced value carries an inline comment with the key and a locator, e.g. `HIGH_FORCE_RATE_TARGETS = {40: 25.3, 80: 41.4} # [deluca1982 p. 118 Table 1]`. Give only locators you have read; write "no locator" otherwise. Each key below resolves to the full citation and DOI.
+
+To add a source: add it to the table in `scripts/build_literature.py`, link it from the parameter index, cite the key in the code, and rerun the script (`--check` fails on unknown keys or a stale file). "Cited at" lists are generated from the code.
+
+**Status:** FT full text · FIG read off a figure · ABS abstract only · SEC known only through a citing paper · CODE source code · LINK named only
+**Use:** ADOPTED sets a default or code value (or is its stated basis) · CONTEXT comparator, bound or background · REJECTED considered and not used
+
+## Parameter index
+
+| Parameter or claim | Sources | Code |
+|---|---|---|
+| sEMG volume conductor | [farina2004](#farina2004) | [SurfaceEMG](../myogen/simulator/core/emg/surface/surface_emg.py), [simulate_fiber](../myogen/simulator/core/emg/surface/simulate_fiber.py) |
+| Skin σ = 1 S/m (effective, fitted) | [roeleveld1997](#roeleveld1997), [blok2002](#blok2002), [farinarainoldi1999](#farinarainoldi1999), [farinamerletti2001](#farinamerletti2001), [farina2004](#farina2004), [lowery2002](#lowery2002), [gabriel1996](#gabriel1996), [itis](#itis), [merlettifarina2016](#merlettifarina2016) | [Muscle](../myogen/simulator/core/muscle/muscle.py), [Roeleveld bench](../benchmarks/semg_roeleveld1997/README.md) |
+| Skin thickness 1.2 mm | [li2018](#li2018), [brodar1960](#brodar1960), [krammer2025](#krammer2025) | [Muscle](../myogen/simulator/core/muscle/muscle.py) |
+| Fat thickness 1.3 mm | [nishikawa2024](#nishikawa2024), [storchle2018](#storchle2018) | [Muscle](../myogen/simulator/core/muscle/muscle.py) |
+| Muscle and fat σ | [botelho2019](#botelho2019), [merlettifarina2016](#merlettifarina2016) | [Muscle](../myogen/simulator/core/muscle/muscle.py) |
+| Source strength σi·π(d/2)² | [lowery2002](#lowery2002), [blok2002](#blok2002), [hodgkinnakajima1972](#hodgkinnakajima1972) | [SurfaceEMG](../myogen/simulator/core/emg/surface/surface_emg.py) |
+| IAP shape and λ (iap_lambda__per_mm) | [rosenfalck1969](#rosenfalck1969), [farinamerletti2001](#farinamerletti2001), [botelho2019](#botelho2019), [dimitrova2006](#dimitrova2006), [dimitrov1998](#dimitrov1998), [rodriguezfalces2013](#rodriguezfalces2013), [petersen2019](#petersen2019) | [simulate_fiber](../myogen/simulator/core/emg/surface/simulate_fiber.py), [SurfaceEMG](../myogen/simulator/core/emg/surface/surface_emg.py) |
+| Fibre CV vs diameter | [nandedkar1983](#nandedkar1983), [blijham2006](#blijham2006), [zalewska2012](#zalewska2012) | [Muscle](../myogen/simulator/core/muscle/muscle.py), [SurfaceFiberDispersion](../myogen/simulator/core/emg/surface/fiber_dispersion.py) |
+| Mean fibre CV 4.2 m/s | [nishizono1990](#nishizono1990), [kouyoumdjian2023](#kouyoumdjian2023) | [Muscle](../myogen/simulator/core/muscle/muscle.py) |
+| MU CV vs threshold; within-MU CV spread | [delvecchio2018](#delvecchio2018), [delvecchio2017](#delvecchio2017) | [SurfaceFiberDispersion](../myogen/simulator/core/emg/surface/fiber_dispersion.py) |
+| End-plate spread (two-layer model) | [akhmadeev2020](#akhmadeev2020), [stalbergkarlsson2001](#stalbergkarlsson2001), [piasecki2020](#piasecki2020), [jahanmiri2015](#jahanmiri2015) | [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) |
+| Nerve-branch delays | [akhmadeev2020](#akhmadeev2020), [lateva2010](#lateva2010) | [nerve_branches](../myogen/simulator/core/emg/nerve_branches.py) |
+| Tendon-end band | [keenan2005](#keenan2005), [farinafosci2002](#farinafosci2002), [lateva2010](#lateva2010) | [SurfaceFiberDispersion](../myogen/simulator/core/emg/surface/fiber_dispersion.py) |
+| FDI fibre length and spread | [jacobson1992](#jacobson1992) | [Muscle](../myogen/simulator/core/muscle/muscle.py) |
+| FDI fibre diameter | [polgar1973](#polgar1973), [hwang2013](#hwang2013) | — |
+| FDI MU number, fibres, CSA | [feinstein1955](#feinstein1955), [zong2021](#zong2021), [infantolino2010](#infantolino2010) | — |
+| FDI recruitment ceiling (--last-recruit-mvc) | [deluca1982](#deluca1982), [skarabot2023](#skarabot2023), [milnerbrown1973b](#milnerbrown1973b), [macefield1996](#macefield1996) | [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py) |
+| Rate-coding targets | [deluca1982](#deluca1982), [delucacontessa2012](#delucacontessa2012) | [fr_optimization_targets](../myogen/utils/fr_optimization_targets.py), [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py) |
+| Twitch and force (ForceModel) | [fuglevand1993](#fuglevand1993), [milnerbrown1973a](#milnerbrown1973a), [milnerbrown1973b](#milnerbrown1973b), [fuglevand1999](#fuglevand1999), [mcnulty2000](#mcnulty2000), [macefield1996](#macefield1996) | [ForceModel](../myogen/simulator/core/force/force_model.py), [07_fdi_spikes_force](../examples/02_finetune/07_fdi_spikes_force.py) |
+| Human-unit force-frequency curve (opt-in) | [macefield1996](#macefield1996), [fuglevand1999](#fuglevand1999), [mcnulty2000](#mcnulty2000) | [force_utils](../myogen/simulator/core/force/force_utils.py), [ForceModel](../myogen/simulator/core/force/force_model.py), [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py) |
+| High-force firing rates (--high-force-rates) | [deluca1982](#deluca1982) | [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py) |
+| Gfluctdv noise, ISI CV | [destexhe2001](#destexhe2001), [moritz2005](#moritz2005) | [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py) |
+| DD input weighting, equal fan-in | [heckmanbinder1991](#heckmanbinder1991), [heckmanbinder1993](#heckmanbinder1993) | [DD optimiser](../examples/02_finetune/06_optimize_dd_multilevel.py) |
+| Motoneuron PICs, ΔF, rate patterns (diagnostics 2026-11-02) | [beauchamp2023](#beauchamp2023), [gorassini2002](#gorassini2002), [afsharipour2020](#afsharipour2020), [powersheckman2017](#powersheckman2017), [revillfuglevand2011](#revillfuglevand2011), [revillfuglevand2017](#revillfuglevand2017), [senneff2019](#senneff2019), [deluca1982](#deluca1982) | [AlphaMN__Pool](../myogen/simulator/neuron/populations/motor_neurons.py), [alpha_mn_default.yaml](../myogen/config/alpha_mn_default.yaml) |
+| FDI sEMG comparators (RMS, MDF, MUAP) | [botelho2019](#botelho2019), [mcmanus2015](#mcmanus2015), [mcmanus2019](#mcmanus2019), [milnerbrown1975](#milnerbrown1975), [hu2013](#hu2013), [zong2021](#zong2021) | — |
+
+## A. Volume conductor and recording (tissue layers, conductivities, source strength, IAP)
+
+| Source | Status | Use | What it supports | Used in |
+|---|---|---|---|---|
+| **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)
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), [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) |
+| **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) |
+| **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) |
+| **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) |
+| **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) |
+| **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. | — |
+| **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)
Cited at: [bioelectric.py:12](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L12), [bioelectric.py:232](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L232) |
+| **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. | — |
+| **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) |
+| **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) |
+| **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) |
+| **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'. | — |
+| **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) |
+| **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)
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), [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) |
+| **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). | — |
+| **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. | — |
+| **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. | — |
+| **petersen2019**: Petersen E, Rostalski P (2019). A comprehensive mathematical model of motor unit pool organization, surface electromyography, and force generation. Front Physiol 10:176. [doi:10.3389/fphys.2019.00176](https://doi.org/10.3389/fphys.2019.00176) · [PMC6418040](https://pmc.ncbi.nlm.nih.gov/articles/PMC6418040/) | FT | CONTEXT | Farina-model implementation with λ = 1. | — |
+| **andreassenrosenfalck1980**: Andreassen & Rosenfalck (1980), muscle conductivities for the intramuscular kernel (details from code comment, unverified). | LINK | ADOPTED | Stated source of the radial and longitudinal muscle conductivities of the intramuscular point-electrode kernel in bioelectric.py (same defaults in fiber_simulation.py). Until 0.13.0 the defaults were 63 / 330 labelled S/m (1000x tissue values); since 2026-11-02 they are 0.063 / 0.33 S/m, the values the code comment attributes to this paper. Not in Zotero. Cited in code; not checked. | —
Cited at: [bioelectric.py:121](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L121), [bioelectric.py:123](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L123), [test_iemg_sfap_units.py:7](../tests/test_iemg_sfap_units.py#L7) |
+| **huigen2002**: Huigen E, Peper A, Grimbergen CA (2002). Investigation into the origin of the noise of surface electrodes. Med Biol Eng Comput 40:332–338. [doi:10.1007/BF02344216](https://doi.org/10.1007/BF02344216) | LINK | ADOPTED | With gondran1996, the stated basis of the 1/f^α baseline-drift form with α in 1.5–2.0 (default α 1.75, the midpoint) in IntramuscularEMG.add_noise and emg_noise.py. Not in Zotero. Cited in code; not checked. | —
Cited at: [intramuscular_emg.py:829](../myogen/simulator/core/emg/intramuscular/intramuscular_emg.py#L829), [emg_noise.py:340](../myogen/utils/emg_noise.py#L340), [emg_noise.py:448](../myogen/utils/emg_noise.py#L448), [emg_noise.py:557](../myogen/utils/emg_noise.py#L557), [emg_noise.py:679](../myogen/utils/emg_noise.py#L679), [emg_noise.py:1126](../myogen/utils/emg_noise.py#L1126) |
+| **gondran1996**: Gondran C, Siebert E, Yacoub S, Novakov E (1996). Noise of surface bio-potential electrodes based on NASICON ceramic and Ag–AgCl. Med Biol Eng Comput 34:460–466. [doi:10.1007/BF02523851](https://doi.org/10.1007/BF02523851) | LINK | ADOPTED | With huigen2002, the stated basis of the 1/f^α baseline-drift form with α in 1.5–2.0 (default α 1.75) in IntramuscularEMG.add_noise and emg_noise.py. Not in Zotero. Cited in code; not checked. | —
Cited at: [intramuscular_emg.py:830](../myogen/simulator/core/emg/intramuscular/intramuscular_emg.py#L830), [emg_noise.py:341](../myogen/utils/emg_noise.py#L341), [emg_noise.py:449](../myogen/utils/emg_noise.py#L449), [emg_noise.py:557](../myogen/utils/emg_noise.py#L557), [emg_noise.py:679](../myogen/utils/emg_noise.py#L679), [emg_noise.py:1126](../myogen/utils/emg_noise.py#L1126) |
+| **deluca2010**: De Luca CJ, Gilmore LD, Kuznetsov M, Roy SH (2010). Filtering the surface EMG signal: Movement artifact and baseline noise contamination. J Biomech 43(8):1573–1579. [doi:10.1016/j.jbiomech.2010.01.027](https://doi.org/10.1016/j.jbiomech.2010.01.027) · [PMID 20206934](https://pubmed.ncbi.nlm.nih.gov/20206934/) | LINK | CONTEXT | Cited in IntramuscularEMG.add_noise for broadband movement artifacts (0–20 Hz), which the baseline-drift model leaves out. Citation from Zotero (no full text attached). Cited in code; not checked. | —
Cited at: [intramuscular_emg.py:836](../myogen/simulator/core/emg/intramuscular/intramuscular_emg.py#L836), [emg_noise.py:345](../myogen/utils/emg_noise.py#L345), [emg_noise.py:456](../myogen/utils/emg_noise.py#L456) |
+
+## B. Muscle fibres and innervation (conduction velocity, end-plates, tendon ends, branch delays)
+
+| Source | Status | Use | What it supports | Used in |
+|---|---|---|---|---|
+| **nandedkar1983**: Nandedkar SD, Stålberg E (1983). Simulation of single muscle fibre action potentials. Med Biol Eng Comput 21:158–165. [doi:10.1007/BF02441531](https://doi.org/10.1007/BF02441531) | FT | ADOPTED | Fibre CV = 0.95 + 0.05·d[µm] m/s (fixes main's 2.0 m/s for every fibre); CV–diameter line for diameter-scaled sources. Read 2026-11-02: a modelling construct, not a measured regression — Stålberg 1966 biceps diameters 25–85 µm and CV 3.7 ± 0.71 m/s, with 25 / 55 / 85 µm pinned to the 2-SD band 2.2 / 3.7 / 5.2 m/s; linear, square and square-root forms gave similar amplitudes (p. 163, Fig. 6). Also eq. 1 of delvecchio2018 (p. 7). At FDI values it is inconsistent: 57.9 µm (polgar1973) → 3.85 m/s, 4.2 m/s (nishizono1990) → 65 µm. | [Muscle](../myogen/simulator/core/muscle/muscle.py), [SurfaceFiberDispersion](../myogen/simulator/core/emg/surface/fiber_dispersion.py), [test_muscle_fibre_properties](../tests/test_muscle_fibre_properties.py), [test_surface_dispersion](../tests/test_surface_dispersion.py), [CHANGELOG](../CHANGELOG.md)
Cited at: [fiber_simulation.py:64](../myogen/simulator/core/emg/fiber_simulation.py#L64), [bioelectric.py:14](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L14), [bioelectric.py:110](../myogen/simulator/core/emg/intramuscular/bioelectric.py#L110) |
+| **nishizono1990**: Nishizono H, Fujimoto T, Ohtake H, Miyashita M (1990). Muscle fiber conduction velocity and contractile properties estimated from surface electrode arrays. Electroencephalogr Clin Neurophysiol 75:75–81. [doi:10.1016/0013-4694(90)90154-C](https://doi.org/10.1016/0013-4694(90)90154-C) · [PMID 1688776](https://pubmed.ncbi.nlm.nih.gov/1688776/) | FT | ADOPTED | FDI motor-unit CV by spike-triggered surface arrays (5 mm IED): 50 units from 4 men, 3.2–5.0 m/s, mean 4.2 ± 0.5 (p. 77) — MyoGen's default mean. Units recruited at 120–930 g, i.e. 5–42 % of the 2.2 kg mean MVC (p. 78); skin temperature held at 31–32 °C (p. 77); CV correlated with threshold force, r = 0.52 (p. 79). Botelho 2019 takes its CV range from this paper. | [Muscle](../myogen/simulator/core/muscle/muscle.py) |
+| **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. | — |
+| **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). | — |
+| **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. | — |
+| **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)
Cited at: [ramp_analysis.py:144](../myogen/utils/ramp_analysis.py#L144), [ramp_analysis.py:152](../myogen/utils/ramp_analysis.py#L152), [test_ramp_recruitment_force.py:43](../tests/test_ramp_recruitment_force.py#L43) |
+| **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. | — |
+| **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)
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), [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) |
+| **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). | — |
+| **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). | — |
+| **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. | — |
+| **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) |
+| **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. | — |
+| **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. | — |
+| **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. | — |
+| **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). | — |
+
+## C. Muscle anatomy (FDI dimensions, layers, motor-unit counts)
+
+| Source | Status | Use | What it supports | Used in |
+|---|---|---|---|---|
+| **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) |
+| **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) |
+| **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) |
+| **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) |
+| **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) |
+| **krammer2025**: Krammer et al. (2025). Diagnostics 15:1600. | FT | CONTEXT | Hand-dorsum dermis 0.52 ± 0.16 mm. | — |
+| **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². | — |
+| **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. | — |
+| **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. | — |
+| **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) |
+| **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. | — |
+| **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). | — |
+| **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). | — |
+| **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. | — |
+
+## D. Recruitment and rate coding (thresholds, firing rates, recruitment ceiling)
+
+| Source | Status | Use | What it supports | Used in |
+|---|---|---|---|---|
+| **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)
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) |
+| **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)
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) |
+| **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)
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) |
+| **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. | —
Cited at: [06_optimize_dd_multilevel.py:61](../examples/02_finetune/06_optimize_dd_multilevel.py#L61) |
+| **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. | —
Cited at: [06_optimize_dd_multilevel.py:61](../examples/02_finetune/06_optimize_dd_multilevel.py#L61) |
+| **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)
+
+| Source | Status | Use | What it supports | Used in |
+|---|---|---|---|---|
+| **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) |
+| **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)
Cited at: [06_optimize_dd_multilevel.py:195](../examples/02_finetune/06_optimize_dd_multilevel.py#L195) |
+| **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. | — |
+| **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) |
+| **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. | — |
+| **afsharipour2020**: Afsharipour B, Manzur N, Duchcherer J, Fenrich KF, Thompson CK, Negro F, et al. (2020). Estimation of self-sustained activity produced by persistent inward currents using firing rate profiles of multiple motor units in humans. J Neurophysiol 124:63–85. [doi:10.1152/jn.00194.2020](https://doi.org/10.1152/jn.00194.2020) · [PMID 32459555](https://pubmed.ncbi.nlm.nih.gov/32459555/) | LINK | CONTEXT | ΔF pairing criteria (control/test timing, rate-rate correlation); check 2 uses the common criteria from memory, not checked against this paper. | — |
+| **powersheckman2017**: Powers RK, Heckman CJ (2017). Synaptic control of the shape of the motoneuron pool input-output function. J Neurophysiol 117:1171–1184. [doi:10.1152/jn.00850.2016](https://doi.org/10.1152/jn.00850.2016) · [PMID 28053245](https://pubmed.ncbi.nlm.nih.gov/28053245/) · [PMC5340877](https://pmc.ncbi.nlm.nih.gov/articles/PMC5340877/) · [link](https://modeldb.science/239582) | ABS | CONTEXT | Cat MG model pool: PICs raise gain and make pool rate a nonlinear function of excitation; neuromodulation and balanced inhibition increase the nonlinearity, push-pull inhibition reduces it; later-recruited units give higher force. Origin of MyoGen's Powers2017 cell model: its mechanisms na3rp, naps, kdrRL, mAHP, gh, L_Ca_inact and Gfluctdv are those of ModelDB 239582, this paper's model (2026-11-02 diff: whitespace, a seed RANGE and comment lines). Cell parameter values not compared. | [AlphaMN__Pool](../myogen/simulator/neuron/populations/motor_neurons.py), [alpha_mn_default.yaml](../myogen/config/alpha_mn_default.yaml)
Cited at: [alpha_mn_default.yaml:10](../myogen/config/alpha_mn_default.yaml#L10), [alpha_mn_default.yaml:89](../myogen/config/alpha_mn_default.yaml#L89), [cells.py:759](../myogen/simulator/neuron/cells.py#L759), [cells.py:878](../myogen/simulator/neuron/cells.py#L878), [motor_neurons.py:27](../myogen/simulator/neuron/populations/motor_neurons.py#L27), [motor_neurons.py:57](../myogen/simulator/neuron/populations/motor_neurons.py#L57) |
+| **revillfuglevand2011**: Revill AL, Fuglevand AJ (2011). Effects of persistent inward currents, accommodation, and adaptation on motor unit behavior: a simulation study. J Neurophysiol 106:1467–1479. [doi:10.1152/jn.00419.2011](https://doi.org/10.1152/jn.00419.2011) · [PMID 21697447](https://pubmed.ncbi.nlm.nih.gov/21697447/) | LINK | CONTEXT | PIC simulation of rate saturation and onion skin. Not read. | — |
+| **revillfuglevand2017**: Revill AL, Fuglevand AJ (2017). Inhibition linearizes firing rate responses in human motor units: implications for the role of persistent inward currents. J Physiol 595:179–191. [doi:10.1113/JP272823](https://doi.org/10.1113/JP272823) · [PMID 27470946](https://pubmed.ncbi.nlm.nih.gov/27470946/) | ABS | CONTEXT | Human MU rates rise steeply at recruitment then level off at modest rates while force keeps rising; reflex inhibition removes the steep rise (attributed to PICs). | — |
+| **senneff2019**: Senneff S, McManus L, Lowery MM (2019). Investigating the effect of persistent inward currents on motor unit firing rates and beta-band coherence in a model of the first dorsal interosseous muscle. Conf Proc IEEE EMBC 2019:2293–2296. [doi:10.1109/EMBC.2019.8857534](https://doi.org/10.1109/EMBC.2019.8857534) · [PMID 31946358](https://pubmed.ncbi.nlm.nih.gov/31946358/) | ABS | CONTEXT | Prior FDI motoneuron-pool model with PICs at 20–40 % MVC: PICs shift mean rates and reduce beta-band coherence. Prior art for the spinal checks. | — |
+| **powersbinder2001**: Powers RK, Binder MD (2001). Input-output functions of mammalian motoneurons. Rev Physiol Biochem Pharmacol 143:137–263. [doi:10.1007/BFb0115594](https://doi.org/10.1007/BFb0115594) | LINK | CONTEXT | Review of motoneuron input-output properties; cited in the AlphaMN cell docstring. Not read here. | —
Cited at: [cells.py:758](../myogen/simulator/neuron/cells.py#L758) |
+| **bui2003**: Bui et al. (2003). Morphological and electrophysiological characterization of Ia inhibitory interneurons in cat spinal cord (details from code comment, unverified). | LINK | ADOPTED | Membrane area (81,390 µm²) and soma geometry of the INgII interneuron model in cells.py. Cited in code; not checked. | —
Cited at: [cells.py:235](../myogen/simulator/neuron/cells.py#L235), [cells.py:256](../myogen/simulator/neuron/cells.py#L256), [cells.py:264](../myogen/simulator/neuron/cells.py#L264), [cells.py:281](../myogen/simulator/neuron/cells.py#L281), [cells.py:285](../myogen/simulator/neuron/cells.py#L285) |
+| **sanches2018**: Sanches et al. (2018), synaptic parameter measurements (details from code comment, unverified). | LINK | ADOPTED | Default Exp2Syn synapse parameters (tau1 0.2 ms, tau2 0.3 ms, weight 0.6) in cells.py. Cited in code; not checked. | —
Cited at: [cells.py:200](../myogen/simulator/neuron/cells.py#L200), [cells.py:206](../myogen/simulator/neuron/cells.py#L206) |
+| **segev1990**: Segev et al. (1990), alpha-function synaptic modelling (details from code comment, unverified). | LINK | CONTEXT | Listed as the reference for the synapse model form in cells.py (the code uses a double-exponential Exp2Syn). Cited in code; not checked. | —
Cited at: [cells.py:207](../myogen/simulator/neuron/cells.py#L207) |
+
+## F. Muscle force (twitch, contraction time, force-frequency)
+
+| Source | Status | Use | What it supports | Used in |
+|---|---|---|---|---|
+| **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. | — |
+| **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) |
+| **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)
Cited at: [force_utils.py:6](../myogen/simulator/core/force/force_utils.py#L6) |
+| **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)
Cited at: [force_utils.py:6](../myogen/simulator/core/force/force_utils.py#L6) |
+| **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)
Cited at: [force_utils.py:6](../myogen/simulator/core/force/force_utils.py#L6) |
+| **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)
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) |
+
+## G. Experimental EMG benchmarks (FDI amplitude and spectra)
+
+| Source | Status | Use | What it supports | Used in |
+|---|---|---|---|---|
+| **milnerbrown1975**: Milner-Brown HS, Stein RB (1975). The relation between the surface electromyogram and muscular force. J Physiol 246:549–569. [doi:10.1113/jphysiol.1975.sp010904](https://doi.org/10.1113/jphysiol.1975.sp010904) | FT | CONTEXT | FDI surface MUAP 0.14–0.40 mV (belly–tendon); duration 1.6–2.2 ms. | — |
+| **hu2013**: Hu X, Rymer WZ, Suresh NL (2013). J Neurophysiol 110:1205; companion J NeuroEng Rehabil (2013). [PMC4073930](https://pmc.ncbi.nlm.nih.gov/articles/PMC4073930/) · [PMC3849018](https://pmc.ncbi.nlm.nih.gov/articles/PMC3849018/) | FT + FIG | CONTEXT | FDI decomposition; Delsys 5-pin geometry; MUAP amplitude linear in threshold. | — |
+| **mcmanus2015**: McManus L et al. (2015). J Neurophysiol 113:3186. [doi:10.1152/jn.00146.2015](https://doi.org/10.1152/jn.00146.2015) | FT | CONTEXT | FDI surface MDF 166.9 Hz. | — |
+| **mcmanus2019**: McManus L, Flood MW, Lowery MM (2019). J Neurophysiol. [doi:10.1152/jn.00228.2019](https://doi.org/10.1152/jn.00228.2019) · [PMC6766730](https://pmc.ncbi.nlm.nih.gov/articles/PMC6766730/) | FT + FIG | CONTEXT | MDF ~120–150 Hz (400 Hz low-pass). 'MUAP amplitude' rejected as a target (metric undefined). | — |
+
+## Corrections and open verification
+
+1. **Skin σ = 1 S/m rationale.** First explained as gel electrodes bypassing the stratum corneum (Farina 2004). Unsourced and withdrawn: the value is Roeleveld et al. 1997's fitted effective skin, inherited by Blok 2002, Farina & Rainoldi 1999, Farina & Merletti 2001 and Farina 2004.
+2. **Roeleveld 1997 DOI** is 10.1016/S1050-6411(97)00009-6 (an earlier (96)00034-8 was wrong).
+3. **'Konstantin et al. 2020'** in MyoGen docstrings (muscle.py, intramuscular_emg.py, nerve_branches.py) is Akhmadeev et al. 2020. Fixed on `fix/semg-amplitude-scale` (028a14362); the `konstantin` threshold mode keeps its name.
+4. **Botelho 2019 simulated amplitude**: median ~0.16 mV; 0.10 mV was one sample trace. Their 24.8–40.6 µm diameters are Hwang 2013's FDI mean ± SD.
+5. **Feinstein 1955** authors are Feinstein, Lindegård, Nyman, Wohlfart; ~40,000 fibres (not 40,500).
+6. **Piasecki 2020 and Stålberg & Karlsson 2001** are one citation chain, not two anchors.
+7. **Jacobson 1992's 2.8 mm** is probably between-specimen spread; MyoGen uses it as a uniform half-width (docstring corrected).
+8. **Unverified, quoted from memory:** Moritz 2005 (ISI CV), Kukulka & Clamann 1981. **Never read in full:** Nandedkar & Stålberg 1983 (adopted CV line), Jacobson 1992 Table 1, Heckman & Binder 1991.
+9. **McManus 2019 'MUAP amplitude'**: metric undefined, not usable as a target.
+10. **ISI CV target 15–22 %** in the optimiser cites enokaduchateau2017 and dideriksen2012; the range was not checked against either paper.
diff --git a/examples/02_finetune/06_optimize_dd_multilevel.py b/examples/02_finetune/06_optimize_dd_multilevel.py
new file mode 100644
index 00000000..835330ee
--- /dev/null
+++ b/examples/02_finetune/06_optimize_dd_multilevel.py
@@ -0,0 +1,1080 @@
+#!/usr/bin/env python3
+"""
+Multilevel Descending-Drive Optimisation on a Trapezoidal Ramp
+==============================================================
+
+Optimizes DD network parameters to match Fuglevand [fuglevand1993] recruitment targets
+across multiple %MVC levels simultaneously. Uses a single trapezoidal
+ramp per trial instead of separate constant-drive simulations.
+
+Usage:
+ cd examples/02_finetune
+ python 06_optimize_dd_multilevel.py --n-mu 200 --n-trials 100
+ python 06_optimize_dd_multilevel.py --n-mu 150 --n-trials 100
+ python 06_optimize_dd_multilevel.py --sweep # runs all 4 MU counts
+"""
+
+import argparse
+import json
+import math
+import os
+import subprocess
+import sys
+import time
+
+os.environ["MPLBACKEND"] = "Agg"
+
+from pathlib import Path
+from types import SimpleNamespace
+
+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
+
+from myogen import set_random_seed
+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
+
+# ─── 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)
+CEILING_RAMP_PER_S = 0.1 # Excitation per second in the MVC ramp (also used by 07_fdi_spikes_force.py)
+MVC_REST_S, MVC_HOLD_S = 0.5, 3.0 # MVC ramp: rest, ramp 0 -> 1 at CEILING_RAMP_PER_S, hold at full excitation
+FORCE_FS_HZ = 2048.0 # ForceModel sampling rate
+FORCE_SMOOTH_S = 0.25 # Moving-average window on the ramp force
+GFLUCTDV_RANGE = (0.05, 0.5) # Optuna search range for Gfluctdv scale (FDI best was 0.243 at the old 0.25 cap)
+ISI_CV_TARGET = (0.15, 0.22) # Literature target: 15-22% at 30 % MVC [enokaduchateau2017] [dideriksen2012] (range not checked against the papers)
+ISI_CV_WEIGHT = 2.0 # Penalty weight for ISI CV outside target
+TIMESTEP_MS = 0.1
+RAMP_MVC_PCT = 30 # Ramp target: 0→30% MVC
+REST_S = 3.0
+RAMP_S = 4.0
+PLATEAU_S = 6.0
+CHECKPOINTS = [5, 10, 20, 25, 30]
+CHECKPOINT_WEIGHTS = {5: 1.0, 10: 3.0, 20: 2.0, 25: 2.0, 30: 1.0}
+# Threshold and target options; set from the command line in main().
+THRESHOLD_MODE = "combined" # RecruitmentThresholds mode for the pool
+RECRUITMENT_RANGE = 30.0 # largest / smallest threshold
+MAX_RECRUITMENT = 1.0 # excitation at which the last unit is recruited (rate-coding targets only)
+LAST_RECRUIT_PCT_MVC: float | None = None # force (% MVC) at which the last unit is recruited; None: no ceiling term
+CEILING_RANK = 0 # score the unit with CEILING_RANK units recruited above it (0: the last unit)
+FORCE_FREQUENCY_SHAPE = FUGLEVAND_FORCE_FREQUENCY # (a, c) of the ForceModel force-frequency saturation
+FORCE_FREQUENCY_SHAPES = {"fuglevand": FUGLEVAND_FORCE_FREQUENCY, "human": HUMAN_UNIT_FORCE_FREQUENCY}
+# Mean FDI peak firing rates (Hz) in 40 and 80 % MVC contractions [deluca1982 p. 118 Table 1].
+HIGH_FORCE_RATE_TARGETS = {40: 25.3, 80: 41.4}
+HIGH_FORCE_RATE_WEIGHT = 1.0 # per level, relative rate error; scored on the MVC ramp
+SCORE_HIGH_FORCE_RATES = False
+DE_LUCA_PARAMS: dict = {} # overrides for de_luca_targets
+RESULTS_TAG = "polygon"
+
+
+def _thresholds(profile, n_mu):
+ """Pool thresholds (normalised to 1) and the excitation thresholds used for the targets.
+
+ ``profile`` is the muscle key; on this branch thresholds come straight from ``RecruitmentThresholds``
+ (the muscle-profile library is not part of main).
+ """
+ thr, _ = RecruitmentThresholds(N=n_mu, recruitment_range__ratio=RECRUITMENT_RANGE, mode=THRESHOLD_MODE)
+ thr = np.asarray(thr, dtype=float)
+ thr = thr / thr.max()
+ return thr, thr * MAX_RECRUITMENT
+# Fuglevand model parameters [fuglevand1993] (back-compat)
+MIN_FR_HZ = 8.0
+MAX_FR_HZ = 17.0
+
+RESULTS_BASE = SCRIPT_DIR.parent.parent / "results" # same root as the other 02_finetune examples
+
+
+def build_ramp_drive(threshold_hz, max_drive_hz, timestep_ms=TIMESTEP_MS):
+ """Build a 0→30% MVC trapezoidal ramp drive signal."""
+ peak_hz = threshold_hz + (RAMP_MVC_PCT / 100.0) * (max_drive_hz - threshold_hz)
+ rest_hz = max(threshold_hz * 0.8, 0)
+
+ ramp_ms = RAMP_S * 1000
+ plat_ms = PLATEAU_S * 1000
+ total_ms = (2 * REST_S + 2 * RAMP_S + PLATEAU_S) * 1000
+ n_steps = int(total_ms / timestep_ms)
+
+ signal = np.full(n_steps, rest_hz)
+ rest1_end = int(REST_S * 1000 / timestep_ms)
+ ramp_up_end = rest1_end + int(ramp_ms / timestep_ms)
+ plat_end = ramp_up_end + int(plat_ms / timestep_ms)
+ ramp_down_end = plat_end + int(ramp_ms / timestep_ms)
+
+ signal[rest1_end:ramp_up_end] = np.linspace(rest_hz, peak_hz, ramp_up_end - rest1_end)
+ signal[ramp_up_end:plat_end] = peak_hz
+ signal[plat_end:ramp_down_end] = np.linspace(peak_hz, rest_hz, ramp_down_end - plat_end)
+
+ rng = np.random.default_rng(42)
+ signal += np.abs(rng.normal(0, peak_hz * 0.01, size=n_steps))
+
+ return signal, total_ms
+
+
+def detect_threshold(thresholds, dd_params, min_fr=5.0):
+ """Binary search for the drive Hz that produces sustained firing."""
+ max_drive = dd_params["dd_drive__Hz"]
+ lo, hi = 0.05, 0.60
+ threshold_hz = hi * max_drive
+ for _ in range(8):
+ mid = (lo + hi) / 2.0
+ drive_hz = mid * max_drive
+ block = run_dd_constant(thresholds, drive_hz, dd_params, duration_ms=3000.0)
+ has_sustained = any(
+ len(st) >= int(min_fr * 3.0)
+ for st in block.segments[0].spiketrains
+ )
+ if has_sustained:
+ threshold_hz = drive_hz
+ hi = mid
+ else:
+ lo = mid
+ 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)}
+
+
+def recruitment_times(block, window_ms=500.0, min_spikes=3):
+ """Recruitment time (s) of each unit: the first spike of the first ``min_spikes`` spikes within ``window_ms``
+ (sustained firing, as in ``extract_ramp_recruitment``). NaN for units never recruited."""
+ out = []
+ for st in block.segments[0].spiketrains:
+ t = np.asarray(st.rescale("s").magnitude)
+ hits = np.nonzero(t[min_spikes - 1:] - t[:len(t) - min_spikes + 1] < window_ms / 1000.0)[0]
+ out.append(float(t[hits[0]]) if len(hits) else np.nan)
+ return np.array(out)
+
+
+def recruitment_excitation(block, excitation, timestep_ms=TIMESTEP_MS, window_ms=500.0, min_spikes=3):
+ """Excitation at each unit's recruitment (see ``recruitment_times``). NaN for units never recruited."""
+ t_rec = recruitment_times(block, window_ms, min_spikes)
+ idx = np.clip(np.round(np.nan_to_num(t_rec) * 1000.0 / timestep_ms).astype(int), 0, len(excitation) - 1)
+ return np.where(np.isnan(t_rec), np.nan, np.asarray(excitation, dtype=float)[idx])
+
+
+def smoothed_force(force, fs):
+ w = max(int(FORCE_SMOOTH_S * fs), 1)
+ return np.convolve(force, np.ones(w) / w, mode="same")
+
+
+def last_recruitment_pct_mvc(t_rec, force, fs, mvc, rank=0):
+ """Force (% MVC, smoothed) at the recruitment of the unit with ``rank`` units recruited after it (0: the last
+ unit), and the number of units never recruited. Unrecruited units count as recruited after it; more of them than
+ ``rank`` push the value past 100 % (100 * (1 + (n_unrecruited - rank) / n))."""
+ t_rec = np.asarray(t_rec, dtype=float)
+ n_unrec = int(np.isnan(t_rec).sum())
+ if n_unrec > rank:
+ return 100.0 * (1.0 + (n_unrec - rank) / len(t_rec)), n_unrec
+ t = np.sort(t_rec[~np.isnan(t_rec)])[::-1][rank - n_unrec]
+ i = min(int(round(t * fs)), len(force) - 1)
+ return float(100.0 * smoothed_force(force, fs)[i] / mvc), n_unrec
+
+
+def rates_at_pct_mvc(trains, force, fs, mvc, levels=tuple(HIGH_FORCE_RATE_TARGETS), window_s=1.0, min_spikes=3):
+ """Mean firing rate (Hz) of the active units (>= min_spikes in window_s) centred where the smoothed force first
+ reaches each level (% MVC). ``trains`` are spike times in s. NaN for a level never reached."""
+ sm = smoothed_force(force, fs)
+ out = {}
+ for lev in levels:
+ hit = np.nonzero(sm >= lev / 100.0 * mvc)[0]
+ if not len(hit):
+ out[lev] = float("nan")
+ continue
+ t0 = hit[0] / fs
+ counts = [np.sum((np.asarray(st) > t0 - window_s / 2) & (np.asarray(st) <= t0 + window_s / 2)) for st in trains]
+ rates = [c / window_s for c in counts if c >= min_spikes]
+ out[lev] = float(np.mean(rates)) if rates else 0.0
+ return out
+
+
+def high_force_rate_penalty(rates):
+ """HIGH_FORCE_RATE_WEIGHT * relative rate error at each HIGH_FORCE_RATE_TARGETS level (1 per missing level)."""
+ err = 0.0
+ for lev, tgt in HIGH_FORCE_RATE_TARGETS.items():
+ r = rates.get(lev, float("nan"))
+ err += 1.0 if not np.isfinite(r) else abs(r - tgt) / tgt
+ return HIGH_FORCE_RATE_WEIGHT * err
+
+
+def ceiling_penalty(pct_last, target_pct=None):
+ """CEILING_WEIGHT * relative miss of the last recruitment force; 0 when no target is set."""
+ if target_pct is None:
+ return 0.0
+ 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,
+ force_frequency_shape=FORCE_FREQUENCY_SHAPE)
+
+
+def simulate_force(thresholds, dd_params, threshold_hz, segments):
+ """Run the drive of ``mvc_drive``; return (excitation, block, force at FORCE_FS_HZ) (see ``make_force_model``)."""
+ e, hz = mvc_drive(segments, threshold_hz, dd_params["dd_drive__Hz"])
+ block = run_dd_ramp(thresholds, hz, dd_params, len(e) * TIMESTEP_MS)
+ return e, block, np.asarray(make_force_model(thresholds).generate_force(block).magnitude)[:, 0]
+
+
+def run_mvc_ramp(thresholds, dd_params, threshold_hz):
+ """MVC ramp: MVC_REST_S rest, excitation 0 -> 1 at CEILING_RAMP_PER_S, MVC_HOLD_S hold at 1.
+ MVC is the mean force over the last 2 s of the hold. Returns a namespace with e, block, force, mvc,
+ t_rec (recruitment times, s), pct_last / n_unrecruited (see ``last_recruitment_pct_mvc``), pct_ranked
+ (the CEILING_RANK unit) and rates (``rates_at_pct_mvc`` at the HIGH_FORCE_RATE_TARGETS levels)."""
+ ramp_s = 1.0 / CEILING_RAMP_PER_S
+ e, block, f = simulate_force(thresholds, dd_params, threshold_hz,
+ [(MVC_REST_S, 0.0, 0.0), (ramp_s, 0.0, 1.0), (MVC_HOLD_S, 1.0, 1.0)])
+ t = np.arange(f.size) / FORCE_FS_HZ
+ t_end = MVC_REST_S + ramp_s + MVC_HOLD_S
+ mvc = float(np.mean(f[(t > t_end - 2.0) & (t < t_end)]))
+ t_rec = recruitment_times(block)
+ pct_last, n_unrec = last_recruitment_pct_mvc(t_rec, f, FORCE_FS_HZ, mvc)
+ pct_ranked, _ = last_recruitment_pct_mvc(t_rec, f, FORCE_FS_HZ, mvc, rank=CEILING_RANK)
+ trains = [np.asarray(st.rescale("s").magnitude) for st in block.segments[0].spiketrains]
+ rates = rates_at_pct_mvc(trains, f, FORCE_FS_HZ, mvc)
+ return SimpleNamespace(e=e, block=block, force=f, mvc=mvc, t_rec=t_rec, pct_last=pct_last,
+ n_unrecruited=n_unrec, pct_ranked=pct_ranked, rates=rates)
+
+
+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
+
+
+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
+
+
+def make_objective(thresholds, fug_targets):
+ """Create the Optuna objective function for a given MU count."""
+
+ def objective(trial):
+ try:
+ dd_neurons = trial.suggest_int("dd_neurons", 50, 500)
+ conn_probability = trial.suggest_float("conn_prob", 0.05, 0.5)
+ dd_drive__Hz = trial.suggest_float("dd_drive", 100.0, 1000.0)
+ gfluctdv_scale = trial.suggest_float(
+ "gfluctdv_scale", *GFLUCTDV_RANGE,
+ )
+ weight_exponent = (trial.suggest_float("weight_exponent", *WEIGHT_EXPONENT_RANGE)
+ if FIT_WEIGHT_EXPONENT else None)
+
+ dd_params = {
+ "dd_neurons": dd_neurons,
+ "conn_probability": conn_probability,
+ "dd_drive__Hz": dd_drive__Hz,
+ "synaptic_weight": SYNAPTIC_WEIGHT,
+ "gamma_shape": GAMMA_SHAPE,
+ "gfluctdv_scale": gfluctdv_scale,
+ "network": network_block(weight_exponent),
+ }
+
+ # Detect threshold for this param set
+ threshold_hz = detect_threshold(thresholds, dd_params)
+
+ # Build and run trapezoidal ramp
+ drive_signal, total_ms = build_ramp_drive(threshold_hz, dd_drive__Hz)
+ block = run_dd_ramp(thresholds, drive_signal, dd_params, total_ms)
+
+ # Extract recruitment curve at checkpoints
+ ramp_data = extract_ramp_recruitment(
+ block, drive_signal, threshold_hz, dd_drive__Hz,
+ checkpoints_pct=CHECKPOINTS,
+ )
+
+ # Score against Fuglevand targets
+ total_error = 0.0
+ for mvc_pct in CHECKPOINTS:
+ sim = ramp_data[mvc_pct]
+ tgt = fug_targets[mvc_pct]
+ w = CHECKPOINT_WEIGHTS.get(mvc_pct, 1.0)
+
+ if tgt["n_active"] == 0:
+ total_error += w * sim["n_active"] * 0.1
+ continue
+
+ n_err = abs(sim["n_active"] - tgt["n_active"]) / max(tgt["n_active"], 1)
+ fr_err = (
+ abs(sim["mean_fr_hz"] - tgt["mean_fr_hz"]) / max(tgt["mean_fr_hz"], 1.0)
+ if sim["n_active"] > 0 else 1.0
+ )
+ total_error += w * (n_err + fr_err)
+
+ # ISI CV measurement at plateau (ARCH-23)
+ plat_start_s = (REST_S + RAMP_S)
+ plat_end_s = plat_start_s + PLATEAU_S
+ isi_cvs = []
+ for st in block.segments[0].spiketrains:
+ times_s = np.array(st.rescale("s").magnitude)
+ plat = times_s[(times_s >= plat_start_s) & (times_s < plat_end_s)]
+ if len(plat) < 5:
+ continue
+ isis = np.diff(plat) * 1000
+ m = float(np.mean(isis))
+ if m > 0:
+ isi_cvs.append(float(np.std(isis)) / m)
+ mean_cv = float(np.mean(isi_cvs)) if isi_cvs else 0.0
+
+ # Penalize ISI CV outside the 15-22% target range
+ cv_lo, cv_hi = ISI_CV_TARGET
+ if mean_cv < cv_lo:
+ total_error += ISI_CV_WEIGHT * (cv_lo - mean_cv) / cv_lo
+ elif mean_cv > cv_hi:
+ total_error += ISI_CV_WEIGHT * (mean_cv - cv_hi) / cv_hi
+
+ # Recruitment ceiling: the ramp above stops at RAMP_MVC_PCT, so the last unit is found on a separate
+ # MVC ramp, scored by the force (% MVC) at its recruitment. Units never recruited push it past 100 %.
+ # The same ramp gives the firing rates at high force (De Luca 1982) [deluca1982], which the ramp above never reaches.
+ pct_last, pct_ranked, n_unrecruited, hf_rates = float("nan"), float("nan"), None, {}
+ if LAST_RECRUIT_PCT_MVC is not None or SCORE_HIGH_FORCE_RATES:
+ mvc_ramp = run_mvc_ramp(thresholds, dd_params, threshold_hz)
+ pct_last, n_unrecruited = mvc_ramp.pct_last, mvc_ramp.n_unrecruited
+ pct_ranked, hf_rates = mvc_ramp.pct_ranked, mvc_ramp.rates
+ total_error += ceiling_penalty(pct_ranked, LAST_RECRUIT_PCT_MVC)
+ if SCORE_HIGH_FORCE_RATES:
+ total_error += high_force_rate_penalty(hf_rates)
+
+ # Store metadata
+ trial.set_user_attr("pct_mvc_last_recruited", pct_last)
+ trial.set_user_attr("pct_mvc_ceiling_ranked", pct_ranked)
+ for lev, r in hf_rates.items():
+ trial.set_user_attr(f"mean_fr_{lev}_mvc_ramp", r)
+ trial.set_user_attr("n_unrecruited", n_unrecruited)
+ trial.set_user_attr("threshold_hz", threshold_hz)
+ trial.set_user_attr("isi_cv_mean", mean_cv)
+ trial.set_user_attr("gfluctdv_scale", gfluctdv_scale)
+ for mvc_pct in CHECKPOINTS:
+ trial.set_user_attr(f"n_active_{mvc_pct}", ramp_data[mvc_pct]["n_active"])
+ trial.set_user_attr(f"mean_fr_{mvc_pct}", ramp_data[mvc_pct]["mean_fr_hz"])
+
+ if trial.number % 5 == 0:
+ active_str = " | ".join(
+ f"{p}%:{ramp_data[p]['n_active']}/{fug_targets[p]['n_active']}"
+ for p in CHECKPOINTS
+ )
+ print(f" Trial {trial.number}: err={total_error:.3f} | {active_str} | "
+ f"drive={dd_drive__Hz:.0f} thr={threshold_hz:.0f} "
+ f"gf={gfluctdv_scale:.3f} cv={mean_cv*100:.1f}% last={pct_last:.0f}%MVC"
+ + (f" p={weight_exponent:.2f}" if weight_exponent is not None else ""))
+
+ return total_error
+
+ except Exception as e:
+ print(f" Trial {trial.number} failed: {e}")
+ return 1000.0
+
+ return objective
+
+
+def optimize_for_mu_count(n_mu, n_trials, profile, target_model="de_luca",
+ fresh=False,
+ target_muscle="vl",
+ sampler_seed=42,
+ save_json=True):
+ """Run the full optimization for one MU count.
+
+ Parameters
+ ----------
+ n_mu : int
+ Motor unit count.
+ n_trials : int
+ Optuna trial count.
+ profile : MuscleProfile
+ Profile to generate thresholds from.
+ target_model : {"fuglevand", "de_luca", "experimental"}
+ Rate coding target model (see fr_optimization_targets).
+ target_muscle : str
+ Muscle name for experimental lookup (ignored by other models).
+ sampler_seed : int
+ TPE sampler seed. Use different values per parallel worker.
+ save_json : bool
+ Whether to save the best-trial JSON. False for parallel workers
+ (the coordinator process saves after all workers finish).
+ """
+ print(f"\n{'='*60}")
+ print(f"Optimizing DD for {n_mu} MUs (density=312/mm2, Lexell anchor)")
+ print(f"Target model: {target_model}")
+ print(f"{'='*60}")
+
+ set_random_seed(42)
+ thresholds, target_thresholds = _thresholds(profile, n_mu)
+
+ # Rate coding targets from selected model
+ if target_model == "fuglevand":
+ fug = compute_fr_targets(
+ target_thresholds, mvc_levels=CHECKPOINTS, model="fuglevand",
+ min_fr_hz=MIN_FR_HZ, max_fr_hz=MAX_FR_HZ,
+ )
+ elif target_model == "experimental":
+ fug = compute_fr_targets(
+ target_thresholds, mvc_levels=CHECKPOINTS, model="experimental",
+ muscle_name=target_muscle,
+ )
+ else: # de_luca (default)
+ fug = compute_fr_targets(
+ target_thresholds, mvc_levels=CHECKPOINTS, model="de_luca", **DE_LUCA_PARAMS,
+ )
+ print(f"Rate coding targets ({target_model}):")
+ for mvc_pct in CHECKPOINTS:
+ t = fug[mvc_pct]
+ print(f" {mvc_pct}% MVC: {t['n_active']} active, {t['mean_fr_hz']:.1f} Hz")
+
+ # Optuna study
+ results_dir = RESULTS_BASE / f"dd_{target_muscle}_{n_mu}mu_{RESULTS_TAG}"
+ results_dir.mkdir(exist_ok=True, parents=True)
+ db_path = results_dir / "optuna_multilevel.db"
+ if fresh and db_path.exists():
+ db_path.unlink()
+ print(f" Removed existing {db_path} (--fresh)")
+ storage = f"sqlite:///{db_path}"
+
+ study = optuna.create_study(
+ direction="minimize",
+ sampler=optuna.samplers.TPESampler(
+ seed=sampler_seed, multivariate=True,
+ warn_independent_sampling=False,
+ ),
+ study_name=f"dd_multilevel_{n_mu}mu",
+ storage=storage,
+ load_if_exists=True,
+ )
+
+ objective_fn = make_objective(thresholds, fug)
+ t0 = time.time()
+ study.optimize(objective_fn, n_trials=n_trials, show_progress_bar=False)
+ elapsed = time.time() - t0
+
+ # Best trial
+ best = study.best_trial
+ print(f"\nBest trial {best.number} (err={best.value:.3f}, {elapsed:.0f}s):")
+ print(f" dd_neurons={best.params['dd_neurons']}, "
+ f"conn_prob={best.params['conn_prob']:.3f}, "
+ f"dd_drive={best.params['dd_drive']:.0f} Hz, "
+ f"gfluctdv={best.params.get('gfluctdv_scale', 0):.3f}, "
+ f"weight exponent={best.params.get('weight_exponent', 1.0):.2f}")
+ print(f" ISI CV: {best.user_attrs.get('isi_cv_mean', 0)*100:.1f}% "
+ f"(target: {ISI_CV_TARGET[0]*100:.0f}-{ISI_CV_TARGET[1]*100:.0f}%)")
+ print(f" Last unit recruited at {best.user_attrs.get('pct_mvc_last_recruited', float('nan')):.1f} % MVC; "
+ f"unit with {CEILING_RANK} above it at {best.user_attrs.get('pct_mvc_ceiling_ranked', float('nan')):.1f} % "
+ f"(target {LAST_RECRUIT_PCT_MVC}), unrecruited {best.user_attrs.get('n_unrecruited', '?')}")
+ print(" MVC-ramp rates: " + ", ".join(
+ f"{lev} % MVC {best.user_attrs.get(f'mean_fr_{lev}_mvc_ramp', float('nan')):.1f}/{tgt:.1f} Hz"
+ for lev, tgt in HIGH_FORCE_RATE_TARGETS.items()))
+ for p in CHECKPOINTS:
+ sim_n = best.user_attrs.get(f"n_active_{p}", "?")
+ sim_fr = best.user_attrs.get(f"mean_fr_{p}", 0)
+ tgt = fug[p]
+ print(f" {p}% MVC: sim={sim_n}/{tgt['n_active']} active, "
+ f"FR={float(sim_fr):.1f}/{tgt['mean_fr_hz']:.1f} Hz")
+
+ # Save results (skip for parallel workers — coordinator saves)
+ result = {
+ "target": {
+ "model": target_model,
+ "threshold_mode": THRESHOLD_MODE,
+ "recruitment_range": RECRUITMENT_RANGE,
+ "max_recruitment": MAX_RECRUITMENT,
+ "last_recruit_pct_mvc": LAST_RECRUIT_PCT_MVC,
+ "ceiling_rank": CEILING_RANK,
+ "force_frequency_shape": list(FORCE_FREQUENCY_SHAPE),
+ "high_force_rate_targets": HIGH_FORCE_RATE_TARGETS if SCORE_HIGH_FORCE_RATES else None,
+ "de_luca_params": DE_LUCA_PARAMS,
+ "targets": {str(k): v for k, v in fug.items()},
+ "min_fr_hz": MIN_FR_HZ,
+ "max_fr_hz": MAX_FR_HZ,
+ },
+ "best_trial": {
+ "trial_number": best.number,
+ "objective_value": best.value,
+ "dd_neurons": best.params["dd_neurons"],
+ "conn_probability": best.params["conn_prob"],
+ "dd_drive__Hz": best.params["dd_drive"],
+ "synaptic_weight": SYNAPTIC_WEIGHT,
+ "gamma_shape": GAMMA_SHAPE,
+ "gfluctdv_scale": best.params.get("gfluctdv_scale", 0.0),
+ "isi_cv_mean": best.user_attrs.get("isi_cv_mean", 0.0),
+ "network": network_block(best.params.get("weight_exponent")),
+ "pct_mvc_last_recruited": best.user_attrs.get("pct_mvc_last_recruited"),
+ "pct_mvc_ceiling_ranked": best.user_attrs.get("pct_mvc_ceiling_ranked"),
+ **{f"mean_fr_{lev}_mvc_ramp": best.user_attrs.get(f"mean_fr_{lev}_mvc_ramp")
+ for lev in HIGH_FORCE_RATE_TARGETS},
+ "n_unrecruited": best.user_attrs.get("n_unrecruited"),
+ "threshold_hz": best.user_attrs.get("threshold_hz"),
+ **{f"n_active_{p}": best.user_attrs.get(f"n_active_{p}") for p in CHECKPOINTS},
+ **{f"mean_fr_{p}": best.user_attrs.get(f"mean_fr_{p}") for p in CHECKPOINTS},
+ },
+ "optimization": {
+ "n_motor_units": n_mu,
+ "n_trials": len(study.trials),
+ "ramp_mvc_pct": RAMP_MVC_PCT,
+ "checkpoints": CHECKPOINTS,
+ "checkpoint_weights": CHECKPOINT_WEIGHTS,
+ "ceiling_weight": CEILING_WEIGHT,
+ "ceiling_ramp_per_s": CEILING_RAMP_PER_S,
+ "mvc_rest_s": MVC_REST_S,
+ "high_force_rate_weight": HIGH_FORCE_RATE_WEIGHT,
+ "mvc_hold_s": MVC_HOLD_S,
+ "elapsed_s": elapsed,
+ },
+ }
+ if save_json:
+ out_path = results_dir / "dd_optimized_params.json"
+ with open(out_path, "w") as f:
+ json.dump(result, f, indent=2, default=float)
+ print(f"Saved: {out_path}")
+
+ return result
+
+
+def _save_best_from_study(study, n_mu, target_muscle, target_model,
+ fug_targets, elapsed):
+ """Extract best trial from an Optuna study and save to JSON."""
+ best = study.best_trial
+ results_dir = RESULTS_BASE / f"dd_{target_muscle}_{n_mu}mu_{RESULTS_TAG}"
+
+ print(f"\nBest trial {best.number} (err={best.value:.3f}, {elapsed:.0f}s):")
+ print(f" dd_neurons={best.params['dd_neurons']}, "
+ f"conn_prob={best.params['conn_prob']:.3f}, "
+ f"dd_drive={best.params['dd_drive']:.0f} Hz, "
+ f"gfluctdv={best.params.get('gfluctdv_scale', 0):.3f}, "
+ f"weight exponent={best.params.get('weight_exponent', 1.0):.2f}")
+ print(f" ISI CV: {best.user_attrs.get('isi_cv_mean', 0)*100:.1f}% "
+ f"(target: {ISI_CV_TARGET[0]*100:.0f}-{ISI_CV_TARGET[1]*100:.0f}%)")
+ print(f" Last unit recruited at {best.user_attrs.get('pct_mvc_last_recruited', float('nan')):.1f} % MVC; "
+ f"unit with {CEILING_RANK} above it at {best.user_attrs.get('pct_mvc_ceiling_ranked', float('nan')):.1f} % "
+ f"(target {LAST_RECRUIT_PCT_MVC}), unrecruited {best.user_attrs.get('n_unrecruited', '?')}")
+ print(" MVC-ramp rates: " + ", ".join(
+ f"{lev} % MVC {best.user_attrs.get(f'mean_fr_{lev}_mvc_ramp', float('nan')):.1f}/{tgt:.1f} Hz"
+ for lev, tgt in HIGH_FORCE_RATE_TARGETS.items()))
+ for p in CHECKPOINTS:
+ sim_n = best.user_attrs.get(f"n_active_{p}", "?")
+ sim_fr = best.user_attrs.get(f"mean_fr_{p}", 0)
+ tgt = fug_targets[p]
+ print(f" {p}% MVC: sim={sim_n}/{tgt['n_active']} active, "
+ f"FR={float(sim_fr):.1f}/{tgt['mean_fr_hz']:.1f} Hz")
+
+ result = {
+ "target": {
+ "model": target_model,
+ "threshold_mode": THRESHOLD_MODE,
+ "recruitment_range": RECRUITMENT_RANGE,
+ "max_recruitment": MAX_RECRUITMENT,
+ "last_recruit_pct_mvc": LAST_RECRUIT_PCT_MVC,
+ "ceiling_rank": CEILING_RANK,
+ "force_frequency_shape": list(FORCE_FREQUENCY_SHAPE),
+ "high_force_rate_targets": HIGH_FORCE_RATE_TARGETS if SCORE_HIGH_FORCE_RATES else None,
+ "de_luca_params": DE_LUCA_PARAMS,
+ "targets": {str(k): v for k, v in fug_targets.items()},
+ "min_fr_hz": MIN_FR_HZ,
+ "max_fr_hz": MAX_FR_HZ,
+ },
+ "best_trial": {
+ "trial_number": best.number,
+ "objective_value": best.value,
+ "dd_neurons": best.params["dd_neurons"],
+ "conn_probability": best.params["conn_prob"],
+ "dd_drive__Hz": best.params["dd_drive"],
+ "synaptic_weight": SYNAPTIC_WEIGHT,
+ "gamma_shape": GAMMA_SHAPE,
+ "gfluctdv_scale": best.params.get("gfluctdv_scale", 0.0),
+ "isi_cv_mean": best.user_attrs.get("isi_cv_mean", 0.0),
+ "network": network_block(best.params.get("weight_exponent")),
+ "pct_mvc_last_recruited": best.user_attrs.get("pct_mvc_last_recruited"),
+ "pct_mvc_ceiling_ranked": best.user_attrs.get("pct_mvc_ceiling_ranked"),
+ **{f"mean_fr_{lev}_mvc_ramp": best.user_attrs.get(f"mean_fr_{lev}_mvc_ramp")
+ for lev in HIGH_FORCE_RATE_TARGETS},
+ "n_unrecruited": best.user_attrs.get("n_unrecruited"),
+ "threshold_hz": best.user_attrs.get("threshold_hz"),
+ **{f"n_active_{p}": best.user_attrs.get(f"n_active_{p}") for p in CHECKPOINTS},
+ **{f"mean_fr_{p}": best.user_attrs.get(f"mean_fr_{p}") for p in CHECKPOINTS},
+ },
+ "optimization": {
+ "n_motor_units": n_mu,
+ "n_trials": len([t for t in study.trials
+ if t.state == optuna.trial.TrialState.COMPLETE]),
+ "ramp_mvc_pct": RAMP_MVC_PCT,
+ "checkpoints": CHECKPOINTS,
+ "checkpoint_weights": CHECKPOINT_WEIGHTS,
+ "ceiling_weight": CEILING_WEIGHT,
+ "ceiling_ramp_per_s": CEILING_RAMP_PER_S,
+ "mvc_rest_s": MVC_REST_S,
+ "high_force_rate_weight": HIGH_FORCE_RATE_WEIGHT,
+ "mvc_hold_s": MVC_HOLD_S,
+ "elapsed_s": elapsed,
+ },
+ }
+ out_path = results_dir / "dd_optimized_params.json"
+ with open(out_path, "w") as f:
+ json.dump(result, f, indent=2, default=float)
+ print(f"Saved: {out_path}")
+ return result
+
+
+def run_parallel(args):
+ """Launch parallel worker processes against a shared Optuna DB.
+
+ Each worker is a separate Python/NEURON process. Optuna coordinates
+ trial assignment through the shared SQLite database. Workers use
+ different TPE sampler seeds for diverse exploration.
+ """
+ from myogen.utils.parallel_budget import estimate_safe_njobs
+
+ n_mu = args.n_mu
+ n_trials = args.n_trials
+ target_muscle = args.target_muscle
+ target_model = args.target_model
+ profile = target_muscle
+
+ # Determine worker count
+ if args.n_workers:
+ n_workers = args.n_workers
+ else:
+ n_workers = estimate_safe_njobs(
+ n_mu=n_mu, dd_neurons=500, # worst-case from search range
+ max_jobs=8,
+ )
+
+ trials_per_worker = math.ceil(n_trials / n_workers)
+
+ # Set up results directory and DB
+ results_dir = RESULTS_BASE / f"dd_{target_muscle}_{n_mu}mu_{RESULTS_TAG}"
+ results_dir.mkdir(exist_ok=True, parents=True)
+ db_path = results_dir / "optuna_multilevel.db"
+ if args.fresh and db_path.exists():
+ db_path.unlink()
+ print(f" Removed existing {db_path} (--fresh)")
+
+ storage = f"sqlite:///{db_path}"
+ study_name = f"dd_multilevel_{n_mu}mu"
+
+ # Pre-create study so workers don't race on creation
+ optuna.logging.set_verbosity(optuna.logging.WARNING)
+ optuna.create_study(
+ direction="minimize",
+ sampler=optuna.samplers.TPESampler(seed=42, multivariate=True),
+ study_name=study_name,
+ storage=storage,
+ load_if_exists=True,
+ )
+
+ # Compute targets for final reporting
+ set_random_seed(42)
+ thresholds, target_thresholds = _thresholds(profile, n_mu)
+ if target_model == "fuglevand":
+ fug = compute_fr_targets(target_thresholds, mvc_levels=CHECKPOINTS,
+ model="fuglevand",
+ min_fr_hz=MIN_FR_HZ, max_fr_hz=MAX_FR_HZ)
+ elif target_model == "experimental":
+ fug = compute_fr_targets(target_thresholds, mvc_levels=CHECKPOINTS,
+ model="experimental",
+ muscle_name=target_muscle)
+ else:
+ fug = compute_fr_targets(target_thresholds, mvc_levels=CHECKPOINTS,
+ model="de_luca", **DE_LUCA_PARAMS)
+
+ print(f"\n{'='*60}")
+ print(f"Parallel DD optimization: {target_muscle} {n_mu} MUs")
+ print(f"Workers: {n_workers}, trials/worker: {trials_per_worker}")
+ print(f"Total trials: ~{n_workers * trials_per_worker} "
+ f"(requested {n_trials})")
+ print(f"Target model: {target_model}")
+ print(f"{'='*60}\n")
+
+ # Spawn workers
+ script = str(Path(__file__).resolve())
+ python = sys.executable
+ log_dir = results_dir / "worker_logs"
+ log_dir.mkdir(exist_ok=True)
+
+ procs = []
+ log_files = []
+ t0 = time.time()
+ for w in range(n_workers):
+ log_path = log_dir / f"worker_{w}.log"
+ cmd = [
+ python, script,
+ "--n-mu", str(n_mu),
+ "--n-trials", str(trials_per_worker),
+ "--target-model", target_model,
+ "--target-muscle", target_muscle,
+ "--threshold-mode", THRESHOLD_MODE,
+ "--recruitment-range", str(RECRUITMENT_RANGE),
+ "--max-recruitment", str(MAX_RECRUITMENT),
+ *(["--last-recruit-mvc", str(LAST_RECRUIT_PCT_MVC)] if LAST_RECRUIT_PCT_MVC is not None else []),
+ "--ceiling-rank", str(CEILING_RANK),
+ "--force-frequency", next(k for k, v in FORCE_FREQUENCY_SHAPES.items() if v == FORCE_FREQUENCY_SHAPE),
+ *(["--high-force-rates"] if SCORE_HIGH_FORCE_RATES else []),
+ *(["--fit-weight-exponent"] if FIT_WEIGHT_EXPONENT else []),
+ "--de-luca-params", json.dumps(DE_LUCA_PARAMS),
+ "--results-tag", RESULTS_TAG,
+ "--_worker-seed", str(42 + w),
+ ]
+ lf = open(log_path, "w")
+ p = subprocess.Popen(cmd, stdout=lf, stderr=subprocess.STDOUT)
+ procs.append((w, p))
+ log_files.append(lf)
+ print(f" Worker {w} started (PID {p.pid}, seed {42 + w})")
+
+ # Wait for all workers with periodic progress reporting
+ target_total = n_workers * trials_per_worker
+ while True:
+ # Check every 30s, but break immediately when all are done
+ for _ in range(30):
+ time.sleep(1)
+ if all(p.poll() is not None for _, p in procs):
+ break
+
+ try:
+ study = optuna.load_study(study_name=study_name, storage=storage)
+ completed = len([t for t in study.trials
+ if t.state == optuna.trial.TrialState.COMPLETE])
+ best_val = study.best_value if completed > 0 else float("inf")
+ elapsed = time.time() - t0
+ rate = completed / elapsed * 3600 if elapsed > 0 else 0
+ print(f" [{elapsed / 60:.0f}m] {completed}/{target_total} trials "
+ f"complete, best err={best_val:.3f}, "
+ f"rate={rate:.0f} trials/hr", flush=True)
+ except Exception:
+ pass
+
+ if all(p.poll() is not None for _, p in procs):
+ break
+
+ # Close log files
+ for lf in log_files:
+ lf.close()
+
+ elapsed = time.time() - t0
+ failures = [(w, p.returncode) for w, p in procs if p.returncode != 0]
+ if failures:
+ print(f"\n WARNING: {len(failures)} worker(s) failed:")
+ for w, rc in failures:
+ print(f" Worker {w}: exit code {rc} "
+ f"(see {log_dir / f'worker_{w}.log'})")
+
+ # Load final study and save best trial
+ study = optuna.load_study(study_name=study_name, storage=storage)
+ completed = len([t for t in study.trials
+ if t.state == optuna.trial.TrialState.COMPLETE])
+ print(f"\nOptimization complete: {completed} trials in {elapsed:.0f}s")
+ _save_best_from_study(study, n_mu, target_muscle, target_model,
+ fug, elapsed)
+
+
+def main():
+ load_nmodl_mechanisms()
+ h.secondorder = 2
+
+ ap = argparse.ArgumentParser()
+ ap.add_argument("--n-mu", type=int, default=200,
+ help="MU count to optimize for")
+ ap.add_argument("--n-trials", type=int, default=100,
+ help="Optuna trials")
+ ap.add_argument("--sweep", action="store_true",
+ help="Sweep all MU counts [150, 200, 250, 300]")
+ ap.add_argument("--fresh", action="store_true",
+ help="Delete existing Optuna database and start fresh")
+ ap.add_argument("--parallel", action="store_true",
+ help="Run optimization with multiple parallel worker "
+ "processes. Each worker gets its own NEURON instance.")
+ ap.add_argument("--n-workers", type=int, default=None,
+ help="Number of parallel workers (default: auto from "
+ "parallel_budget, capped at 8)")
+ ap.add_argument("--target-model", type=str, default="de_luca",
+ choices=["fuglevand", "de_luca", "experimental"],
+ help="Rate coding target model. de_luca (default) is "
+ "more realistic than fuglevand linear. experimental "
+ "uses a muscle-specific lookup table.")
+ ap.add_argument("--target-muscle", type=str, default="vl",
+ help="Muscle profile key. Used for results directory, "
+ "experimental FR lookup, and profile selection.")
+ ap.add_argument("--_worker-seed", type=int, default=None,
+ help=argparse.SUPPRESS) # internal: per-worker sampler seed
+ ap.add_argument("--threshold-mode", default="combined",
+ choices=["fuglevand", "deluca", "konstantin", "combined"])
+ ap.add_argument("--recruitment-range", type=float, default=30.0)
+ ap.add_argument("--max-recruitment", type=float, default=1.0,
+ help="Excitation fraction at which the last unit is recruited (rate-coding targets).")
+ ap.add_argument("--last-recruit-mvc", type=float, default=None,
+ help="Force (%% MVC) at which the last unit should be recruited, scored on a separate MVC ramp. "
+ "Default: no ceiling term.")
+ ap.add_argument("--ceiling-rank", type=int, default=0,
+ help="Score the unit with K units recruited above it instead of the last unit (default 0).")
+ ap.add_argument("--force-frequency", choices=sorted(FORCE_FREQUENCY_SHAPES), default="fuglevand",
+ help="ForceModel force-frequency saturation: fuglevand (1993, default) [fuglevand1993] or human (fit to human "
+ "single motor units, force_utils.HUMAN_UNIT_FORCE_FREQUENCY).")
+ ap.add_argument("--fit-weight-exponent", action="store_true",
+ help="Also fit the DD weight exponent p in [0, 1]: weights scale as (k_max / k_i) ** p "
+ "(default p = 1, equal current per input at rest; see build_dd_mn_network).")
+ ap.add_argument("--high-force-rates", action="store_true",
+ help="Also score MVC-ramp firing rates at 40 and 80 %% MVC against De Luca et al. 1982 [deluca1982].")
+ ap.add_argument("--de-luca-params", type=str, default="{}",
+ help="JSON overrides for de_luca_targets, e.g. fr_onset_hz.")
+ ap.add_argument("--results-tag", type=str, default="polygon")
+ args = ap.parse_args()
+ global THRESHOLD_MODE, RECRUITMENT_RANGE, MAX_RECRUITMENT, LAST_RECRUIT_PCT_MVC, DE_LUCA_PARAMS, RESULTS_TAG
+ global CEILING_RANK, FORCE_FREQUENCY_SHAPE, SCORE_HIGH_FORCE_RATES, FIT_WEIGHT_EXPONENT
+ THRESHOLD_MODE, RECRUITMENT_RANGE = args.threshold_mode, args.recruitment_range
+ MAX_RECRUITMENT, RESULTS_TAG = args.max_recruitment, args.results_tag
+ DE_LUCA_PARAMS = json.loads(args.de_luca_params)
+ LAST_RECRUIT_PCT_MVC = args.last_recruit_mvc
+ CEILING_RANK, SCORE_HIGH_FORCE_RATES = args.ceiling_rank, args.high_force_rates
+ FIT_WEIGHT_EXPONENT = args.fit_weight_exponent
+ FORCE_FREQUENCY_SHAPE = FORCE_FREQUENCY_SHAPES[args.force_frequency]
+ if LAST_RECRUIT_PCT_MVC is None:
+ print("Recruitment ceiling term off (pass --last-recruit-mvc PCT to score the last unit's recruitment force).")
+
+ if args.parallel:
+ run_parallel(args)
+ return
+
+ profile = args.target_muscle
+ is_worker = args._worker_seed is not None
+
+ if args.sweep:
+ for n_mu in [150, 200, 250, 300]:
+ optimize_for_mu_count(n_mu, args.n_trials, profile,
+ target_model=args.target_model,
+ target_muscle=args.target_muscle,
+ fresh=args.fresh)
+ else:
+ optimize_for_mu_count(
+ args.n_mu, args.n_trials, profile,
+ target_model=args.target_model,
+ target_muscle=args.target_muscle,
+ fresh=args.fresh if not is_worker else False,
+ sampler_seed=args._worker_seed if is_worker else 42,
+ save_json=not is_worker,
+ )
+
+
+if __name__ == "__main__":
+ main()
diff --git a/examples/02_finetune/07_fdi_spikes_force.py b/examples/02_finetune/07_fdi_spikes_force.py
new file mode 100644
index 00000000..d4081155
--- /dev/null
+++ b/examples/02_finetune/07_fdi_spikes_force.py
@@ -0,0 +1,118 @@
+"""
+FDI Spike Trains and Force from an Optimised Descending Drive
+=============================================================
+
+Usage: python 07_fdi_spikes_force.py [--n-mu N] [--max-recruitment R]
+
+1. MVC: slow ramp of the drive from rest to 100 % (10 s) and a 3 s hold; force from the
+ Fuglevand (1993) [fuglevand1993] twitch model; MVC = mean force over the last 2 s of the hold.
+2. The excitation (drive fraction above threshold, as in the optimiser) at which the ramp
+ force first reaches 25 % MVC.
+3. Benchmark trapezoids with the sEMG harness timing (0.5 s rest, 2 s ramp, 4 s hold, 2 s ramp,
+ 0.5 s rest): one held at the 25 % MVC excitation, one at 25 % excitation. Each is saved as
+ spike times per unit (s) with its force.
+Thresholds, drive mapping and the network come from 06_optimize_dd_multilevel.py itself.
+"""
+import argparse
+import importlib.util
+import json
+import sys
+from pathlib import Path
+
+import numpy as np
+import quantities as pq
+
+here = Path(__file__).resolve().parent
+spec = importlib.util.spec_from_file_location("opt", here / "06_optimize_dd_multilevel.py")
+opt = importlib.util.module_from_spec(spec)
+spec.loader.exec_module(opt)
+
+from myogen import set_random_seed # noqa: E402
+from myogen.utils.nmodl import load_nmodl_mechanisms # noqa: E402
+
+load_nmodl_mechanisms()
+
+ap = argparse.ArgumentParser()
+ap.add_argument("results_dir", type=Path)
+ap.add_argument("out_dir", type=Path)
+ap.add_argument("--n-mu", type=int, default=120)
+ap.add_argument("--threshold-mode", default=None, help="override the results file (older files lack it)")
+ap.add_argument("--recruitment-range", type=float, default=None)
+ap.add_argument("--max-recruitment", type=float, default=None)
+ap.add_argument("--extra-trap", action="append", default=[], metavar="NAME:LEVEL",
+ help="additional trapezoid at a fixed excitation level (e.g. from an external force model)")
+args = ap.parse_args()
+args.out_dir.mkdir(parents=True, exist_ok=True)
+
+res = json.load(open(args.results_dir / "dd_optimized_params.json"))
+tgt, best = res["target"], res["best_trial"]
+opt.THRESHOLD_MODE = args.threshold_mode or tgt.get("threshold_mode", "combined")
+opt.RECRUITMENT_RANGE = args.recruitment_range or float(tgt.get("recruitment_range", 30.0))
+opt.MAX_RECRUITMENT = args.max_recruitment or float(tgt.get("max_recruitment", 1.0))
+opt.FORCE_FREQUENCY_SHAPE = tuple(tgt.get("force_frequency_shape", opt.FORCE_FREQUENCY_SHAPE)) # as calibrated
+opt.CEILING_RANK = int(tgt.get("ceiling_rank", 0))
+set_random_seed(42)
+thresholds, _ = opt._thresholds("fdi", args.n_mu)
+dd = {k: best[k] for k in ("dd_neurons", "conn_probability", "dd_drive__Hz", "synaptic_weight", "gamma_shape",
+ "gfluctdv_scale")}
+dd["network"] = best.get("network") # absent in pre-2026-11-01 files: run_dd_ramp refuses them
+thr_hz = float(best["threshold_hz"])
+dt = opt.TIMESTEP_MS
+FS_FORCE = opt.FORCE_FS_HZ
+
+
+def simulate(segments):
+ """Drive, spikes and force as in the optimiser (``opt.simulate_force``)."""
+ return opt.simulate_force(thresholds, dd, thr_hz, segments)
+
+
+def save(name, e, block, force, mvc):
+ trains = [np.asarray(st.rescale(pq.s).magnitude) for st in block.segments[0].spiketrains]
+ np.savez(args.out_dir / f"{name}.npz", spikes=np.array(trains, dtype=object), force=force / mvc,
+ excitation=e[:: int(round(1000 / dt / FS_FORCE)) or 1], fs_force=FS_FORCE, thresholds=thresholds,
+ dd=json.dumps(dd), threshold_hz=thr_hz, results_dir=str(args.results_dir))
+
+
+# 1. MVC and the force-excitation curve (the optimiser's MVC ramp, ``opt.run_mvc_ramp``)
+r = opt.run_mvc_ramp(thresholds, dd, thr_hz)
+e, block, f, mvc = r.e, r.block, r.force, r.mvc
+t = np.arange(f.size) / FS_FORCE
+e_t = np.interp(t, np.arange(e.size) * dt / 1000, e)
+ramp = (t > opt.MVC_REST_S) & (t < opt.MVC_REST_S + 1.0 / opt.CEILING_RAMP_PER_S)
+smooth = opt.smoothed_force(f, FS_FORCE)
+curve = {int(p): float(e_t[ramp][np.argmax(smooth[ramp] >= p / 100 * mvc)]) for p in (5, 10, 20, 25, 30, 50)}
+save("mvc_ramp", e, block, f, mvc)
+print("MVC force (model units):", mvc)
+print("excitation at % MVC:", curve)
+print(f"force-frequency shape {opt.FORCE_FREQUENCY_SHAPE}; last unit recruited at {r.pct_last:.1f} % MVC, "
+ f"unit with {opt.CEILING_RANK} above at {r.pct_ranked:.1f} %, unrecruited {r.n_unrecruited}; rates {r.rates}",
+ flush=True)
+
+# 2. Benchmark trapezoids
+summary = {"mvc": mvc, "excitation_at_pct_mvc": curve, "pct_mvc_last_recruited": r.pct_last,
+ "n_unrecruited": r.n_unrecruited, "pct_mvc_ceiling_ranked": r.pct_ranked, "ceiling_rank": opt.CEILING_RANK,
+ "force_frequency_shape": list(opt.FORCE_FREQUENCY_SHAPE), "rates_at_pct_mvc": r.rates}
+
+
+def trapezoid(name, level):
+ e, block, f = simulate([(0.5, 0.0, 0.0), (2.0, 0.0, level), (4.0, level, level), (2.0, level, 0.0), (0.5, 0.0, 0.0)])
+ t = np.arange(f.size) / FS_FORCE
+ hold = (t > 3.0) & (t < 6.0)
+ rates = [np.sum((st.magnitude > 3.0) & (st.magnitude < 6.0)) / 3.0 for st in block.segments[0].spiketrains]
+ active = [r for r in rates if r >= 1.0]
+ summary[name] = {"excitation": level, "hold_force_pct_mvc": float(100 * np.mean(f[hold]) / mvc),
+ "n_active": len(active), "mean_rate_hz": float(np.mean(active)) if active else 0.0}
+ save(name, e, block, f, mvc)
+ print(name, summary[name], flush=True)
+ return summary[name]["hold_force_pct_mvc"]
+
+
+f_exc = trapezoid("trap_25exc", 0.25)
+f_mvc = trapezoid("trap_25mvc", curve[25])
+if abs(f_mvc - 25.0) > 1.0 and f_mvc != f_exc:
+ # One secant step on the hold force (the ramp estimate lags the isometric hold).
+ trapezoid("trap_25mvc", curve[25] + (25.0 - f_mvc) * (curve[25] - 0.25) / (f_mvc - f_exc))
+for spec in args.extra_trap:
+ name, level = spec.split(":")
+ trapezoid(name, float(level))
+json.dump(summary, open(args.out_dir / "summary.json", "w"), indent=1)
diff --git a/examples/02_finetune/README.md b/examples/02_finetune/README.md
index a232c963..bcdd5c92 100644
--- a/examples/02_finetune/README.md
+++ b/examples/02_finetune/README.md
@@ -9,3 +9,10 @@ results saved by its predecessor.
- **02** uses those optimised parameters to compute the resulting muscle force.
- **03** re-optimises the drive to hit a specific force target (e.g. 30 % MVC).
- **04** extracts inter-spike intervals and their coefficient of variation across contraction ramps.
+
+Examples 06 and 07 are command-line tools and stand apart from that chain:
+
+- **06** fits the descending-drive network (and the Gfluctdv noise scale) to rate-coding targets at several
+ force levels from one trapezoidal ramp per trial, sequentially or with parallel Optuna workers.
+- **07** maps excitation to % MVC with `ForceModel` on an MVC ramp and writes spike trains for force-matched
+ trapezoids from a parameter set fitted by 06.
diff --git a/myogen/config/alpha_mn_default.yaml b/myogen/config/alpha_mn_default.yaml
index 8dfb5df9..d941daed 100644
--- a/myogen/config/alpha_mn_default.yaml
+++ b/myogen/config/alpha_mn_default.yaml
@@ -7,7 +7,7 @@ mode: "active" # Options: "active", "passive"
axon_velocities: [50, 65] # [min, max] conduction velocities (m/s)
axon_length: 0.6 # Length of the axon (mm)
gamma: 0.2 # Neuromodulation level (a.u.)
-lambda_factor: 1.0 # Lambda factor for Powers2017 model persistent sodium scaling
+lambda_factor: 1.0 # Lambda factor for Powers2017 [powersheckman2017] model persistent sodium scaling
initial_voltage__mV: -67 # Initial membrane voltage (mV)
spike_threshold__mV: 50.0 # Spike detection threshold (mV)
@@ -86,7 +86,7 @@ nerlab:
r_alpha_v_offset: 55.0 # mV
r_alpha_k: 4.0 # mV
-# Powers2017 Model Parameters
+# Powers2017 Model Parameters (cell model of [powersheckman2017], ModelDB 239582; values not checked against it)
# All ranges are specified as [min, max, curve]
powers2017:
# Soma parameters
diff --git a/myogen/simulator/core/emg/electrodes.py b/myogen/simulator/core/emg/electrodes.py
index ef4d457b..d7323a83 100644
--- a/myogen/simulator/core/emg/electrodes.py
+++ b/myogen/simulator/core/emg/electrodes.py
@@ -186,7 +186,7 @@ def _create_electrode_grid(self) -> None:
- self._inter_electrode_distance__mm / self._bending_radius__mm
)
- ## Rotated detection system (Farina, 2004), eq (36)
+ ## Rotated detection system (Farina, 2004), eq (36) [farina2004 eq. 36]
# Must use original (z, θ) for both transformations — save before overwriting
displacement = self._center_point__mm_deg[0] * np.ones(_pos_z.shape)
pos_z_orig = _pos_z.copy()
@@ -320,24 +320,24 @@ def get_H_sf(
elif self.differentiation_mode == "bipolar_longitudinal":
# Differential along muscle fiber direction (z-axis)
- # Apply coordinate transformation for rotation (Farina 2004, eq 38)
+ # Apply coordinate transformation for rotation (Farina 2004, eq 38) [farina2004 eq. 38]
alpha_rad = self._rotation_angle__deg * np.pi / 180
kz_new = ktheta_mesh_kzktheta / self._bending_radius__mm * np.sin(
alpha_rad
) + kz_mesh_kzktheta * np.cos(alpha_rad)
- # Spatial filter for longitudinal differential (Farina 2004, eq 30)
+ # Spatial filter for longitudinal differential (Farina 2004, eq 30) [farina2004 eq. 30]
# Two electrodes at ±d_z/2: H_sf = exp(j·k'_z·d_z/2) - exp(-j·k'_z·d_z/2)
half_ied = self._inter_electrode_distance__mm / 2
H_sf = np.exp(1j * kz_new * half_ied) - np.exp(-1j * kz_new * half_ied)
elif self.differentiation_mode == "bipolar_transversal":
# Differential around muscle circumference (theta-axis)
- # Apply coordinate transformation for rotation (Farina 2004, eq 38)
+ # Apply coordinate transformation for rotation (Farina 2004, eq 38) [farina2004 eq. 38]
alpha_rad = self._rotation_angle__deg * np.pi / 180
ktheta_new = ktheta_mesh_kzktheta * np.cos(
alpha_rad
) - kz_mesh_kzktheta * self._bending_radius__mm * np.sin(alpha_rad)
- # Spatial filter for transversal differential (Farina 2004, eq 30-31)
+ # Spatial filter for transversal differential (Farina 2004, eq 30-31) [farina2004 eq. 30-31]
# d_theta = IED / R_ele in radians
half_d_theta = self._inter_electrode_distance__mm / (2 * self._bending_radius__mm)
H_sf = np.exp(1j * ktheta_new * half_d_theta) - np.exp(
diff --git a/myogen/simulator/core/emg/fiber_simulation.py b/myogen/simulator/core/emg/fiber_simulation.py
index be665c10..8264b1d6 100644
--- a/myogen/simulator/core/emg/fiber_simulation.py
+++ b/myogen/simulator/core/emg/fiber_simulation.py
@@ -2,13 +2,13 @@
Unified fiber simulation module.
Combines iEMG and sEMG fiber simulation under a single Rosenfalck source
-model (D1=96, z in physical mm) with volume conductor kernels computed
+model (D1=96, z in physical mm) [rosenfalck1969] with volume conductor kernels computed
as 1D spatial impulse responses.
References
----------
-[1] Rosenfalck, P., 1969. Intra- and extracellular potential fields of active nerve and muscle fibres. Acta Physiol. Scand. Suppl. 321, 1-168.
-[2] Farina, D. et al., 2004. A surface EMG generation model with multilayer cylindrical description of the volume conductor. IEEE TBME 51(3), 415-426.
+[1] Rosenfalck, P., 1969. Intra- and extracellular potential fields of active nerve and muscle fibres. Acta Physiol. Scand. Suppl. 321, 1-168. [rosenfalck1969]
+[2] Farina, D. et al., 2004. A surface EMG generation model with multilayer cylindrical description of the volume conductor. IEEE TBME 51(3), 415-426. [farina2004]
"""
from __future__ import annotations
@@ -37,7 +37,7 @@ def rosenfalck_dVm_dz(z: np.ndarray, D1: float = 96.0) -> np.ndarray:
z : np.ndarray
Spatial coordinates along fiber in mm.
D1 : float, default=96.0
- Amplitude parameter in mV/mm^3 (Rosenfalck 1969).
+ Amplitude parameter in mV/mm^3 (Rosenfalck 1969) [rosenfalck1969].
Returns
-------
@@ -61,7 +61,7 @@ def compute_intramuscular_kernel(
"""
Compute the intramuscular volume conductor kernel using the Green's function.
- For an anisotropic infinite medium (Nandedkar & Stalberg 1983):
+ For an anisotropic infinite medium (Nandedkar & Stalberg 1983) [nandedkar1983]:
b(z) = 1/(4*pi*sigma_r) / sqrt(sigma_z/sigma_r * r^2 + (z - z_elec)^2)
Parameters
@@ -112,7 +112,7 @@ def compute_surface_kernel(
"""
Compute the surface volume conductor spatial kernel b(z).
- Extracts the Bessel volume conductor computation (Farina 2004) from the
+ Extracts the Bessel volume conductor computation (Farina 2004) [farina2004] from the
existing sEMG pipeline and adds an IFFT step to produce the spatial kernel.
Steps:
@@ -284,7 +284,7 @@ def simulate_fiber_hybrid(
"""
Simulate a single fiber using the hybrid approach.
- Time-domain Rosenfalck source (D1=96, physical mm) combined with the
+ Time-domain Rosenfalck source [rosenfalck1969] (D1=96, physical mm) combined with the
frequency-domain Bessel volume conductor. Same math as the old path
but with the correct, unified source parameterization.
@@ -327,7 +327,7 @@ def simulate_fiber_hybrid(
fiber_length__mm : float or None
Physical fiber length for IAP kernel evaluation in mm.
D1 : float, optional
- Rosenfalck amplitude parameter, default 96.0.
+ Rosenfalck amplitude parameter [rosenfalck1969], default 96.0.
Returns
-------
@@ -386,7 +386,7 @@ def simulate_fiber_hybrid(
else:
z = np.linspace(-(N - 1) * (v / Fs) / 2, (N - 1) * (v / Fs) / 2, N)
- # Raw Rosenfalck — D1=96, z in physical mm, NO z/=2 scaling
+ # Raw Rosenfalck [rosenfalck1969] — D1=96, z in physical mm, NO z/=2 scaling
source = rosenfalck_dVm_dz(z, D1=D1)
# Convert dVm/dz (mV/mm) to transmembrane line current (A) by scaling
@@ -415,7 +415,7 @@ def simulate_fiber_hybrid(
ones = np.ones((len(PSI_conj[0, :]), 1))
PSI_mesh_conj = np.dot(ones, PSI_conj)
- # Build I(kz, kt) using Farina 2004 eq (3)
+ # Build I(kz, kt) using Farina 2004 eq (3) [farina2004 eq. 3]
I_kzkt = 1j * np.multiply(
kz_mesh_kzkt / v,
np.multiply(PSI_mesh_conj, np.exp(-1j * kz_mesh_kzkt * zi)),
@@ -550,7 +550,7 @@ def simulate_fiber_unified(
duration_ms : float
Signal duration in ms.
D1 : float, optional
- Rosenfalck amplitude parameter, default 96.0.
+ Rosenfalck amplitude parameter [rosenfalck1969], default 96.0.
Returns
-------
@@ -565,7 +565,7 @@ def simulate_fiber_unified(
n_samples = int(duration_ms * Fs)
t = np.arange(n_samples) / Fs # time in ms
- # Rosenfalck spatial extent (effectively zero beyond ~10 mm)
+ # Rosenfalck [rosenfalck1969] spatial extent (effectively zero beyond ~10 mm)
W = 10.0
# Compute the source: dVm/dz on the kernel grid
diff --git a/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py b/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py
index 938d0d46..1182bf64 100644
--- a/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py
+++ b/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py
@@ -67,16 +67,16 @@ class IntramuscularEMG:
By default, the endplate is located at the center of the muscle (50% of the muscle length).
nmj_jitter__s : Quantity__s, default=35e-6 * pq.s
Standard deviation of neuromuscular junction jitter in seconds.
- Default is set to 35e-6 s as determined by Akhmadeev et al. 2020 [1].
+ Default is set to 35e-6 s as determined by Akhmadeev et al. 2020 [akhmadeev2020].
branch_cvs__m_per_s : tuple[Quantity__m_per_s, Quantity__m_per_s], default=(5.0 * pq.m / pq.s, 2.0 * pq.m / pq.s)
Conduction velocities for the two-layer model of the neuromuscular junction in m/s.
- Default is set to (5.0, 2.0) m/s as determined by Akhmadeev et al. 2020 [1].
+ Default is set to (5.0, 2.0) m/s as determined by Akhmadeev et al. 2020 [akhmadeev2020].
Note: The two-layer model is a simplification of the actual arborization pattern, but it is a good approximation for the purposes of this simulation.
- Follows the implementation of Akhmadeev et al. 2020 [1].
+ Follows the implementation of Akhmadeev et al. 2020 [akhmadeev2020].
nmj_coefficients : {"published", "reference_code"}, default="published"
z spreads of the two-layer end-plate model: the values printed in
- Akhmadeev et al. 2020 [1], eq. 5-6, or those of their reference simulator
+ Akhmadeev et al. 2020 [1], eq. 5-6 [akhmadeev2020 eq. 5-6], or those of their reference simulator
code (1.5-2 times wider; the default up to 0.13.0).
MUs_to_simulate : list[int], optional
Indices of motor units to simulate. If None, all motor units are simulated.
@@ -97,7 +97,7 @@ class IntramuscularEMG:
References
----------
- [1] 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 Transactions on Biomedical Engineering 67, 2005–2014. https://doi.org/10.1109/TBME.2019.2953680
+ [1] 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 Transactions on Biomedical Engineering 67, 2005–2014. https://doi.org/10.1109/TBME.2019.2953680 [akhmadeev2020]
"""
def __init__(
@@ -277,7 +277,7 @@ def _simulate_neuromuscular_junctions(self, verbose: bool = True) -> None:
Simulate neuromuscular junction distributions for all motor units.
Two-layer branch model of Akhmadeev et al. (2020) with the z spreads of
- their eq. 5-6, growing with motor-unit size.
+ their eq. 5-6 [akhmadeev2020 eq. 5-6], growing with motor-unit size.
"""
if not self._motor_units:
raise ValueError("Must call _initialize_motor_units() first")
@@ -826,14 +826,14 @@ def add_noise(
Ignored when ``noise_type="gaussian"``.
baseline_drift_rms_uv : float, default=0.0
Target RMS (post-HPF) of a band-limited 1/f^α baseline
- drift representing electrode/interface noise (Huigen 2002,
- Gondran 1996). Same units as the EMG block. Set to 0 (the
+ drift representing electrode/interface noise (Huigen 2002 [huigen2002],
+ Gondran 1996 [gondran1996]). Same units as the EMG block. Set to 0 (the
default) to disable, preserving legacy behaviour.
The spectral form is paper-constrained; the amplitude
defaults are *not* validated for intramuscular EMG —
calibrate against real recordings via
`myogen.utils.calibrate_baseline_drift_profile`.
- Broadband movement artifacts (0–20 Hz, De Luca 2010) are
+ Broadband movement artifacts (0–20 Hz, De Luca 2010 [deluca2010]) are
a separate phenomenon out of scope here.
Ignored when ``noise_type="gaussian"``.
baseline_drift_alpha : float, default=1.75
diff --git a/myogen/simulator/core/emg/intramuscular/motor_unit_sim.py b/myogen/simulator/core/emg/intramuscular/motor_unit_sim.py
index ded55f12..93adab6c 100644
--- a/myogen/simulator/core/emg/intramuscular/motor_unit_sim.py
+++ b/myogen/simulator/core/emg/intramuscular/motor_unit_sim.py
@@ -5,7 +5,7 @@
junction modeling, single fiber action potential (SFAP) calculation, and
motor unit action potential (MUAP) generation with realistic jitter.
-Based on the MU_Sim class from the MATLAB iemg_simulator.
+Based on the MU_Sim class from the MATLAB iemg_simulator [akhmadeev2020].
"""
import os
diff --git a/myogen/simulator/core/emg/nerve_branches.py b/myogen/simulator/core/emg/nerve_branches.py
index 2a0a5763..594a434c 100644
--- a/myogen/simulator/core/emg/nerve_branches.py
+++ b/myogen/simulator/core/emg/nerve_branches.py
@@ -1,4 +1,4 @@
-"""Two-layer motor-nerve branching model (Akhmadeev et al. 2020).
+"""Two-layer motor-nerve branching model (Akhmadeev et al. 2020) [akhmadeev2020].
Single authority for where a motor unit's end-plates sit along its fibres and
how far the action potential travels along the axon branches to reach them.
@@ -18,8 +18,8 @@
# (branch intercept, branch slope, arborisation intercept, arborisation slope), mm.
-# "reference_code": Akhmadeev's reference simulator (iemg_simulator, s08_cl_init_muaps.m).
-# "published": Akhmadeev et al. 2020, eq. 5-6 (a_mu = 1, b_mu = 2.5, a_sigma = 0.25, b_sigma = 1).
+# "reference_code": Akhmadeev's reference simulator (iemg_simulator, s08_cl_init_muaps.m) [akhmadeev2020].
+# "published": Akhmadeev et al. 2020, eq. 5-6 (a_mu = 1, b_mu = 2.5, a_sigma = 0.25, b_sigma = 1) [akhmadeev2020 eq. 5-6].
NMJ_SCHEDULE_COEFFICIENTS = {
"reference_code": (1.5, 4.0, 0.5, 1.5),
"published": (1.0, 2.5, 0.25, 1.0),
diff --git a/myogen/simulator/core/emg/surface/volume_conductor.py b/myogen/simulator/core/emg/surface/volume_conductor.py
index eee64ee7..8d8b2dc3 100644
--- a/myogen/simulator/core/emg/surface/volume_conductor.py
+++ b/myogen/simulator/core/emg/surface/volume_conductor.py
@@ -1,4 +1,4 @@
-"""Layered-cylinder volume conductor transfer function (Farina & Merletti 2004).
+"""Layered-cylinder volume conductor transfer function (Farina & Merletti 2004) [farina2004].
Single authority for the spatial transfer function ``H_vc(k_z, k_theta)`` of
the four-layer (bone / anisotropic muscle / fat / skin) cylindrical volume
@@ -94,7 +94,7 @@ def layered_cylinder_transfer(
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Surface potential transfer function of a fibre in a layered cylinder.
- Solves the Farina 2004 boundary-condition system for a line source at
+ Solves the Farina 2004 [farina2004] boundary-condition system for a line source at
radius ``R`` inside the anisotropic muscle layer and returns the skin
potential per unit source, ``H_vc(k_z, k_theta)``, on the full 2-D
frequency grid.
@@ -137,7 +137,7 @@ def layered_cylinder_transfer(
k_z = np.asarray(k_z, dtype=np.float64)
k_theta = np.asarray(k_theta, dtype=np.float64)
- # Model radii -> (Farina, 2004), Figure 1 - b)
+ # Model radii -> (Farina, 2004), Figure 1 - b) [farina2004 Fig. 1b]
th_muscle = r - th_fat - th_skin - r_bone
a = r_bone
b = r_bone + th_muscle
diff --git a/myogen/simulator/core/force/force_model.py b/myogen/simulator/core/force/force_model.py
index 4853dadd..871b3b01 100644
--- a/myogen/simulator/core/force/force_model.py
+++ b/myogen/simulator/core/force/force_model.py
@@ -17,15 +17,15 @@
SPIKE_TRAIN__Block,
)
-from .force_utils import fuglevand_gain
+from .force_utils import FUGLEVAND_FORCE_FREQUENCY, fuglevand_gain
@beartowertype
class ForceModel:
"""
- Force model based on Fuglevand et al. (1993) [1].
+ Force model based on Fuglevand et al. (1993) [fuglevand1993].
- This class implements the Fuglevand force generation model for motor unit pools,
+ This class implements the Fuglevand [fuglevand1993] force generation model for motor unit pools,
computing individual motor unit twitch responses and their nonlinear gain modulation
based on discharge rates. The model generates realistic force outputs from spike trains
using physiologically-based parameters.
@@ -48,17 +48,21 @@ class ForceModel:
twitch_force_range : float or None, default=None
Twitch-force range (RP in _[1]) used in equations 13 and 14. Must be greater
than 1. If None, defaults to the last/first recruitment-threshold ratio for
- backward compatibility. Fuglevand et al. (1993) use RP = 100 independently
+ backward compatibility. Fuglevand et al. (1993) [fuglevand1993] use RP = 100 independently
of the recruitment range RR = 30.
+ force_frequency_shape : tuple of float, default=(2.0, 3.0)
+ (a, c) of the force-frequency saturation S(r) = 1 - exp(-a r^c), r = contraction time / ISI.
+ The default is Fuglevand et al. (1993) eq. 17 [fuglevand1993 eq. 17]. ``force_utils.HUMAN_UNIT_FORCE_FREQUENCY``
+ (0.93, 1.29) is fitted to human single motor units, which keep gaining force above ~20 Hz.
Attributes
----------
peak_twitch_forces__unitless : np.ndarray
Peak twitch forces for each motor unit (unitless). Available after initialization.
- Computed according to Fuglevand equation 13.
+ Computed according to Fuglevand equation 13 [fuglevand1993 eq. 13].
contraction_times__samples : np.ndarray
Contraction times for each motor unit in samples of recording_frequency__Hz.
- Available after initialization. Computed according to Fuglevand equation 14.
+ Available after initialization. Computed according to Fuglevand equation 14 [fuglevand1993 eq. 14].
twitch_mat : np.ndarray
Twitch matrix sampled up to 10 times the longest contraction time, for plotting only.
Available after initialization.
@@ -87,7 +91,7 @@ class ForceModel:
----------
[1] Fuglevand, A. J., Winter, D. A., & Patla, A. E. (1993). Models of recruitment
and rate coding organization in motor-unit pools. J Neurophysiol, 70(6), 2470–2488.
- https://doi.org/10.1152/jn.1993.70.6.2470.
+ https://doi.org/10.1152/jn.1993.70.6.2470. [fuglevand1993]
Examples
--------
@@ -107,6 +111,7 @@ def __init__(
longest_duration_rise_time__ms: Quantity__ms = 90.0 * pq.ms,
contraction_time_range_factor: float = 3.0,
twitch_force_range: float | None = None,
+ force_frequency_shape: tuple[float, float] = FUGLEVAND_FORCE_FREQUENCY,
) -> None:
# Input validation
if len(recruitment_thresholds) == 0:
@@ -145,7 +150,13 @@ def __init__(
if twitch_force_range is not None and twitch_force_range <= 1.0:
raise ValueError(
f"twitch_force_range must be greater than 1.0, got {twitch_force_range}. "
- "This parameter determines the spread of peak twitch forces. Fuglevand et al. (1993) use 100."
+ "This parameter determines the spread of peak twitch forces. Fuglevand et al. (1993) [fuglevand1993] use 100."
+ )
+
+ force_frequency_shape = tuple(float(v) for v in force_frequency_shape)
+ if len(force_frequency_shape) != 2 or min(force_frequency_shape) <= 0:
+ raise ValueError(
+ f"force_frequency_shape must be two positive numbers (a, c), got {force_frequency_shape}."
)
# Immutable public access
@@ -154,6 +165,7 @@ def __init__(
self.longest_duration_rise_time__ms = longest_duration_rise_time__ms
self.contraction_time_range_factor = contraction_time_range_factor
self.twitch_force_range = twitch_force_range
+ self.force_frequency_shape = force_frequency_shape
# Private copies for internal modifications
self._recruitment_thresholds = recruitment_thresholds.copy()
@@ -186,10 +198,10 @@ def __init__(
def _compute_twitch_parameters(self) -> None:
"""
- Compute peak twitch forces and contraction times based on Fuglevand model.
+ Compute peak twitch forces and contraction times based on Fuglevand model [fuglevand1993].
This method calculates the physiological parameters for each motor unit
- according to the Fuglevand et al. (1993) model equations.
+ according to the Fuglevand et al. (1993) [fuglevand1993] model equations.
Results are stored in the `peak_twitch_forces` and `contraction_times`
properties after execution.
@@ -251,7 +263,7 @@ def _initialize_twitches(self) -> None:
def generate_force(self, spike_train__Block: SPIKE_TRAIN__Block, verbose: bool = True) -> FORCE__AnalogSignal:
"""
- Generate force output from motor unit spike trains using the Fuglevand model.
+ Generate force output from motor unit spike trains using the Fuglevand model [fuglevand1993].
This method simulates muscle force by converting spike trains into force output
through individual motor unit twitches with nonlinear gain modulation based on
@@ -363,7 +375,7 @@ def _generate_force(
P = self._peak_twitch_forces__unitless[n]
r = T / (np.diff(idx) * spiketrain_timestep__s)
x = np.zeros(L)
- x[idx] = np.concatenate(([1.0], fuglevand_gain(r)))
+ x[idx] = np.concatenate(([1.0], fuglevand_gain(r, *self.force_frequency_shape)))
# Exact sampled twitch impulse response, with no tail truncation.
a = np.exp(-spiketrain_timestep__s / T)
@@ -388,7 +400,7 @@ def peak_twitch_forces__unitless(self) -> np.ndarray:
Returns
-------
np.ndarray
- Array of peak twitch forces according to Fuglevand model equation 13.
+ Array of peak twitch forces according to Fuglevand model equation 13 [fuglevand1993 eq. 13].
Raises
------
@@ -410,7 +422,7 @@ def contraction_times__samples(self) -> np.ndarray:
Returns
-------
np.ndarray
- Array of contraction times according to Fuglevand model equation 14.
+ Array of contraction times according to Fuglevand model equation 14 [fuglevand1993 eq. 14].
Raises
------
diff --git a/myogen/simulator/core/force/force_model_vectorized.py b/myogen/simulator/core/force/force_model_vectorized.py
index d911b365..0b107b4c 100644
--- a/myogen/simulator/core/force/force_model_vectorized.py
+++ b/myogen/simulator/core/force/force_model_vectorized.py
@@ -28,5 +28,5 @@ class ForceModelVectorized(ForceModel):
----------
[1] Fuglevand, A. J., Winter, D. A., & Patla, A. E. (1993). Models of recruitment
and rate coding organization in motor-unit pools. J Neurophysiol, 70(6), 2470–2488.
- https://doi.org/10.1152/jn.1993.70.6.2470.
+ https://doi.org/10.1152/jn.1993.70.6.2470. [fuglevand1993]
"""
diff --git a/myogen/simulator/core/force/force_utils.py b/myogen/simulator/core/force/force_utils.py
index 2a88dc32..6b76a540 100644
--- a/myogen/simulator/core/force/force_utils.py
+++ b/myogen/simulator/core/force/force_utils.py
@@ -1,26 +1,35 @@
import numpy as np
+# Force-frequency saturation S(r) = 1 - exp(-a r^c). Keys in [brackets] are entries of docs/literature.md.
+FUGLEVAND_FORCE_FREQUENCY = (2.0, 3.0) # [fuglevand1993 p. 2475 eq. 16]
+# Fitted to human single motor units stimulated through their axons (mean plateau force vs contraction time x rate):
+# [macefield1996 Fig. 6C; CT Table 1], [fuglevand1999 Table 1, Fig. 6A], [mcnulty2000 Table 2; CT Table 1].
+# rms error 0.05 (Fuglevand 1993 [fuglevand1993]: 0.13); predicts twitch:tetanus 0.23 (measured 0.22-0.36; Fuglevand 1993: 0.11 [fuglevand1993 p. 2475]).
+HUMAN_UNIT_FORCE_FREQUENCY = (0.93, 1.29)
-def fuglevand_gain(r: np.ndarray) -> np.ndarray:
+
+def fuglevand_gain(r: np.ndarray, a: float = 2.0, c: float = 3.0) -> np.ndarray:
"""
- Compute the rate-dependent twitch gain from Fuglevand et al. (1993), eq. 17.
+ Compute the rate-dependent twitch gain from Fuglevand et al. (1993), eq. 17 [fuglevand1993 p. 2475 eqs. 16-17].
Parameters
----------
r : np.ndarray
Dimensionless contraction time divided by the preceding inter-spike interval.
Both times must be expressed in the same units.
+ a, c : float, default 2.0, 3.0
+ Shape of the saturation S(r) = 1 - exp(-a * r**c); Fuglevand et al. (1993) [fuglevand1993 p. 2475 eq. 16] use a = 2, c = 3.
+ ``HUMAN_UNIT_FORCE_FREQUENCY`` is a fit to human single motor units.
Returns
-------
np.ndarray
- Gain of 1 for r <= 0.4, otherwise (S(r) / r) / (S(0.4) / 0.4),
- where S(r) = 1 - exp(-2 * r**3).
+ Gain of 1 for r <= 0.4, otherwise (S(r) / r) / (S(0.4) / 0.4).
"""
r = np.asarray(r, dtype=float)
gain = np.ones_like(r)
mask = r > 0.4
- gain[mask] = ((1 - np.exp(-2 * r[mask] ** 3)) / r[mask]) / (
- (1 - np.exp(-2 * 0.4**3)) / 0.4
+ gain[mask] = ((1 - np.exp(-a * r[mask] ** c)) / r[mask]) / (
+ (1 - np.exp(-a * 0.4**c)) / 0.4
)
return gain
diff --git a/myogen/simulator/core/physiological_distribution.py b/myogen/simulator/core/physiological_distribution.py
index 0afb6a19..eff2c321 100644
--- a/myogen/simulator/core/physiological_distribution.py
+++ b/myogen/simulator/core/physiological_distribution.py
@@ -32,9 +32,9 @@ def __init__(
recruitment range (RR) and, for some models, additional parameters.
Following models are available:
- - Fuglevand et al. (1993) [1]
- - De Luca & Contessa (2012) [2]
- - Akhmadeev et al. (2020) [3]
+ - Fuglevand et al. (1993) [fuglevand1993]
+ - De Luca & Contessa (2012) [delucacontessa2012]
+ - Akhmadeev et al. (2020) [akhmadeev2020]
- Combined model
Parameters
@@ -70,25 +70,25 @@ def __init__(
References
----------
- [1] Fuglevand, A.J., Winter, D.A., Patla, A.E., 1993. Models of recruitment and rate coding organization in motor-unit pools. Journal of Neurophysiology 70, 2470-2488. https://doi.org/10.1152/jn.1993.70.6.2470
- [2] De Luca, C.J., Contessa, P., 2012. Hierarchical control of motor units in voluntary contractions. Journal of Neurophysiology 107, 178-195. https://doi.org/10.1152/jn.00961.2010
- [3] 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 Transactions on Biomedical Engineering 67, 2005-2014. https://doi.org/10.1109/TBME.2019.2953680
+ [1] Fuglevand, A.J., Winter, D.A., Patla, A.E., 1993. Models of recruitment and rate coding organization in motor-unit pools. Journal of Neurophysiology 70, 2470-2488. https://doi.org/10.1152/jn.1993.70.6.2470 [fuglevand1993]
+ [2] De Luca, C.J., Contessa, P., 2012. Hierarchical control of motor units in voluntary contractions. Journal of Neurophysiology 107, 178-195. https://doi.org/10.1152/jn.00961.2010 [delucacontessa2012]
+ [3] 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 Transactions on Biomedical Engineering 67, 2005-2014. https://doi.org/10.1109/TBME.2019.2953680 [akhmadeev2020]
Notes
-----
- **fuglevand** : Fuglevand et al. (1993) [1] exponential model
+ **fuglevand** : Fuglevand et al. (1993) [fuglevand1993] exponential model
$$rt(i) = \exp\left( \frac{i \cdot \ln(RR)}{N} \right) / 100$$
where $i = 1, 2, \ldots, N$
- **deluca** : De Luca & Contessa (2012) [2] model with slope correction
+ **deluca** : De Luca & Contessa (2012) [delucacontessa2012] model with slope correction
$$rt(i) = \frac{b \cdot i}{N} \cdot \exp\left(\frac{i \cdot \ln(RR / b)}{N}\right) / 100$$
where $b$ = ``deluca__slope``, $i = 1, 2, \ldots, N$
- **konstantin** : Akhmadeev et al. (2020) [3] model (named after the first author's given name) allowing explicit maximum threshold control
+ **konstantin** : Akhmadeev et al. (2020) [akhmadeev2020] model (named after the first author's given name) allowing explicit maximum threshold control
$$
\begin{aligned}
diff --git a/myogen/simulator/neuron/cells.py b/myogen/simulator/neuron/cells.py
index 8847d377..0b6caf92 100644
--- a/myogen/simulator/neuron/cells.py
+++ b/myogen/simulator/neuron/cells.py
@@ -197,14 +197,14 @@ def create_synapses(
Notes
-----
Default parameters (tau1=0.2, tau2=0.3, weight=0.6) are based on
- experimental data from Sanches (2018). For NERLab model compatibility,
+ experimental data from Sanches (2018) [sanches2018]. For NERLab model compatibility,
the reversal potential is automatically set to 70 mV when the cell
has the "model" attribute set to "NERLab".
References
----------
- - Sanches et al. (2018): Synaptic parameter measurements
- - Segev et al. (1990): Alpha function synaptic modeling
+ - Sanches et al. (2018): Synaptic parameter measurements [sanches2018]
+ - Segev et al. (1990): Alpha function synaptic modeling [segev1990]
"""
assert decay_time_constant__ms > rise_time_constant__ms, (
"decay_time_constant__ms must be greater than tau1 for proper kinetics"
@@ -232,7 +232,7 @@ class INgII(_Cell):
are part of the spinal reflex circuitry and contribute to muscle tone regulation
and reciprocal inhibition between antagonist muscles.
- The cell model is based on experimental data from Bui et al. (2003) for
+ The cell model is based on experimental data from Bui et al. (2003) [bui2003] for
Ia inhibitory interneurons, with single-compartment soma containing
physiologically realistic ion channels including sodium, potassium, h-current,
and calcium-dependent afterhyperpolarization channels.
@@ -253,7 +253,7 @@ class INgII(_Cell):
Notes
-----
The geometry and biophysics are based on cat spinal cord measurements:
- - Membrane area: ~81,390 um² (Bui et al., 2003)
+ - Membrane area: ~81,390 um² (Bui et al., 2003) [bui2003]
- Input resistance: ~20 MΩ (from g_pas = 5e-5 S/cm²)
- Resting potential: -71 mV
- Active channels: Na+, K+, h-current, Ca²⁺-dependent K+
@@ -261,7 +261,7 @@ class INgII(_Cell):
References
----------
Bui et al. (2003). Morphological and electrophysiological characterization
- of Ia inhibitory interneurons in cat spinal cord.
+ of Ia inhibitory interneurons in cat spinal cord. [bui2003]
"""
_ids2 = itertools.count(0)
@@ -278,11 +278,11 @@ def _define_geometry(self):
"""
Set soma geometry based on experimental morphometric data.
- Implements measurements from Bui et al. (2003) for Ia inhibitory
+ Implements measurements from Bui et al. (2003) [bui2003] for Ia inhibitory
interneurons, calculating diameter from measured membrane surface area
assuming spherical geometry.
"""
- # Data from Bui et al., 2003 Ia Inhibitory interneurons
+ # Data from Bui et al., 2003 Ia Inhibitory interneurons [bui2003]
Amu = 81390 + 3113 # Mean membrane area (um²)
Aci = 1.96 * (891.5 + 46.141) / np.sqrt(8) # Area confidence interval
i = 0 # Use lower bound of area distribution
@@ -755,8 +755,8 @@ class AlphaMN(_Cell):
References
----------
- - Powers & Binder (2001): Motor neuron biophysical characterization
- - Powers et al. (2017): Updated ion channel models and parameters
+ - Powers & Binder (2001): Motor neuron biophysical characterization [powersbinder2001]
+ - Powers & Heckman (2017): ion channel mechanisms of their model, ModelDB 239582 [powersheckman2017]
"""
_ids2 = itertools.count(0)
@@ -875,7 +875,7 @@ def _define_biophysics(self, gamma=1):
Insert ion channels and set biophysical parameters for motor neuron.
Implements detailed voltage-gated ion channel models based on experimental
- characterization. Supports both Powers2017 and NERLab parameter sets with
+ characterization. Supports both Powers2017 [powersheckman2017] and NERLab parameter sets with
different channel implementations and calibrations.
The Powers2017 model includes the complete complement of motor neuron
diff --git a/myogen/simulator/neuron/network.py b/myogen/simulator/neuron/network.py
index 666d7af3..6e76ca24 100644
--- a/myogen/simulator/neuron/network.py
+++ b/myogen/simulator/neuron/network.py
@@ -7,6 +7,7 @@
from typing import Optional
+import numpy as np
import quantities as pq
from beartype.typing import Callable
from neuron import h
@@ -1107,6 +1108,67 @@ def connect(
return netcons
+ def connect_fan_in(
+ self,
+ source: str,
+ target: str,
+ n_inputs: int,
+ weights__uS,
+ delay__ms: Quantity__ms = DEFAULT_SYNAPTIC_DELAY,
+ ) -> list:
+ """
+ Give every target neuron exactly ``n_inputs`` excitatory inputs from ``source``.
+
+ For each target, in population order, ``n_inputs`` source neurons are drawn without
+ replacement from MyoGen's random generator.
+
+ Parameters
+ ----------
+ source, target : str
+ Population names.
+ n_inputs : int
+ Inputs per target neuron, 1 to the size of ``source``.
+ weights__uS : Quantity
+ One weight for all targets, or one per target neuron.
+ delay__ms : Quantity, optional
+ Synaptic delay, by default 1 ms.
+
+ Returns
+ -------
+ list[h.NetCon]
+ The connections, grouped by target neuron.
+ """
+ if source not in self.populations or target not in self.populations:
+ raise ValueError(f"Populations '{source}' and '{target}' must both exist")
+ sources = list(self.populations[source])
+ targets = list(self.populations[target])
+ if not 1 <= n_inputs <= len(sources):
+ raise ValueError(f"n_inputs must be 1 to {len(sources)}, got {n_inputs}")
+ weights = np.broadcast_to(
+ np.asarray(pq.Quantity(weights__uS).rescale(pq.uS).magnitude, dtype=float), (len(targets),))
+ delay_value = _to_float(delay__ms, pq.ms)
+
+ netcons = []
+ for target_neuron, weight in zip(targets, weights):
+ for j in get_random_generator().choice(len(sources), size=n_inputs, replace=False):
+ netcon = _create_netcon(sources[j], _select_synapse(target_neuron, inhibitory=False),
+ neuron_id=sources[j].global__ID, synaptic_delay=delay_value)
+ netcon.weight[0] = weight
+ netcons.append(netcon)
+
+ self.connections.append(
+ {
+ "type": "fan_in",
+ "source": source,
+ "target": target,
+ "n_inputs": n_inputs,
+ "weight__uS": weights.tolist(),
+ "delay__ms": delay_value,
+ }
+ )
+ self._netcons_by_connection[(source, target)] = netcons
+ return netcons
+
def connect_to_muscle(
self,
source: str,
@@ -1417,6 +1479,12 @@ def print_network(self):
f" {i + 1}. {conn['source']} → {conn['target']} [1-to-1] "
f"(p={conn['probability']}, w={conn['weight__uS']}uS)"
)
+ elif conn["type"] == "fan_in":
+ w = conn["weight__uS"]
+ print(
+ f" {i + 1}. {conn['source']} → {conn['target']} [fan-in {conn['n_inputs']}] "
+ f"(w={min(w):.4g}-{max(w):.4g}uS)"
+ )
if __name__ == "__main__":
diff --git a/myogen/simulator/neuron/populations/motor_neurons.py b/myogen/simulator/neuron/populations/motor_neurons.py
index 8b28aa1c..2897e399 100644
--- a/myogen/simulator/neuron/populations/motor_neurons.py
+++ b/myogen/simulator/neuron/populations/motor_neurons.py
@@ -24,7 +24,7 @@ class AlphaMN__Pool(_Pool):
Container for a population of alpha motor neurons.
Manages a collection of AlphaMN (alpha motor neuron) cells with different
- biophysical models: ModALS or Powers2017. These cells form the final
+ biophysical models: ModALS or Powers2017 [powersheckman2017]. These cells form the final
common pathway for motor control.
Parameters
@@ -54,7 +54,7 @@ class AlphaMN__Pool(_Pool):
cell_index : Optional[int], optional
Specific cell index to create (creates only one cell), by default None.
lambda_factor : float, optional
- Lambda factor for Powers2017 model persistent sodium scaling, by default 1.0.
+ Lambda factor for Powers2017 [powersheckman2017] model persistent sodium scaling, by default 1.0.
initial_voltage__mV : float or list[float], optional
Initial membrane voltage (mV), by default -67.
spike_threshold__mV : float, optional
diff --git a/myogen/utils/emg_noise.py b/myogen/utils/emg_noise.py
index 0cbcff1b..c8802daf 100644
--- a/myogen/utils/emg_noise.py
+++ b/myogen/utils/emg_noise.py
@@ -337,12 +337,12 @@ def _baseline_drift(
Spectral form is paper-constrained: PSD ∝ ``1/f^alpha`` over an
explicit band ``[low_hz, high_hz]`` and zero outside, consistent
- with Huigen, Peper & Grimbergen 2002 (DOI 10.1007/BF02344216) and
- Gondran et al. 1996 (DOI 10.1007/BF02523851), which report α in
+ with Huigen, Peper & Grimbergen 2002 (DOI 10.1007/BF02344216) [huigen2002] and
+ Gondran et al. 1996 (DOI 10.1007/BF02523851) [gondran1996], which report α in
the 1.5–2.0 range for electrode/interface noise. The hard band
limit (rather than a Lorentzian rolloff) keeps the model scoped
to baseline wander; broadband movement artifacts (0–20 Hz per
- De Luca et al. 2010, DOI 10.1016/j.jbiomech.2010.01.027) are a
+ De Luca et al. 2010, DOI 10.1016/j.jbiomech.2010.01.027 [deluca2010]) are a
separate phenomenon and require a different contaminant model.
The drift represents the residual that survives the amplifier's
@@ -445,15 +445,15 @@ def generate_realistic_noise(
Spectral form is paper-constrained: PSD ∝ 1/f^α over an
explicit band ``[baseline_drift_low_hz, baseline_drift_high_hz]``,
consistent with electrode/interface noise studies (Huigen,
- Peper & Grimbergen 2002, DOI 10.1007/BF02344216; Gondran
- et al. 1996, DOI 10.1007/BF02523851). The default α=1.75 and
+ Peper & Grimbergen 2002, DOI 10.1007/BF02344216 [huigen2002]; Gondran
+ et al. 1996, DOI 10.1007/BF02523851 [gondran1996]). The default α=1.75 and
``high_hz=1.0`` are the midpoint and upper bound of the
reported electrode-noise regime, not validated physiological
constants for intramuscular EMG. Calibrate against real
recordings via `calibrate_baseline_drift_profile` if
you need amplitude-accurate simulation. Movement artifacts
(broadband, 0–20 Hz per De Luca et al. 2010, DOI
- 10.1016/j.jbiomech.2010.01.027) are a separate phenomenon
+ 10.1016/j.jbiomech.2010.01.027 [deluca2010]) are a separate phenomenon
and would need their own contaminant model.
Note: **SNR contract**: drift is **additive on top** of the
@@ -554,7 +554,7 @@ def generate_realistic_noise(
noise *= noise_rms / current_rms
# 5b. Baseline drift — slow low-frequency wander attributed to
- # the electrode/interface (Huigen 2002; Gondran 1996). Additive
+ # the electrode/interface (Huigen 2002; Gondran 1996) [huigen2002] [gondran1996]. Additive
# in absolute units, post-HPF (i.e., what actually survives into
# the recorded signal). Off by default.
if baseline_drift_rms_uv > 0:
@@ -676,7 +676,7 @@ def add_realistic_noise(
is not modelled.
baseline_drift_alpha : float, default 1.75
Drift PSD slope α. Midpoint of the [1.5, 2.0] electrode-noise
- regime (Huigen 2002, Gondran 1996).
+ regime (Huigen 2002, Gondran 1996) [huigen2002] [gondran1996].
baseline_drift_low_hz : float or None, default None
Lower edge of the drift band (lowest nonzero FFT bin when None).
baseline_drift_high_hz : float, default 1.0
@@ -1123,7 +1123,7 @@ def calibrate_baseline_drift_profile(
Note: The model's spectral form (PSD ∝ 1/f^α over a bounded band)
is paper-constrained for electrode/interface noise (Huigen
- 2002, Gondran 1996). The fitted parameters become validated
+ 2002, Gondran 1996) [huigen2002] [gondran1996]. The fitted parameters become validated
only when this calibration is run on real recordings — they
are not pre-baked physiological constants.
diff --git a/myogen/utils/fr_optimization_targets.py b/myogen/utils/fr_optimization_targets.py
new file mode 100644
index 00000000..2604396b
--- /dev/null
+++ b/myogen/utils/fr_optimization_targets.py
@@ -0,0 +1,362 @@
+"""
+Firing rate optimization targets for DD network tuning.
+
+Provides closed-form predictions for number of active motor units and
+mean firing rate at any %MVC level. Multiple target models are supported
+so the DD optimizer can be tuned against different physiological
+assumptions.
+
+Available models
+----------------
+- ``fuglevand``: Linear rate coding (Fuglevand et al. 1993) [fuglevand1993].
+ Simple, closed-form. Tends to underestimate real rate coding at
+ higher %MVC because it assumes linear FR vs excitation.
+
+- ``de_luca``: Exponential approach to peak (De Luca & Contessa 2012) [delucacontessa2012].
+ Captures rapid rise after recruitment (PIC effect) and plateau at
+ peak rate. Supports onion-skin (earliest recruited fire fastest).
+ Recommended default for new work.
+
+- ``experimental``: Lookup table from published decomposed MU data.
+ Muscle-specific, most accurate where data is available. Requires
+ a target muscle name.
+
+All models return the same dict format:
+
+ {mvc_pct: {"n_active": int, "mean_fr_hz": float,
+ "fr_range": (min, max)}}
+
+Example
+-------
+>>> from myogen.utils.fr_optimization_targets import compute_fr_targets
+>>> thresholds = profile.generate_thresholds(mu_count_method="fixed", n=200)
+>>> targets = compute_fr_targets(thresholds, model="de_luca")
+>>> targets[10]["mean_fr_hz"]
+11.2
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+
+# ─── Model 1: Fuglevand linear (back-compat) ─────────────────────
+
+
+def fuglevand_targets(
+ thresholds: np.ndarray,
+ mvc_levels: list[float] | None = None,
+ min_fr_hz: float = 8.0,
+ max_fr_hz: float = 17.0,
+) -> dict[float, dict]:
+ """Fuglevand 1993 [fuglevand1993] linear rate coding targets.
+
+ FR_i(E) = min_fr + (max_fr - min_fr) * (E - threshold_i)
+ clipped to [min_fr, max_fr].
+
+ Parameters
+ ----------
+ thresholds : ndarray
+ Normalized recruitment thresholds (0 to 1).
+ mvc_levels : list of float, optional
+ %MVC levels to compute targets for. Default [5, 10, 20, 30, 50, 100].
+ min_fr_hz : float
+ Onset firing rate for newly recruited MUs.
+ max_fr_hz : float
+ Mean firing rate at 100% MVC.
+ """
+ if mvc_levels is None:
+ mvc_levels = [5, 10, 20, 30, 50, 100]
+
+ thresholds = np.sort(thresholds)
+ gain = max_fr_hz - min_fr_hz
+
+ targets = {}
+ for mvc in mvc_levels:
+ excitation = mvc / 100.0
+ active_mask = thresholds <= excitation
+ n_active = int(np.sum(active_mask))
+
+ if n_active == 0:
+ targets[mvc] = {"n_active": 0, "mean_fr_hz": 0.0, "fr_range": (0.0, 0.0)}
+ continue
+
+ active_thresholds = thresholds[active_mask]
+ frs = min_fr_hz + gain * (excitation - active_thresholds)
+ frs = np.clip(frs, min_fr_hz, max_fr_hz)
+
+ targets[mvc] = {
+ "n_active": n_active,
+ "mean_fr_hz": float(np.mean(frs)),
+ "fr_range": (float(frs.min()), float(frs.max())),
+ }
+
+ return targets
+
+
+# ─── Model 2: De Luca & Contessa exponential ────────────────────
+
+
+def de_luca_targets(
+ thresholds: np.ndarray,
+ mvc_levels: list[float] | None = None,
+ fr_onset_hz: float = 8.0,
+ fr_peak_early_hz: float = 16.0,
+ fr_peak_late_hz: float = 12.0,
+ steepness: float = 3.0,
+) -> dict[float, dict]:
+ """De Luca & Contessa 2012 [delucacontessa2012] exponential rate coding targets.
+
+ Each MU's firing rate follows an exponential approach from
+ ``fr_onset_hz`` at recruitment to its peak rate at full excitation:
+
+ FR_i(E) = FR_peak_i - (FR_peak_i - fr_onset_hz)
+ * exp(-steepness * (E - threshold_i))
+
+ For E < threshold_i, the MU is silent.
+
+ Onion-skin rate coding: the first-recruited MU has a higher peak
+ rate (``fr_peak_early_hz``) than the last-recruited MU
+ (``fr_peak_late_hz``). Peak rates interpolate linearly between
+ these two endpoints by threshold rank.
+
+ This produces the characteristic nonlinear rate coding observed
+ in real muscles: rapid rise after recruitment (matching PIC
+ amplification in biophysical MN models), then saturation at peak
+ (matching Ca²⁺ saturation and AHP dynamics).
+
+ Parameters
+ ----------
+ thresholds : ndarray
+ Normalized recruitment thresholds (0 to 1).
+ mvc_levels : list of float, optional
+ %MVC levels to compute targets for. Default [5, 10, 20, 30, 50, 100].
+ fr_onset_hz : float
+ Firing rate at recruitment (both onset and minimum). Typical 6-8 Hz.
+ fr_peak_early_hz : float
+ Peak rate (at MVC) for the first-recruited (lowest threshold) MU.
+ Default 16 Hz tuned to approach Skarabot 2023 [skarabot2023] untrained VL mean
+ (14 Hz). For trained VL use ~24 Hz, for FDI ~26 Hz.
+ fr_peak_late_hz : float
+ Peak rate (at MVC) for the last-recruited (highest threshold) MU.
+ Typically 70-85% of fr_peak_early_hz (onion-skin principle).
+ Default 12 Hz pairs with fr_peak_early_hz=16 for untrained VL.
+ steepness : float
+ Exponential approach rate. Higher = faster saturation.
+ Typical 2-5 for human limb muscles.
+ """
+ if mvc_levels is None:
+ mvc_levels = [5, 10, 20, 30, 50, 100]
+
+ thresholds = np.sort(thresholds)
+ n_mu = len(thresholds)
+
+ # Per-MU peak rates: interpolate from early to late by rank
+ if n_mu > 1:
+ ranks = np.arange(n_mu) / (n_mu - 1) # 0 → 1
+ else:
+ ranks = np.array([0.0])
+ fr_peaks = fr_peak_early_hz + ranks * (fr_peak_late_hz - fr_peak_early_hz)
+
+ targets = {}
+ for mvc in mvc_levels:
+ excitation = mvc / 100.0
+ active_mask = thresholds <= excitation
+ n_active = int(np.sum(active_mask))
+
+ if n_active == 0:
+ targets[mvc] = {"n_active": 0, "mean_fr_hz": 0.0, "fr_range": (0.0, 0.0)}
+ continue
+
+ active_thresholds = thresholds[active_mask]
+ active_peaks = fr_peaks[active_mask]
+ # Exponential approach: FR = peak - (peak - onset) * exp(-k * (E - thr))
+ delta_e = excitation - active_thresholds
+ frs = active_peaks - (active_peaks - fr_onset_hz) * np.exp(-steepness * delta_e)
+ frs = np.clip(frs, fr_onset_hz, active_peaks)
+
+ targets[mvc] = {
+ "n_active": n_active,
+ "mean_fr_hz": float(np.mean(frs)),
+ "fr_range": (float(frs.min()), float(frs.max())),
+ }
+
+ return targets
+
+
+# ─── Model 3: Experimental lookup ──────────────────────────────
+
+
+# Per-muscle lookup: {mvc_pct: (mean_fr_hz, std_fr_hz)}
+#
+# Values are monotonically increasing toward the MVC peak rate from
+# Skarabot 2023 [skarabot2023]. Sub-MVC values approximate the rise to peak based on
+# De Luca-style exponential approach.
+#
+# TODO: refine with direct measurements from decomposed MU recordings
+# (MyoGen paper Fig 4, Skarabot raw data, etc.). Current values are
+# best-estimate placeholders.
+EXPERIMENTAL_FR_TARGETS: dict[str, dict[float, tuple[float, float]]] = {
+ "vl": {
+ # Vastus lateralis, untrained.
+ # Peak: Skarabot 2023 [skarabot2023] untrained VL (14.0 ± 4.5 pps at MVC).
+ # Sub-MVC: exponential approach from onset ~8 Hz to MVC peak.
+ 5: (8.5, 2.0),
+ 10: (9.5, 2.3),
+ 20: (11.0, 2.8),
+ 30: (12.0, 3.2),
+ 50: (13.0, 3.8),
+ 100: (14.0, 4.5), # Skarabot untrained VL mean [skarabot2023]
+ },
+ "vl_trained": {
+ # Resistance-trained VL. Peak: Skarabot 2023 [skarabot2023] (21.1 ± 4.1 pps).
+ 5: (9.0, 2.2),
+ 10: (11.0, 2.6),
+ 20: (14.0, 3.0),
+ 30: (16.5, 3.4),
+ 50: (19.0, 3.8),
+ 100: (21.1, 4.1), # Skarabot trained VL mean [skarabot2023]
+ },
+ "fdi": {
+ # First dorsal interosseous. Higher peak rates than VL
+ # (hand muscles tend to fire faster). Placeholder based on
+ # MyoGen paper Fig 4 FDI data.
+ 5: (9.0, 2.0),
+ 10: (11.0, 2.3),
+ 20: (14.0, 2.8),
+ 30: (17.0, 3.2),
+ 50: (20.0, 3.8),
+ 100: (24.0, 4.5),
+ },
+ "ta": {
+ # Tibialis anterior. Placeholder — fewer decomposition studies
+ # than VL, should refine with data from Desmedt and others.
+ 5: (9.0, 2.0),
+ 10: (10.5, 2.3),
+ 20: (12.5, 2.8),
+ 30: (14.0, 3.2),
+ 50: (16.0, 3.8),
+ 100: (19.0, 4.5),
+ },
+}
+
+
+def experimental_targets(
+ thresholds: np.ndarray,
+ mvc_levels: list[float] | None = None,
+ muscle_name: str = "vl",
+) -> dict[float, dict]:
+ """Experimental firing rate targets from published decomposed MU data.
+
+ Uses a muscle-specific lookup table keyed by %MVC. Since the
+ experimental data gives population-level mean and std, all active
+ MUs at a given level are assigned the same target (the optimizer
+ matches mean/std, not per-MU rates).
+
+ Parameters
+ ----------
+ thresholds : ndarray
+ Normalized recruitment thresholds (0 to 1). Used only to
+ determine n_active at each %MVC.
+ mvc_levels : list of float, optional
+ %MVC levels. Must be keys in the lookup table.
+ muscle_name : str
+ Key into EXPERIMENTAL_FR_TARGETS. Default "vl".
+
+ Raises
+ ------
+ KeyError
+ If muscle_name is not in the lookup table.
+ ValueError
+ If a requested mvc_level is not in the table for that muscle.
+ """
+ if mvc_levels is None:
+ mvc_levels = [5, 10, 20, 30, 50, 100]
+
+ if muscle_name not in EXPERIMENTAL_FR_TARGETS:
+ raise KeyError(
+ f"No experimental targets for muscle '{muscle_name}'. "
+ f"Available: {list(EXPERIMENTAL_FR_TARGETS.keys())}"
+ )
+
+ lookup = EXPERIMENTAL_FR_TARGETS[muscle_name]
+ thresholds = np.sort(thresholds)
+
+ targets = {}
+ for mvc in mvc_levels:
+ if mvc not in lookup:
+ raise ValueError(
+ f"No experimental data for {mvc}% MVC in muscle '{muscle_name}'. "
+ f"Available levels: {sorted(lookup.keys())}"
+ )
+ mean_fr, std_fr = lookup[mvc]
+
+ excitation = mvc / 100.0
+ n_active = int(np.sum(thresholds <= excitation))
+
+ targets[mvc] = {
+ "n_active": n_active,
+ "mean_fr_hz": float(mean_fr),
+ "fr_range": (float(mean_fr - std_fr), float(mean_fr + std_fr)),
+ }
+
+ return targets
+
+
+# ─── Dispatcher ────────────────────────────────────────────────
+
+
+def compute_fr_targets(
+ thresholds: np.ndarray,
+ mvc_levels: list[float] | None = None,
+ model: str = "de_luca",
+ muscle_name: str = "vl",
+ **model_params,
+) -> dict[float, dict]:
+ """Compute firing rate optimization targets using the selected model.
+
+ Parameters
+ ----------
+ thresholds : ndarray
+ Normalized recruitment thresholds (0 to 1).
+ mvc_levels : list of float, optional
+ %MVC levels to compute targets for.
+ model : {"fuglevand", "de_luca", "experimental"}
+ Target model to use.
+
+ - "fuglevand": linear rate coding (Fuglevand 1993) [fuglevand1993]. Params:
+ ``min_fr_hz``, ``max_fr_hz``.
+ - "de_luca": exponential approach to peak (De Luca & Contessa
+ 2012) [delucacontessa2012]. Params: ``fr_onset_hz``, ``fr_peak_early_hz``,
+ ``fr_peak_late_hz``, ``steepness``.
+ - "experimental": muscle-specific lookup from published data.
+ Requires ``muscle_name``.
+ muscle_name : str
+ Muscle identifier for experimental lookup. Ignored by other models.
+ **model_params
+ Additional parameters passed to the selected model function.
+
+ Returns
+ -------
+ dict
+ Mapping mvc_pct -> {"n_active": int, "mean_fr_hz": float,
+ "fr_range": (min, max)}
+
+ Raises
+ ------
+ ValueError
+ If ``model`` is not recognized.
+ """
+ if model == "fuglevand":
+ return fuglevand_targets(thresholds, mvc_levels=mvc_levels, **model_params)
+ elif model == "de_luca":
+ return de_luca_targets(thresholds, mvc_levels=mvc_levels, **model_params)
+ elif model == "experimental":
+ return experimental_targets(
+ thresholds, mvc_levels=mvc_levels, muscle_name=muscle_name, **model_params
+ )
+ else:
+ raise ValueError(
+ f"Unknown rate coding model '{model}'. "
+ f"Options: fuglevand, de_luca, experimental"
+ )
diff --git a/myogen/utils/parallel_budget.py b/myogen/utils/parallel_budget.py
new file mode 100644
index 00000000..c1f9a6d6
--- /dev/null
+++ b/myogen/utils/parallel_budget.py
@@ -0,0 +1,192 @@
+"""Estimate safe parallelism for NEURON-based DD optimization trials.
+
+Each Optuna trial spawns a full NEURON simulation with AlphaMN__Pool +
+DescendingDrive__Pool + Network. NEURON uses global state, so parallel
+trials must run in separate **processes** (not threads). Each process
+needs its own memory for:
+
+- AlphaMN pool: ~1.5 MB per MN (compartmental model + ion channels)
+- DescendingDrive pool: ~0.1 MB per DD neuron
+- Network connectivity: ~0.5 MB per MN (synaptic connections)
+- Spike recorders + drive signal: ~10 MB overhead
+- Python/NEURON runtime: ~200 MB base per process
+
+Usage::
+
+ from myogen.utils.parallel_budget import estimate_safe_njobs
+
+ n_jobs = estimate_safe_njobs(n_mu=300, dd_neurons=280)
+ study.optimize(objective, n_trials=100, n_jobs=n_jobs)
+"""
+
+import os
+import platform
+
+
+def _get_system_memory_gb() -> float:
+ """Return total system RAM in GB."""
+ try:
+ import psutil
+ return psutil.virtual_memory().total / (1024**3)
+ except ImportError:
+ pass
+
+ # Fallback for Windows without psutil
+ if platform.system() == "Windows":
+ try:
+ import ctypes
+ kernel32 = ctypes.windll.kernel32
+ c_ulonglong = ctypes.c_ulonglong
+ class MEMORYSTATUSEX(ctypes.Structure):
+ _fields_ = [
+ ("dwLength", ctypes.c_ulong),
+ ("dwMemoryLoad", ctypes.c_ulong),
+ ("ullTotalPhys", c_ulonglong),
+ ("ullAvailPhys", c_ulonglong),
+ ("ullTotalPageFile", c_ulonglong),
+ ("ullAvailPageFile", c_ulonglong),
+ ("ullTotalVirtual", c_ulonglong),
+ ("ullAvailVirtual", c_ulonglong),
+ ("ullAvailExtendedVirtual", c_ulonglong),
+ ]
+ stat = MEMORYSTATUSEX()
+ stat.dwLength = ctypes.sizeof(stat)
+ kernel32.GlobalMemoryStatusEx(ctypes.byref(stat))
+ return stat.ullTotalPhys / (1024**3)
+ except Exception:
+ pass
+
+ # Last resort: assume 16 GB
+ return 16.0
+
+
+def _get_available_memory_gb() -> float:
+ """Return available (free) RAM in GB."""
+ try:
+ import psutil
+ return psutil.virtual_memory().available / (1024**3)
+ except ImportError:
+ pass
+
+ if platform.system() == "Windows":
+ try:
+ import ctypes
+ kernel32 = ctypes.windll.kernel32
+ c_ulonglong = ctypes.c_ulonglong
+ class MEMORYSTATUSEX(ctypes.Structure):
+ _fields_ = [
+ ("dwLength", ctypes.c_ulong),
+ ("dwMemoryLoad", ctypes.c_ulong),
+ ("ullTotalPhys", c_ulonglong),
+ ("ullAvailPhys", c_ulonglong),
+ ("ullTotalPageFile", c_ulonglong),
+ ("ullAvailPageFile", c_ulonglong),
+ ("ullTotalVirtual", c_ulonglong),
+ ("ullAvailVirtual", c_ulonglong),
+ ("ullAvailExtendedVirtual", c_ulonglong),
+ ]
+ stat = MEMORYSTATUSEX()
+ stat.dwLength = ctypes.sizeof(stat)
+ kernel32.GlobalMemoryStatusEx(ctypes.byref(stat))
+ return stat.ullAvailPhys / (1024**3)
+ except Exception:
+ pass
+
+ return 8.0
+
+
+def estimate_trial_memory_gb(
+ n_mu: int = 300,
+ dd_neurons: int = 280,
+) -> float:
+ """Estimate memory footprint of one DD optimization trial in GB.
+
+ Based on empirical measurements of NEURON AlphaMN pool simulations.
+
+ Parameters
+ ----------
+ n_mu : int
+ Number of motor units (AlphaMN cells).
+ dd_neurons : int
+ Number of descending drive neurons.
+
+ Returns
+ -------
+ float
+ Estimated memory per trial process in GB.
+ """
+ base_process_gb = 0.25 # Python + NEURON runtime overhead
+ per_mn_mb = 1.5 # compartmental model + channels + synapses
+ per_dd_mb = 0.1 # DD gamma cell (lightweight)
+ overhead_mb = 50 # spike recorders, drive signal, network state
+
+ trial_mb = (
+ base_process_gb * 1024
+ + n_mu * per_mn_mb
+ + dd_neurons * per_dd_mb
+ + overhead_mb
+ )
+ return trial_mb / 1024
+
+
+def estimate_safe_njobs(
+ n_mu: int = 300,
+ dd_neurons: int = 280,
+ reserve_gb: float = 4.0,
+ max_jobs: int | None = None,
+ verbose: bool = True,
+) -> int:
+ """Estimate the number of parallel Optuna trials that can run safely.
+
+ Checks available RAM and CPU cores, estimates per-trial memory,
+ and returns a conservative ``n_jobs`` value that leaves headroom
+ for the OS and other applications.
+
+ Parameters
+ ----------
+ n_mu : int
+ Motor units per trial.
+ dd_neurons : int
+ DD neurons per trial (used in memory estimate).
+ reserve_gb : float
+ RAM to reserve for OS + other applications (default 4 GB).
+ max_jobs : int or None
+ Hard cap on parallelism. None = no cap beyond hardware.
+ verbose : bool
+ Print the estimation breakdown.
+
+ Returns
+ -------
+ int
+ Recommended ``n_jobs`` for ``study.optimize()``. Always >= 1.
+ """
+ total_ram = _get_system_memory_gb()
+ available_ram = _get_available_memory_gb()
+ cpu_cores = os.cpu_count() or 4
+ trial_mem = estimate_trial_memory_gb(n_mu, dd_neurons)
+
+ # Use available RAM minus reserve, not total
+ usable_ram = max(available_ram - reserve_gb, trial_mem)
+ n_by_ram = int(usable_ram / trial_mem)
+
+ # Leave at least 1 core for OS + responsiveness
+ n_by_cpu = max(cpu_cores - 1, 1)
+
+ n_jobs = max(1, min(n_by_ram, n_by_cpu))
+
+ if max_jobs is not None:
+ n_jobs = min(n_jobs, max_jobs)
+
+ if verbose:
+ print(f" Parallel budget estimation:")
+ print(f" System RAM: {total_ram:.1f} GB total, "
+ f"{available_ram:.1f} GB available")
+ print(f" CPU cores: {cpu_cores}")
+ print(f" Per-trial memory: ~{trial_mem:.2f} GB "
+ f"({n_mu} MNs + {dd_neurons} DD neurons)")
+ print(f" Reserve: {reserve_gb:.1f} GB")
+ print(f" Limit by RAM: {n_by_ram} trials")
+ print(f" Limit by CPU: {n_by_cpu} trials")
+ print(f" => Recommended n_jobs = {n_jobs}")
+
+ return n_jobs
diff --git a/myogen/utils/ramp_analysis.py b/myogen/utils/ramp_analysis.py
new file mode 100644
index 00000000..8f43c103
--- /dev/null
+++ b/myogen/utils/ramp_analysis.py
@@ -0,0 +1,179 @@
+# myogen/utils/ramp_analysis.py
+"""
+Extract recruitment curves from trapezoidal ramp simulations.
+
+Analyzes spike trains produced during a time-varying drive signal
+to count active MUs and compute mean firing rates at %MVC checkpoints.
+"""
+
+import numpy as np
+from neo import Block
+
+
+def extract_ramp_recruitment(
+ spike_block: Block,
+ drive_signal: np.ndarray,
+ threshold_hz: float,
+ max_drive_hz: float,
+ timestep_ms: float = 0.1,
+ checkpoints_pct: list[float] | None = None,
+ window_ms: float = 500.0,
+ min_spikes: int = 3,
+) -> dict[float, dict]:
+ """Extract recruitment data at %MVC checkpoints during ramp.
+
+ Parameters
+ ----------
+ spike_block : neo.Block
+ Block with one segment containing spike trains from the ramp sim.
+ drive_signal : ndarray
+ The time-varying drive signal (Hz) used in the simulation.
+ threshold_hz : float
+ DD recruitment threshold (Hz).
+ max_drive_hz : float
+ Maximum DD drive (Hz) corresponding to 100% MVC.
+ timestep_ms : float
+ Simulation timestep in ms.
+ checkpoints_pct : list of float
+ %MVC levels to extract data at. Default [5, 10, 20, 30].
+ window_ms : float
+ Sliding window width for counting spikes (ms).
+ min_spikes : int
+ Minimum spikes in window to count a MU as "active".
+
+ Returns
+ -------
+ dict
+ Mapping mvc_pct -> {"n_active": int, "mean_fr_hz": float}
+ """
+ if checkpoints_pct is None:
+ checkpoints_pct = [5, 10, 20, 30]
+
+ spiketrains = spike_block.segments[0].spiketrains
+ window_s = window_ms / 1000.0
+
+ def mvc_to_hz(pct):
+ return threshold_hz + (pct / 100.0) * (max_drive_hz - threshold_hz)
+
+ results = {}
+ for mvc_pct in checkpoints_pct:
+ target_hz = mvc_to_hz(mvc_pct)
+
+ # Find the time window where drive is closest to target_hz
+ # (during the ramp-up phase — first occurrence)
+ drive_diff = np.abs(drive_signal - target_hz)
+ midpoint = len(drive_signal) // 2
+ ramp_up_diff = drive_diff[:midpoint]
+ best_step = int(np.argmin(ramp_up_diff))
+ center_s = best_step * timestep_ms / 1000.0
+
+ win_start = center_s - window_s / 2
+ win_end = center_s + window_s / 2
+
+ n_active = 0
+ firing_rates = []
+ for st in spiketrains:
+ spike_times = np.asarray(st.magnitude)
+ in_window = spike_times[(spike_times >= win_start) & (spike_times < win_end)]
+ if len(in_window) >= min_spikes:
+ n_active += 1
+ if len(in_window) >= 2:
+ isis = np.diff(in_window)
+ fr = 1.0 / np.mean(isis)
+ firing_rates.append(fr)
+
+ mean_fr = float(np.mean(firing_rates)) if firing_rates else 0.0
+ results[mvc_pct] = {"n_active": n_active, "mean_fr_hz": mean_fr}
+
+ return results
+
+
+def recruitment_force_pct_mvc(
+ spike_times_s,
+ force_pct_mvc: np.ndarray,
+ force_fs_hz: float,
+ window_s: float = 0.5,
+ min_spikes: int = 3,
+) -> np.ndarray:
+ """Recruitment force of each motor unit on an MVC ramp, in % MVC.
+
+ A unit is recruited at its first spike that starts ``min_spikes`` spikes
+ within ``window_s`` (the rule of ``extract_ramp_recruitment``); its
+ recruitment force is the force at that spike. Units that never meet the
+ rule get NaN.
+
+ Parameters
+ ----------
+ spike_times_s : sequence of array_like
+ Spike times (s) of each unit, in pool order.
+ force_pct_mvc : ndarray
+ Force during the ramp in % MVC (e.g. ``ForceModel`` output divided by
+ its MVC), sampled at ``force_fs_hz`` from t = 0.
+ force_fs_hz : float
+ Sampling rate of ``force_pct_mvc`` (Hz).
+ """
+ force = np.asarray(force_pct_mvc, dtype=float)
+ out = np.full(len(spike_times_s), np.nan)
+ for unit, times in enumerate(spike_times_s):
+ t = np.sort(np.asarray(times, dtype=float))
+ n = min_spikes - 1
+ if len(t) <= n:
+ continue
+ sustained = np.flatnonzero(t[n:] - t[:-n] <= window_s)
+ if len(sustained):
+ index = min(int(round(t[sustained[0]] * force_fs_hz)), len(force) - 1)
+ out[unit] = force[index]
+ return out
+
+
+def cv_from_recruitment_force(
+ recruitment_pct_mvc: np.ndarray,
+ low_band_pct: tuple[float, float] = (0.0, 30.0),
+ high_band_pct: tuple[float, float] = (50.0, 70.0),
+ cv_low__m_s: float = 3.89,
+ cv_high__m_s: float = 5.12,
+ unrecruited_pct_mvc: float = 100.0,
+ residual_sd__m_s: float = 0.39,
+ rng: np.random.Generator | None = None,
+) -> tuple[np.ndarray, tuple[float, float]]:
+ """Motor-unit conduction velocity linear in recruitment force.
+
+ Del Vecchio et al. (2018) measured CV rising linearly with recruitment
+ threshold in tibialis anterior with high-density surface EMG: 3.89 m/s for
+ units recruited at 0-30 % MVC and 5.12 m/s at 50-70 % MVC
+ [delvecchio2018 p. 2]. The line ``cv = a + b * rt`` passes through those
+ values at the mean recruitment force of this pool's units in each band
+ (a derivation: the paper gives no band-mean thresholds). Units with NaN
+ recruitment force (not recruited on the ramp) are placed at
+ ``unrecruited_pct_mvc``.
+
+ Each unit then scatters about the line by ``N(0, residual_sd__m_s)``. The
+ default 0.39 m/s is derived, not measured: CV is linear in threshold with
+ R² 0.70 ± 0.09 and SD 0.71 m/s over units [delvecchio2018 p. 2], so
+ sqrt(1 - 0.70) * 0.71 of it is left unexplained. Pass 0 for the line alone.
+ ``rng`` defaults to myogen's global generator (``set_random_seed``).
+
+ Returns
+ -------
+ cv : ndarray
+ CV of each unit (m/s), e.g. for ``SurfaceFiberDispersion(fiber_cv_per_unit__m_s=...)``.
+ (a, b) : tuple of float
+ Intercept (m/s) and slope (m/s per % MVC).
+ """
+ rt = np.asarray(recruitment_pct_mvc, dtype=float)
+ rt = np.where(np.isnan(rt), unrecruited_pct_mvc, rt)
+ low = (rt >= low_band_pct[0]) & (rt < low_band_pct[1])
+ high = (rt >= high_band_pct[0]) & (rt <= high_band_pct[1])
+ if not low.any() or not high.any():
+ raise ValueError("cv_from_recruitment_force needs recruited units in both bands")
+ x_low, x_high = rt[low].mean(), rt[high].mean()
+ b = (cv_high__m_s - cv_low__m_s) / (x_high - x_low)
+ a = cv_low__m_s - b * x_low
+ cv = a + b * rt
+ if residual_sd__m_s > 0:
+ if rng is None:
+ from myogen import get_random_generator
+
+ rng = get_random_generator()
+ cv = cv + rng.normal(0.0, residual_sd__m_s, len(cv))
+ return cv, (float(a), float(b))
diff --git a/properdocs.yml b/properdocs.yml
index 3e881a3f..95bc47f0 100644
--- a/properdocs.yml
+++ b/properdocs.yml
@@ -122,6 +122,7 @@ nav:
- Home: index.md
- Getting Started: getting-started.md
- Neo blocks: neo-blocks.md
+ - Literature: literature.md
- Examples:
- Basics: auto_examples/01_basic/index.md
- Fine-tuning: auto_examples/02_finetune/index.md
diff --git a/scripts/build_docs.py b/scripts/build_docs.py
index fe835d34..7f0b0820 100644
--- a/scripts/build_docs.py
+++ b/scripts/build_docs.py
@@ -50,7 +50,7 @@
EXECUTE_SUBDIRS = ["01_basic", "02_finetune"]
# Source-only CLI example: requires arguments and trusted local input files.
# Keep in sync with filename_pattern in docs/gallery_conf.py.
-SOURCE_ONLY = {"16_adaptive_cpu_execution"}
+SOURCE_ONLY = {"16_adaptive_cpu_execution", "06_optimize_dd_multilevel", "07_fdi_spikes_force"}
# keep in sync with ignore_pattern in docs/gallery_conf.py
IGNORE = re.compile(
r"(14_calibrate_noise_from_real|_oscillating_dc_helpers|_optimize_dc_worker|_pic_protocols)\.py"
diff --git a/scripts/build_literature.py b/scripts/build_literature.py
new file mode 100644
index 00000000..2fa2df65
--- /dev/null
+++ b/scripts/build_literature.py
@@ -0,0 +1,485 @@
+"""Build docs/literature.md, the literature register, from the table below.
+
+Usage:
+ python scripts/build_literature.py [--nexus NOTE.md] # write docs/literature.md (and an Obsidian mirror)
+ python scripts/build_literature.py --check # exit 1 if the file is stale or code cites unknown keys
+
+Code cites a source by its key and a locator in a comment, e.g. ``# [deluca1982 p. 118 Table 1]``. The code is
+scanned for such citations (CITE_DIRS); each source lists where it is cited, and unknown keys are errors.
+"""
+import argparse
+import re
+import sys
+from pathlib import Path
+
+REPO = Path(__file__).resolve().parents[1]
+_ap = argparse.ArgumentParser(description=__doc__.splitlines()[0])
+_ap.add_argument("--nexus", type=Path, default=None, help="also write the Obsidian mirror note here")
+_ap.add_argument("--check", action="store_true", help="write nothing; fail if docs/literature.md is stale or a key is unknown")
+ARGS = _ap.parse_args()
+NEXUS = ARGS.nexus
+DATE = "2026-11-02"
+
+# Repo targets (relative to docs/)
+R = {
+ "semg": ("SurfaceEMG", "../myogen/simulator/core/emg/surface/surface_emg.py"),
+ "simfib": ("simulate_fiber", "../myogen/simulator/core/emg/surface/simulate_fiber.py"),
+ "muscle": ("Muscle", "../myogen/simulator/core/muscle/muscle.py"),
+ "disp": ("SurfaceFiberDispersion", "../myogen/simulator/core/emg/surface/fiber_dispersion.py"),
+ "nerve": ("nerve_branches", "../myogen/simulator/core/emg/nerve_branches.py"),
+ "iemg": ("IntramuscularEMG", "../myogen/simulator/core/emg/intramuscular/intramuscular_emg.py"),
+ "force": ("ForceModel", "../myogen/simulator/core/force/force_model.py"),
+ "frt": ("fr_optimization_targets", "../myogen/utils/fr_optimization_targets.py"),
+ "opt": ("DD optimiser", "../examples/02_finetune/06_optimize_dd_multilevel.py"),
+ "fsf": ("07_fdi_spikes_force", "../examples/02_finetune/07_fdi_spikes_force.py"),
+ "futil": ("force_utils", "../myogen/simulator/core/force/force_utils.py"),
+ "mnpool": ("AlphaMN__Pool", "../myogen/simulator/neuron/populations/motor_neurons.py"),
+ "mncfg": ("alpha_mn_default.yaml", "../myogen/config/alpha_mn_default.yaml"),
+ "roe": ("Roeleveld bench", "../benchmarks/semg_roeleveld1997/README.md"),
+ "dispb": ("dispersion bench", "../benchmarks/semg_dispersion/README.md"),
+ "t_amp": ("test_surface_amplitude", "../tests/test_surface_amplitude.py"),
+ "t_fib": ("test_muscle_fibre_properties", "../tests/test_muscle_fibre_properties.py"),
+ "t_disp": ("test_surface_dispersion", "../tests/test_surface_dispersion.py"),
+ "t_gfl": ("test_gfluctdv_noise", "../tests/test_gfluctdv_noise.py"),
+ "t_dd": ("test_dd_optimizer_network", "../tests/test_dd_optimizer_network.py"),
+ "cl": ("CHANGELOG", "../CHANGELOG.md"),
+}
+MAIN = "2026-10-31 sEMG amplitude fixes and fibre dispersion"
+BANDS = "2026-10-31 FDI step-4 literature bands (draft)"
+# Nexus anchors: (note, heading or None). Headings with ':' or backticks are linked at note level.
+N = {
+ "partA": (MAIN, None, "Part A"), "partB": (MAIN, None, "Part B"),
+ "values": (MAIN, "Values chosen and why", None),
+ "recheck": (MAIN, "Re-check on the default (FDI) cylinder", None),
+ "skinlit": (MAIN, "Literature on FDI skin, fat and skin conductivity (2026-10-31)", None),
+ "layers": (MAIN, "FDI with literature layer values (2026-10-31)", None),
+ "sigma": (MAIN, "Where does skin σ = 1 S/m come from? (search 2026-10-31)", None),
+ "fig6": (MAIN, "Roeleveld 1997 Fig. 6 check (exponent-free) and source strength (2026-10-31)", None),
+ "hwang": (MAIN, "Hwang et al. 2013 read in full (FDI fibre diameter)", None),
+ "polgar": (MAIN, "Polgar et al. 1973 read in full (FDI fibre diameter)", None),
+ "diam": (MAIN, "FDI fibre-diameter bracket, simulated (2026-10-31)", None),
+ "rates": (MAIN, None, "Firing rate: De Luca 1982 onion-skin rates"),
+ "pool": (MAIN, "Pool size (FDI, 57.9 µm, De Luca 1982 rates)", None),
+ "botcorr": (MAIN, None, "Correction: what the Botelho 2019 reference is"),
+ "rebase": (MAIN, "Rebase onto main and re-test (2026-11-01 afternoon)", None),
+ "mdf": (MAIN, "MDF gap (2026-11-01 evening)", None),
+ "iapdec": (MAIN, None, "Decision: IAP length"),
+ "late": (MAIN, "Late recruitment of the largest FDI units (2026-11-01 night)", None),
+ "ceilwhy": (MAIN, "Why the ceiling sits at ~70 % MVC (2026-11-02)", None),
+ "spinal": (MAIN, "Are spinal mechanisms behind the low high-force rates? (2026-11-02)", None),
+ "calib": (MAIN, "Calibrations (120 MUs, Fuglevand 50-fold thresholds, De Luca 1982 targets)", None),
+ "b1": (BANDS, "1. Single-MU surface MUAP peak-to-peak amplitude", None),
+ "b2": (BANDS, "2. MUAP duration", None), "b3": (BANDS, "3. Global sEMG amplitude", None),
+ "b4": (BANDS, "4. Spectrum", None), "b5": (BANDS, "5. Anatomy and physiology", None),
+}
+
+# key, section, citation, identifiers, status, use, supports, repo targets, nexus anchors
+S = []
+def src(*a):
+ S.append(dict(zip(("key", "sec", "cite", "ids", "status", "use", "what", "repo", "nexus"), a)))
+
+# ── A. Volume conductor, tissue, source strength, IAP ──
+src("farina2004", "A", "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", "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.",
+ ["semg", "simfib", "muscle", "t_amp", "cl"], ["partA", "sigma"])
+src("roeleveld1997", "A", "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", "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", "roe", "t_amp", "cl"], ["sigma", "fig6", "rebase"])
+src("roeleveld1997b", "A", "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.", ["roe"], ["sigma"])
+src("blok2002", "A", "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", "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"], ["sigma", "fig6"])
+src("lowery2002", "A", "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", "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", "semg"], ["sigma"])
+src("farinarainoldi1999", "A", "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", "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.",
+ [], ["sigma"])
+src("farinamerletti2001", "A", "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; PMID 11396594", "FT", "CONTEXT",
+ "Same conductivity ratios; IAP 96 mV/mm³ with λ = 1 / mm.", ["simfib"], ["sigma"])
+src("merlettifarina2016", "A", "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.", [], ["sigma"])
+src("botelho2019", "A", "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; 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", "semg", "cl"], ["botcorr", "mdf", "b3", "b4", "layers"])
+src("gabriel1996", "A", "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"], ["sigma"])
+src("itis", "A", "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",
+ "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"], ["skinlit", "sigma"])
+src("dideriksen2019", "A", "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'.", [], ["sigma"])
+src("hodgkinnakajima1972", "A", "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", "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.", ["semg"], ["fig6"])
+src("rosenfalck1969", "A", "Rosenfalck P (1969). Intra- and extracellular potential fields of active nerve and muscle fibres. Acta Physiol Scand Suppl 321.",
+ "PMID 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.",
+ ["simfib", "semg", "t_amp", "cl"], ["iapdec", "mdf"])
+src("dimitrova2006", "A", "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).", [], ["mdf"])
+src("dimitrov1998", "A", "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.", [], ["mdf"])
+src("rodriguezfalces2013", "A", "Rodriguez-Falces J et al. (2013). Med Biol Eng Comput 50:447 (IAP / dipole metrics; title unverified).",
+ "doi:10.1007/s11517-012-0879-7", "FT", "CONTEXT", "IAP duration metrics for the λ question.", [], ["mdf"])
+src("petersen2019", "A", "Petersen E, Rostalski P (2019). A comprehensive mathematical model of motor unit pool organization, surface electromyography, and force generation. Front Physiol 10:176.",
+ "doi:10.3389/fphys.2019.00176; PMC6418040", "FT", "CONTEXT", "Farina-model implementation with λ = 1.", [], ["mdf"])
+src("andreassenrosenfalck1980", "A", "Andreassen & Rosenfalck (1980), muscle conductivities for the intramuscular kernel (details from code comment, unverified).",
+ "", "LINK", "ADOPTED", "Stated source of the radial and longitudinal muscle conductivities of the intramuscular point-electrode kernel in bioelectric.py (same defaults in fiber_simulation.py). Until 0.13.0 the defaults were 63 / 330 labelled S/m (1000x tissue values); since 2026-11-02 they are 0.063 / 0.33 S/m, the values the code comment attributes to this paper. Not in Zotero. Cited in code; not checked.", [], [])
+src("huigen2002", "A", "Huigen E, Peper A, Grimbergen CA (2002). Investigation into the origin of the noise of surface electrodes. Med Biol Eng Comput 40:332–338.",
+ "doi:10.1007/BF02344216", "LINK", "ADOPTED", "With gondran1996, the stated basis of the 1/f^α baseline-drift form with α in 1.5–2.0 (default α 1.75, the midpoint) in IntramuscularEMG.add_noise and emg_noise.py. Not in Zotero. Cited in code; not checked.", [], [])
+src("gondran1996", "A", "Gondran C, Siebert E, Yacoub S, Novakov E (1996). Noise of surface bio-potential electrodes based on NASICON ceramic and Ag–AgCl. Med Biol Eng Comput 34:460–466.",
+ "doi:10.1007/BF02523851", "LINK", "ADOPTED", "With huigen2002, the stated basis of the 1/f^α baseline-drift form with α in 1.5–2.0 (default α 1.75) in IntramuscularEMG.add_noise and emg_noise.py. Not in Zotero. Cited in code; not checked.", [], [])
+src("deluca2010", "A", "De Luca CJ, Gilmore LD, Kuznetsov M, Roy SH (2010). Filtering the surface EMG signal: Movement artifact and baseline noise contamination. J Biomech 43(8):1573–1579.",
+ "doi:10.1016/j.jbiomech.2010.01.027; PMID 20206934", "LINK", "CONTEXT", "Cited in IntramuscularEMG.add_noise for broadband movement artifacts (0–20 Hz), which the baseline-drift model leaves out. Citation from Zotero (no full text attached). Cited in code; not checked.", [], [])
+
+# ── B. Fibre CV, end-plates, tendon ends, branch delays ──
+src("nandedkar1983", "B", "Nandedkar SD, Stålberg E (1983). Simulation of single muscle fibre action potentials. Med Biol Eng Comput 21:158–165.",
+ "doi:10.1007/BF02441531", "FT", "ADOPTED",
+ "Fibre CV = 0.95 + 0.05·d[µm] m/s (fixes main's 2.0 m/s for every fibre); CV–diameter line for diameter-scaled sources. Read 2026-11-02: a modelling construct, not a measured regression — Stålberg 1966 biceps diameters 25–85 µm and CV 3.7 ± 0.71 m/s, with 25 / 55 / 85 µm pinned to the 2-SD band 2.2 / 3.7 / 5.2 m/s; linear, square and square-root forms gave similar amplitudes (p. 163, Fig. 6). Also eq. 1 of delvecchio2018 (p. 7). At FDI values it is inconsistent: 57.9 µm (polgar1973) → 3.85 m/s, 4.2 m/s (nishizono1990) → 65 µm.",
+ ["muscle", "disp", "t_fib", "t_disp", "cl"], ["values"])
+src("nishizono1990", "B", "Nishizono H, Fujimoto T, Ohtake H, Miyashita M (1990). Muscle fiber conduction velocity and contractile properties estimated from surface electrode arrays. Electroencephalogr Clin Neurophysiol 75:75–81.",
+ "doi:10.1016/0013-4694(90)90154-C; PMID 1688776", "FT", "ADOPTED", "FDI motor-unit CV by spike-triggered surface arrays (5 mm IED): 50 units from 4 men, 3.2–5.0 m/s, mean 4.2 ± 0.5 (p. 77) — MyoGen's default mean. Units recruited at 120–930 g, i.e. 5–42 % of the 2.2 kg mean MVC (p. 78); skin temperature held at 31–32 °C (p. 77); CV correlated with threshold force, r = 0.52 (p. 79). Botelho 2019 takes its CV range from this paper.", ["muscle"], ["values", "b5"])
+src("kouyoumdjian2023", "B", "Kouyoumdjian JA et al. (2023). Muscle fiber conduction velocity in situ revisited. Front Neurol 14:1118510.",
+ "doi:10.3389/fneur.2023.1118510; PMC9996126", "FT", "CONTEXT", "TA single-fibre CV 4.10 ± 0.66 m/s.", [], ["values"])
+src("zalewska2012", "B", "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; PMC3514689", "FT", "CONTEXT", "Alternative CV–diameter line v = 2.2 + 0.05(d − 25).", [], ["values"])
+src("blijham2006", "B", "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", "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.", [], ["values"])
+src("delvecchio2018", "B", "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", "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.",
+ ["disp", "cl", "muscle"], ["values"])
+src("delvecchio2017", "B", "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", "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.",
+ [], ["values"])
+src("akhmadeev2020", "B", "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; code: 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", "iemg", "disp", "t_disp", "cl"], ["partB", "values"])
+src("stalbergkarlsson2001", "B", "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", "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).", [], ["values"])
+src("piasecki2020", "B", "Piasecki M et al. (2020/2021). Near-fibre EMG (medRxiv; Clin Neurophysiol 132:1089).",
+ "doi:10.1101/2020.10.07.20208348", "FT", "CONTEXT", "Healthy end-plate SD 1 mm (from Stålberg & Karlsson).", [], ["values"])
+src("lateva2010", "B", "Lateva ZC, McGill KC, Johanson ME (2010). J Appl Physiol 108:1530.",
+ "PMC2886675", "FT", "CONTEXT", "Myelinated branch ~50 m/s bounds fast branch delays; fibre-end spread.", [], ["values"])
+src("keenan2005", "B", "Keenan KG, Farina D et al. (2005). J Appl Physiol 98:120 (FDI sEMG simulation).",
+ "doi:10.1152/japplphysiol.00894.2004", "FT", "ADOPTED",
+ "Convention for the opt-in 5 mm tendon-end band (no measurement of fibre-end stagger found).", ["disp", "cl"], ["values"])
+src("farinafosci2002", "B", "Farina D, Fosci M, Merletti R (2002). J Appl Physiol 92:235.",
+ "doi:10.1152/jappl.2002.92.1.235", "FT", "CONTEXT", "IZ / tendon spread 10 mm used as the sensitivity bound.", [], ["values"])
+src("jahanmiri2015", "B", "Jahanmiri-Nezhad F et al. (2015). Front Hum Neurosci 9:239.",
+ "doi:10.3389/fnhum.2015.00239", "FT", "CONTEXT", "Single-MU IZ spread 9–18 mm.", [], ["values"])
+
+# ── C. FDI anatomy and layers ──
+src("jacobson1992", "C", "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", "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", "cl"], ["recheck"])
+src("li2018", "C", "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", "FT", "ADOPTED", "Skin thickness 1.2 mm: hand dorsum over MCP II–III, median of 31 controls (0.79–1.50).",
+ ["muscle", "t_amp", "cl"], ["skinlit", "b5"])
+src("nishikawa2024", "C", "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", "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", "t_amp", "cl"], ["skinlit", "b5"])
+src("storchle2018", "C", "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", "FT", "REJECTED", "Old fat default 0.3 mm (whole-hand average, excludes skin).", ["muscle"], ["skinlit"])
+src("brodar1960", "C", "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"], ["skinlit"])
+src("krammer2025", "C", "Krammer et al. (2025). Diagnostics 15:1600.", "", "FT", "CONTEXT", "Hand-dorsum dermis 0.52 ± 0.16 mm.", [], ["skinlit"])
+src("infantolino2010", "C", "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; PMC2849523", "FT", "CONTEXT", "Volume 5.2 mL, fascicle 29.9 mm; derived CSA ~89 mm² vs model 150 mm².", [], ["pool"])
+src("martin1988", "B", "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", "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.",
+ [], ["values"])
+src("andreassen1987", "B", "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", "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).",
+ [], ["values"])
+src("burketsairis1973", "C", "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", "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.",
+ [], ["values"])
+src("polgar1973", "C", "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", "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.", [], ["polgar", "diam"])
+src("garnett1979", "C", "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; 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"], [])
+src("hwang2013", "C", "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", "FT", "REJECTED",
+ "FDI ~31 µm (elderly, paraffin, diameter undefined): judged an underestimate. Source of Botelho's 24.8–40.6 µm.", [], ["hwang", "diam"])
+src("feinstein1955", "C", "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).", [], ["pool", "b5"])
+src("zong2021", "C", "Zong Y et al. (2021). IEEE Trans Neural Syst Rehabil Eng 29:1199.",
+ "doi:10.1109/TNSRE.2021.3088061; PMC8780215", "FT", "CONTEXT", "FDI MUNE 108 ± 21; CMAP 16.75 mV (ceiling).", [], ["b3", "b5"])
+src("brooke1970", "C", "Brooke MH, Kaiser KK (1970). Muscle fiber types: how many and what kind? Arch Neurol 23:369–379.",
+ "doi: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.", [], [])
+
+# ── D–G. Recruitment, motor-neuron model, force, EMG benchmarks (section letter per entry; titles in TOPICS) ──
+src("deluca1982", "D", "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", "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).",
+ ["opt", "frt"], ["rates", "calib"])
+src("delucacontessa2012", "D", "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.", ["frt", "opt"], ["rates"])
+src("skarabot2023", "D", "Škarabot J et al. (2023). J Physiol (FDI/VL motor-unit decomposition).",
+ "PMC10952962", "FT", "CONTEXT", "FDI recruitment 0.8–59 % MVF (ceiling); VL peak rates in fr_optimization_targets.", ["frt"], ["b5"])
+src("milnerbrown1973a", "F", "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", "FT", "CONTEXT", "Methods; twitch shape.", [], ["rebase"])
+src("milnerbrown1973b", "F", "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; 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).",
+ ["force"], ["rebase"])
+src("milnerbrown1975", "G", "Milner-Brown HS, Stein RB (1975). The relation between the surface electromyogram and muscular force. J Physiol 246:549–569.",
+ "doi:10.1113/jphysiol.1975.sp010904", "FT", "CONTEXT", "FDI surface MUAP 0.14–0.40 mV (belly–tendon); duration 1.6–2.2 ms.", [], ["b1", "b2"])
+src("fuglevand1999", "F", "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", "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).", ["force", "futil", "opt"], ["ceilwhy"])
+src("mcnulty2000", "F", "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", "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.", ["force", "futil", "opt"], ["ceilwhy"])
+src("macefield1996", "F", "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", "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).", ["force", "futil", "opt"], ["ceilwhy"])
+src("fuglevand1993", "F", "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", "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.",
+ ["force", "fsf", "opt"], ["rebase"])
+src("destexhe2001", "E", "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.", ["t_gfl", "opt", "cl"], [])
+src("heckmanbinder1991", "E", "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", "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.",
+ ["opt", "t_dd", "cl"], ["late"])
+src("heckmanbinder1993", "E", "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.", [], ["late"])
+src("beauchamp2023", "E", "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; PMID 36626825; 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.", ["mnpool", "mncfg"], ["spinal"])
+src("gorassini2002", "E", "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; PMID 11929906", "LINK", "CONTEXT", "Origin of the paired motor-unit ΔF method used in the spinal check 2. Not read.", [], ["spinal"])
+src("afsharipour2020", "E", "Afsharipour B, Manzur N, Duchcherer J, Fenrich KF, Thompson CK, Negro F, et al. (2020). Estimation of self-sustained activity produced by persistent inward currents using firing rate profiles of multiple motor units in humans. J Neurophysiol 124:63–85.",
+ "doi:10.1152/jn.00194.2020; PMID 32459555", "LINK", "CONTEXT", "ΔF pairing criteria (control/test timing, rate-rate correlation); check 2 uses the common criteria from memory, not checked against this paper.", [], ["spinal"])
+src("powersheckman2017", "E", "Powers RK, Heckman CJ (2017). Synaptic control of the shape of the motoneuron pool input-output function. J Neurophysiol 117:1171–1184.",
+ "doi:10.1152/jn.00850.2016; PMID 28053245; PMC5340877; https://modeldb.science/239582", "ABS", "CONTEXT",
+ "Cat MG model pool: PICs raise gain and make pool rate a nonlinear function of excitation; neuromodulation and balanced inhibition increase the nonlinearity, push-pull inhibition reduces it; later-recruited units give higher force. Origin of MyoGen's Powers2017 cell model: its mechanisms na3rp, naps, kdrRL, mAHP, gh, L_Ca_inact and Gfluctdv are those of ModelDB 239582, this paper's model (2026-11-02 diff: whitespace, a seed RANGE and comment lines). Cell parameter values not compared.", ["mnpool", "mncfg"], ["spinal"])
+src("revillfuglevand2011", "E", "Revill AL, Fuglevand AJ (2011). Effects of persistent inward currents, accommodation, and adaptation on motor unit behavior: a simulation study. J Neurophysiol 106:1467–1479.",
+ "doi:10.1152/jn.00419.2011; PMID 21697447", "LINK", "CONTEXT", "PIC simulation of rate saturation and onion skin. Not read.", [], ["spinal"])
+src("revillfuglevand2017", "E", "Revill AL, Fuglevand AJ (2017). Inhibition linearizes firing rate responses in human motor units: implications for the role of persistent inward currents. J Physiol 595:179–191.",
+ "doi:10.1113/JP272823; PMID 27470946", "ABS", "CONTEXT", "Human MU rates rise steeply at recruitment then level off at modest rates while force keeps rising; reflex inhibition removes the steep rise (attributed to PICs).", [], ["spinal"])
+src("senneff2019", "E", "Senneff S, McManus L, Lowery MM (2019). Investigating the effect of persistent inward currents on motor unit firing rates and beta-band coherence in a model of the first dorsal interosseous muscle. Conf Proc IEEE EMBC 2019:2293–2296.",
+ "doi:10.1109/EMBC.2019.8857534; PMID 31946358", "ABS", "CONTEXT", "Prior FDI motoneuron-pool model with PICs at 20–40 % MVC: PICs shift mean rates and reduce beta-band coherence. Prior art for the spinal checks.", [], ["spinal"])
+src("powersbinder2001", "E", "Powers RK, Binder MD (2001). Input-output functions of mammalian motoneurons. Rev Physiol Biochem Pharmacol 143:137–263.",
+ "doi:10.1007/BFb0115594", "LINK", "CONTEXT",
+ "Review of motoneuron input-output properties; cited in the AlphaMN cell docstring. Not read here.", [], [])
+src("bui2003", "E", "Bui et al. (2003). Morphological and electrophysiological characterization of Ia inhibitory interneurons in cat spinal cord (details from code comment, unverified).",
+ "", "LINK", "ADOPTED", "Membrane area (81,390 µm²) and soma geometry of the INgII interneuron model in cells.py. Cited in code; not checked.", [], [])
+src("sanches2018", "E", "Sanches et al. (2018), synaptic parameter measurements (details from code comment, unverified).",
+ "", "LINK", "ADOPTED", "Default Exp2Syn synapse parameters (tau1 0.2 ms, tau2 0.3 ms, weight 0.6) in cells.py. Cited in code; not checked.", [], [])
+src("segev1990", "E", "Segev et al. (1990), alpha-function synaptic modelling (details from code comment, unverified).",
+ "", "LINK", "CONTEXT", "Listed as the reference for the synapse model form in cells.py (the code uses a double-exponential Exp2Syn). Cited in code; not checked.", [], [])
+src("enokaduchateau2017", "D", "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", "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.", [], [])
+src("dideriksen2012", "D", "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", "LINK", "CONTEXT",
+ "Simulation of synaptic noise and force steadiness; cited with enokaduchateau2017 for the ISI CV target. Not read here.", [], [])
+src("moritz2005", "D", "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**.", [], [])
+src("hu2013", "G", "Hu X, Rymer WZ, Suresh NL (2013). J Neurophysiol 110:1205; companion J NeuroEng Rehabil (2013).",
+ "PMC4073930; PMC3849018", "FT + FIG", "CONTEXT", "FDI decomposition; Delsys 5-pin geometry; MUAP amplitude linear in threshold.", [], ["b1", "b2"])
+src("mcmanus2015", "G", "McManus L et al. (2015). J Neurophysiol 113:3186.", "doi:10.1152/jn.00146.2015", "FT", "CONTEXT", "FDI surface MDF 166.9 Hz.", [], ["b4"])
+src("mcmanus2019", "G", "McManus L, Flood MW, Lowery MM (2019). J Neurophysiol.", "doi:10.1152/jn.00228.2019; PMC6766730", "FT + FIG", "CONTEXT",
+ "MDF ~120–150 Hz (400 Hz low-pass). 'MUAP amplitude' rejected as a target (metric undefined).", [], ["b1", "b4"])
+
+TOPICS = {"A": "Volume conductor and recording (tissue layers, conductivities, source strength, IAP)",
+ "B": "Muscle fibres and innervation (conduction velocity, end-plates, tendon ends, branch delays)",
+ "C": "Muscle anatomy (FDI dimensions, layers, motor-unit counts)",
+ "D": "Recruitment and rate coding (thresholds, firing rates, recruitment ceiling)",
+ "E": "Motor-neuron model and spinal mechanics (synaptic input, PICs, inhibition, noise)",
+ "F": "Muscle force (twitch, contraction time, force-frequency)",
+ "G": "Experimental EMG benchmarks (FDI amplitude and spectra)"}
+
+INDEX = [
+ ("sEMG volume conductor", ["farina2004"], ["semg", "simfib"]),
+ ("Skin σ = 1 S/m (effective, fitted)", ["roeleveld1997", "blok2002", "farinarainoldi1999", "farinamerletti2001", "farina2004", "lowery2002", "gabriel1996", "itis", "merlettifarina2016"], ["muscle", "roe"]),
+ ("Skin thickness 1.2 mm", ["li2018", "brodar1960", "krammer2025"], ["muscle"]),
+ ("Fat thickness 1.3 mm", ["nishikawa2024", "storchle2018"], ["muscle"]),
+ ("Muscle and fat σ", ["botelho2019", "merlettifarina2016"], ["muscle"]),
+ ("Source strength σi·π(d/2)²", ["lowery2002", "blok2002", "hodgkinnakajima1972"], ["semg"]),
+ ("IAP shape and λ (iap_lambda__per_mm)", ["rosenfalck1969", "farinamerletti2001", "botelho2019", "dimitrova2006", "dimitrov1998", "rodriguezfalces2013", "petersen2019"], ["simfib", "semg"]),
+ ("Fibre CV vs diameter", ["nandedkar1983", "blijham2006", "zalewska2012"], ["muscle", "disp"]),
+ ("Mean fibre CV 4.2 m/s", ["nishizono1990", "kouyoumdjian2023"], ["muscle"]),
+ ("MU CV vs threshold; within-MU CV spread", ["delvecchio2018", "delvecchio2017"], ["disp"]),
+ ("End-plate spread (two-layer model)", ["akhmadeev2020", "stalbergkarlsson2001", "piasecki2020", "jahanmiri2015"], ["nerve", "iemg", "disp"]),
+ ("Nerve-branch delays", ["akhmadeev2020", "lateva2010"], ["nerve"]),
+ ("Tendon-end band", ["keenan2005", "farinafosci2002", "lateva2010"], ["disp"]),
+ ("FDI fibre length and spread", ["jacobson1992"], ["muscle"]),
+ ("FDI fibre diameter", ["polgar1973", "hwang2013"], []),
+ ("FDI MU number, fibres, CSA", ["feinstein1955", "zong2021", "infantolino2010"], []),
+ ("FDI recruitment ceiling (--last-recruit-mvc)", ["deluca1982", "skarabot2023", "milnerbrown1973b", "macefield1996"], ["opt"]),
+ ("Rate-coding targets", ["deluca1982", "delucacontessa2012"], ["frt", "opt"]),
+ ("Twitch and force (ForceModel)", ["fuglevand1993", "milnerbrown1973a", "milnerbrown1973b", "fuglevand1999", "mcnulty2000", "macefield1996"], ["force", "fsf"]),
+ ("Human-unit force-frequency curve (opt-in)", ["macefield1996", "fuglevand1999", "mcnulty2000"], ["futil", "force", "opt"]),
+ ("High-force firing rates (--high-force-rates)", ["deluca1982"], ["opt"]),
+ ("Gfluctdv noise, ISI CV", ["destexhe2001", "moritz2005"], ["opt"]),
+ ("DD input weighting, equal fan-in", ["heckmanbinder1991", "heckmanbinder1993"], ["opt"]),
+ ("Motoneuron PICs, ΔF, rate patterns (diagnostics 2026-11-02)", ["beauchamp2023", "gorassini2002", "afsharipour2020", "powersheckman2017", "revillfuglevand2011", "revillfuglevand2017", "senneff2019", "deluca1982"], ["mnpool", "mncfg"]),
+ ("FDI sEMG comparators (RMS, MDF, MUAP)", ["botelho2019", "mcmanus2015", "mcmanus2019", "milnerbrown1975", "hu2013", "zong2021"], []),
+]
+
+CORRECTIONS = [
+ "**Skin σ = 1 S/m rationale.** First explained as gel electrodes bypassing the stratum corneum (Farina 2004). Unsourced and withdrawn: the value is Roeleveld et al. 1997's fitted effective skin, inherited by Blok 2002, Farina & Rainoldi 1999, Farina & Merletti 2001 and Farina 2004.",
+ "**Roeleveld 1997 DOI** is 10.1016/S1050-6411(97)00009-6 (an earlier (96)00034-8 was wrong).",
+ "**'Konstantin et al. 2020'** in MyoGen docstrings (muscle.py, intramuscular_emg.py, nerve_branches.py) is Akhmadeev et al. 2020. Fixed on `fix/semg-amplitude-scale` (028a14362); the `konstantin` threshold mode keeps its name.",
+ "**Botelho 2019 simulated amplitude**: median ~0.16 mV; 0.10 mV was one sample trace. Their 24.8–40.6 µm diameters are Hwang 2013's FDI mean ± SD.",
+ "**Feinstein 1955** authors are Feinstein, Lindegård, Nyman, Wohlfart; ~40,000 fibres (not 40,500).",
+ "**Piasecki 2020 and Stålberg & Karlsson 2001** are one citation chain, not two anchors.",
+ "**Jacobson 1992's 2.8 mm** is probably between-specimen spread; MyoGen uses it as a uniform half-width (docstring corrected).",
+ "**Unverified, quoted from memory:** Moritz 2005 (ISI CV), Kukulka & Clamann 1981. **Never read in full:** Nandedkar & Stålberg 1983 (adopted CV line), Jacobson 1992 Table 1, Heckman & Binder 1991.",
+ "**McManus 2019 'MUAP amplitude'**: metric undefined, not usable as a target.",
+ "**ISI CV target 15–22 %** in the optimiser cites enokaduchateau2017 and dideriksen2012; the range was not checked against either paper.",
+]
+
+STATUS = "FT full text · FIG read off a figure · ABS abstract only · SEC known only through a citing paper · CODE source code · LINK named only"
+USE = "ADOPTED sets a default or code value (or is its stated basis) · CONTEXT comparator, bound or background · REJECTED considered and not used"
+
+
+def idlinks(ids):
+ out = []
+ for part in [p.strip() for p in ids.split(";") if p.strip()]:
+ if part.startswith("doi:"):
+ d = part[4:]
+ out.append(f"[doi:{d}](https://doi.org/{d})")
+ elif part.startswith("PMC"):
+ out.append(f"[{part}](https://pmc.ncbi.nlm.nih.gov/articles/{part}/)")
+ elif part.startswith("PMID"):
+ out.append(f"[{part}](https://pubmed.ncbi.nlm.nih.gov/{part.split()[1]}/)")
+ elif part.startswith("http"):
+ out.append(f"[link]({part})")
+ elif part.startswith("code:"):
+ out.append(f"[code](https://{part[5:].strip()})")
+ else:
+ out.append(part)
+ return " · ".join(out)
+
+
+def repo_links(keys):
+ return ", ".join(f"[{R[k][0]}]({R[k][1]})" for k in keys) or "—"
+
+
+def nexus_links(keys, table=True):
+ out = []
+ bar = "\\|" if table else "|" # a bare | would split a table cell
+ for k in keys:
+ note, head, label = N[k]
+ out.append(f"[[{note}#{head}{bar}{head}]]" if head else f"[[{note}]] › {label}")
+ return "
".join(out) or "—"
+
+
+by_key = {s["key"]: s for s in S}
+
+# In-code citations: [key] or [key locator], key = letters + year (+ optional letter), e.g. [mcnulty2000 Table 2].
+CITE = re.compile(r"\[([a-z]+\d{4}[a-z]?)((?:[ ,;:][^\]\n]*)?)\]")
+CITE_DIRS = ["myogen", "examples/02_finetune", "benchmarks", "tests", "scripts"]
+CITE_SKIP = ("_archive", "__pycache__", "build_literature.py")
+
+
+def scan_code():
+ """{key: [(path relative to the repo, line, citation text)]} for every citation in CITE_DIRS."""
+ hits = {}
+ for d in CITE_DIRS:
+ for f in sorted((REPO / d).rglob("*")):
+ if f.suffix not in {".py", ".yaml", ".yml"} or any(x in f.parts or x == f.name for x in CITE_SKIP):
+ continue
+ for i, line in enumerate(f.read_text(encoding="utf-8", errors="replace").splitlines(), 1):
+ for m in CITE.finditer(line):
+ hits.setdefault(m.group(1), []).append((f.relative_to(REPO).as_posix(), i, m.group(0)))
+ return hits
+
+
+CITED = scan_code()
+UNKNOWN = {k: v for k, v in CITED.items() if k not in by_key}
+
+
+def cited_links(key):
+ return ", ".join(f"[{Path(p).name}:{i}](../{p}#L{i})" for p, i, _ in CITED.get(key, []))
+
+
+def render(nexus):
+ L = []
+ if nexus:
+ L += ["---", f"created: {DATE}", "tags:", " - semg", " - literature", " - register", "---", "",
+ f"> Mirror of the repo register `docs/literature.md` (branch `calib/fdi-dd-main`), generated {DATE} from one table. "
+ "The repo file is canonical for code links; this note adds links into the session notes. "
+ f"Parent: [[README - Start Here]]. Full transcript sweep (≈120 sources incl. passing mentions and local copies): [[Literature sweep 2026-10-28 to 2026-11-01]]. "
+ "Generator: `scripts/build_literature.py` in the repo (edit its table, rerun with `--nexus`, commit `docs/literature.md`).", ""]
+ L += ["# Literature register", ""]
+ if not nexus:
+ L += ["Sources behind MyoGen's surface-EMG defaults, the fibre-dispersion options and the FDI descending-drive calibration "
+ f"(work of 2026-10-30 to {DATE}). Each parameter links to its sources; each source links to the code, tests and notes that use it. "
+ "Every parameter needs a literature basis or must be labelled a modelling assumption.", "",
+ "**Citing in code:** a sourced value carries an inline comment with the key and a locator, e.g. "
+ "`HIGH_FORCE_RATE_TARGETS = {40: 25.3, 80: 41.4} # [deluca1982 p. 118 Table 1]`. Give only locators you have read; "
+ "write \"no locator\" otherwise. Each key below resolves to the full citation and DOI.", "",
+ "To add a source: add it to the table in `scripts/build_literature.py`, link it from the parameter index, cite the key in the code, "
+ "and rerun the script (`--check` fails on unknown keys or a stale file). \"Cited at\" lists are generated from the code.", ""]
+ L += [f"**Status:** {STATUS} ", f"**Use:** {USE}", "", "## Parameter index", "",
+ "| Parameter or claim | Sources | " + ("Session notes" if nexus else "Code") + " |", "|---|---|---|"]
+ for name, keys, code in INDEX:
+ srcs = ", ".join((f"[[#{by_key[k]['key']}]]" if nexus else f"[{k}](#{k})") for k in keys)
+ if nexus:
+ notes = sorted({n for k in keys for n in by_key[k]["nexus"]}, key=list(N).index)
+ col = nexus_links(notes)
+ else:
+ col = repo_links(code)
+ L.append(f"| {name} | {srcs} | {col} |")
+ for sec, title in TOPICS.items():
+ L += ["", f"## {sec}. {title}", ""]
+ if nexus:
+ for s in [x for x in S if x["sec"] == sec]:
+ L += [f"### {s['key']}", f"{s['cite']} {idlinks(s['ids'])}", "",
+ f"- **{s['use']}**, {s['status']}: {s['what']}",
+ f"- Notes: {nexus_links(s['nexus'], table=False).replace('
', '; ')}",
+ f"- Repo: {', '.join(R[k][0] for k in s['repo']) or '—'}", ""]
+ else:
+ L += ["| Source | Status | Use | What it supports | Used in |", "|---|---|---|---|---|"]
+ for s in [x for x in S if x["sec"] == sec]:
+ ids = idlinks(s["ids"])
+ used = repo_links(s["repo"]) + (f"
Cited at: {cited_links(s['key'])}" if s["key"] in CITED else "")
+ L.append(f"| **{s['key']}**: {s['cite']}{' ' + ids if ids else ''} | {s['status']} | {s['use']} | {s['what']} | {used} |")
+ L += ["", "## Corrections and open verification", ""] + [f"{i}. {c}" for i, c in enumerate(CORRECTIONS, 1)]
+ if nexus:
+ L += ["", "## Local copies", "",
+ "PDF/text copies (git-ignored), relative to the repository root: `sources/literature/2026-10_semg_fdi/lit/` "
+ "(subfolders fdi/, fdi_layers/, fdi_pool/, iap/, mb73/, skin_sigma/, source_strength/, iemg_simulator/), with hwang2013/, polgar1973/ and roeleveld/ next to lit/. "
+ "Inventory: [[Literature sweep 2026-10-28 to 2026-11-01]] §G."]
+ return "\n".join(L) + "\n"
+
+
+OUT = REPO / "docs" / "literature.md"
+missing = [k for _, ks, _ in INDEX for k in ks if k not in by_key]
+for k, v in UNKNOWN.items():
+ for p, i, c in v:
+ print(f"unknown citation key {k!r} at {p}:{i}: {c}", file=sys.stderr)
+n_cites = sum(len(v) for v in CITED.values())
+if ARGS.check:
+ stale = not OUT.exists() or OUT.read_text(encoding="utf-8") != render(False)
+ if stale:
+ print("docs/literature.md is stale: run python scripts/build_literature.py", file=sys.stderr)
+ sys.exit(1 if (stale or UNKNOWN or missing) else 0)
+OUT.write_text(render(False), encoding="utf-8")
+if NEXUS is not None:
+ NEXUS.write_text(render(True), encoding="utf-8")
+print(len(S), "sources;", len(INDEX), "index rows;", n_cites, "code citations of", len(CITED), "keys; missing index keys:",
+ missing, "; unknown code keys:", sorted(UNKNOWN))
+sys.exit(1 if (UNKNOWN or missing) else 0)
diff --git a/tests/test_dd_optimizer_network.py b/tests/test_dd_optimizer_network.py
new file mode 100644
index 00000000..fd2d6b14
--- /dev/null
+++ b/tests/test_dd_optimizer_network.py
@@ -0,0 +1,221 @@
+"""DD -> motor neuron network of the multilevel DD optimiser (examples/02_finetune/06_optimize_dd_multilevel.py).
+
+Every motor neuron gets the same number of DD inputs, and each input's weight is scaled so that it delivers the
+same current to a clamped soma whatever the neuron's size (equal effective synaptic current per input, the
+equal-share case of Heckman & Binder 1991). Parameter files without a ``network`` block are refused.
+"""
+
+import importlib.util
+from pathlib import Path
+
+import numpy as np
+import pytest
+import quantities as pq
+from neo import Block, Segment, SpikeTrain
+from neuron import h
+
+from myogen import load_nmodl_mechanisms, set_random_seed
+from myogen.simulator.neuron.populations import AlphaMN__Pool, DescendingDrive__Pool
+
+OPT_PATH = Path(__file__).resolve().parents[1] / "examples" / "02_finetune" / "06_optimize_dd_multilevel.py"
+
+
+@pytest.fixture(scope="module")
+def opt():
+ load_nmodl_mechanisms()
+ spec = importlib.util.spec_from_file_location("dd_multilevel_opt", OPT_PATH)
+ module = importlib.util.module_from_spec(spec)
+ spec.loader.exec_module(module)
+ return module
+
+
+def _pools(opt, n_dd=40):
+ set_random_seed(7)
+ thresholds = np.linspace(0.02, 1.0, 12)
+ mn = AlphaMN__Pool(recruitment_thresholds__array=thresholds, config_file="alpha_mn_default.yaml")
+ dd = DescendingDrive__Pool(n=n_dd, timestep__ms=opt.TIMESTEP_MS * pq.ms, process_type="gamma", shape=3.0)
+ return thresholds, mn, dd
+
+
+def _params(**network):
+ return {"dd_neurons": 40, "conn_probability": 0.3, "synaptic_weight": 0.05,
+ "network": {"fan_in": "equal", "weight_scaling": "dendritic_transfer", **network}}
+
+
+def _owner(mn):
+ owner = {}
+ for i, cell in enumerate(mn):
+ for sec in [cell.soma, *cell.dend]:
+ owner[sec.hname()] = i
+ return owner
+
+
+def test_every_motor_neuron_gets_the_same_number_of_excitatory_dd_inputs(opt):
+ _, mn, dd = _pools(opt)
+ net = opt.build_dd_mn_network(mn, dd, _params())
+ owner = _owner(mn)
+ counts = np.zeros(len(mn), int)
+ for nc in net.dd_to_mn:
+ syn = nc.syn()
+ assert syn.e > -40, "DD input landed on an inhibitory synapse"
+ counts[owner[syn.get_segment().sec.hname()]] += 1
+ assert np.all(counts == round(0.3 * 40)), counts
+
+
+def test_input_weights_equalise_the_current_reaching_the_soma(opt):
+ _, mn, dd = _pools(opt)
+ net = opt.build_dd_mn_network(mn, dd, _params())
+ k = np.asarray(net.transfer)
+ assert np.all(np.diff(k) < 0), "dendrite-to-soma transfer should fall with neuron size"
+ owner = _owner(mn)
+ w = np.zeros(len(mn))
+ for nc in net.dd_to_mn:
+ w[owner[nc.syn().get_segment().sec.hname()]] = nc.weight[0]
+ np.testing.assert_allclose(w * k, w[0] * k[0], rtol=1e-9)
+ assert w[0] == pytest.approx(0.05), "smallest unit keeps the nominal weight"
+
+
+def test_parameter_sets_without_a_network_block_are_refused(opt):
+ _, mn, dd = _pools(opt)
+ params = _params()
+ del params["network"]
+ with pytest.raises(ValueError, match="network"):
+ opt.build_dd_mn_network(mn, dd, params)
+
+
+def _weights(opt, net, mn):
+ owner = _owner(mn)
+ w = np.zeros(len(mn))
+ for nc in net.dd_to_mn:
+ w[owner[nc.syn().get_segment().sec.hname()]] = nc.weight[0]
+ return w
+
+
+@pytest.mark.parametrize("p", [0.0, 0.5])
+def test_weight_exponent_sets_how_much_larger_neurons_are_compensated(opt, p):
+ _, mn, dd = _pools(opt)
+ net = opt.build_dd_mn_network(mn, dd, _params(weight_exponent=p))
+ k = np.asarray(net.transfer)
+ w = _weights(opt, net, mn)
+ np.testing.assert_allclose(w, 0.05 * (k.max() / k) ** p, rtol=1e-9)
+
+
+def test_weight_exponent_defaults_to_equal_current_per_input(opt):
+ _, mn, dd = _pools(opt)
+ a = _weights(opt, opt.build_dd_mn_network(mn, dd, _params()), mn)
+ _, mn2, dd2 = _pools(opt)
+ b = _weights(opt, opt.build_dd_mn_network(mn2, dd2, _params(weight_exponent=1.0)), mn2)
+ np.testing.assert_allclose(a, b, rtol=1e-12)
+
+
+@pytest.mark.parametrize("bad", [-0.1, float("nan"), "1"])
+def test_invalid_weight_exponent_is_refused(opt, bad):
+ _, mn, dd = _pools(opt)
+ with pytest.raises(ValueError, match="weight_exponent must be"):
+ opt.build_dd_mn_network(mn, dd, _params(weight_exponent=bad))
+
+
+def test_network_block_records_the_weight_exponent(opt):
+ assert opt.network_block() == opt.NETWORK
+ assert opt.network_block(0.4) == {**opt.NETWORK, "weight_exponent": 0.4}
+
+
+def test_last_recruitment_excitation_uses_sustained_firing(opt):
+ dt_ms = 1.0
+ e = np.linspace(0.0, 1.0, 10001) # 10 s ramp, 1 ms steps
+ t_stop = 10.0 * pq.s
+ early = SpikeTrain(np.arange(1.0, 10.0, 0.05) * pq.s, t_stop=t_stop) # sustained from e = 0.1
+ late = SpikeTrain(np.r_[0.5, np.arange(6.0, 10.0, 0.05)] * pq.s, t_stop=t_stop) # one stray spike, then e = 0.6
+ block = Block()
+ block.segments.append(Segment())
+ block.segments[0].spiketrains = [early, late]
+ e_rec = opt.recruitment_excitation(block, e, timestep_ms=dt_ms)
+ np.testing.assert_allclose(e_rec, [0.1, 0.6], atol=0.03)
+
+ never = SpikeTrain(np.array([]) * pq.s, t_stop=t_stop)
+ block.segments[0].spiketrains = [early, never]
+ assert np.isnan(opt.recruitment_excitation(block, e, timestep_ms=dt_ms)[1])
+
+
+def test_recruitment_times_use_sustained_firing(opt):
+ t_stop = 10.0 * pq.s
+ early = SpikeTrain(np.arange(1.0, 10.0, 0.05) * pq.s, t_stop=t_stop)
+ late = SpikeTrain(np.r_[0.5, np.arange(6.0, 10.0, 0.05)] * pq.s, t_stop=t_stop)
+ never = SpikeTrain(np.array([]) * pq.s, t_stop=t_stop)
+ block = Block()
+ block.segments.append(Segment())
+ block.segments[0].spiketrains = [early, late, never]
+ t_rec = opt.recruitment_times(block)
+ np.testing.assert_allclose(t_rec[:2], [1.0, 6.0])
+ assert np.isnan(t_rec[2])
+
+
+def test_last_recruitment_force_is_read_from_the_smoothed_force_curve(opt):
+ fs = 1000.0
+ t = np.arange(0, 10.0, 1 / fs)
+ force = 50.0 * t # 0 -> 500 model units; MVC 500
+ pct, n_unrec = opt.last_recruitment_pct_mvc(np.array([1.0, 4.4, 2.0]), force, fs, mvc=500.0)
+ assert n_unrec == 0
+ assert pct == pytest.approx(100 * 50.0 * 4.4 / 500.0, abs=0.1)
+
+
+def test_unrecruited_units_push_the_ceiling_past_mvc(opt):
+ fs = 1000.0
+ force = np.full(10000, 100.0)
+ pct, n_unrec = opt.last_recruitment_pct_mvc(np.array([1.0, np.nan, np.nan, 2.0]), force, fs, mvc=100.0)
+ assert n_unrec == 2
+ assert pct == pytest.approx(150.0)
+
+
+def test_ceiling_target_is_given_in_percent_mvc(opt):
+ assert opt.ceiling_penalty(55.0, 55.0) == 0.0
+ assert opt.ceiling_penalty(66.0, 55.0) == pytest.approx(opt.CEILING_WEIGHT * 11.0 / 55.0)
+ assert opt.ceiling_penalty(80.0, None) == 0.0, "no target: the ceiling term is off"
+
+
+def test_ceiling_rank_scores_the_unit_with_k_units_above_it(opt):
+ fs = 1000.0
+ force = 10.0 * np.arange(0, 10.0, 1 / fs) # MVC 100 at t = 10 s
+ t_rec = np.array([1.0, 6.0, 3.0, 5.0])
+ assert opt.last_recruitment_pct_mvc(t_rec, force, fs, mvc=100.0)[0] == pytest.approx(60.0, abs=0.2)
+ assert opt.last_recruitment_pct_mvc(t_rec, force, fs, mvc=100.0, rank=1)[0] == pytest.approx(50.0, abs=0.2)
+
+
+def test_unrecruited_units_count_as_above_the_ranked_unit(opt):
+ fs = 1000.0
+ force = 10.0 * np.arange(0, 10.0, 1 / fs)
+ t_rec = np.array([1.0, np.nan, 3.0, 5.0])
+ pct, n_unrec = opt.last_recruitment_pct_mvc(t_rec, force, fs, mvc=100.0, rank=1)
+ assert n_unrec == 1
+ assert pct == pytest.approx(50.0, abs=0.2), "one unit (unrecruited) above: the highest recruited unit is scored"
+ pct, _ = opt.last_recruitment_pct_mvc(np.array([1.0, np.nan, np.nan, 5.0]), force, fs, mvc=100.0, rank=1)
+ assert pct == pytest.approx(125.0), "more unrecruited units than the rank allows: past 100 % as before"
+
+
+def test_rates_are_read_where_the_force_reaches_each_level(opt):
+ fs = 1000.0
+ force = 10.0 * np.arange(0, 10.0, 1 / fs) # 40 % MVC at 4 s, 80 % at 8 s
+ trains = [np.arange(0.0, 10.0, 1 / 20.0), np.r_[np.arange(0.0, 6.0, 1 / 10.0), np.arange(6.0, 10.0, 1 / 30.0)],
+ np.array([9.5])] # third unit not active at either level
+ rates = opt.rates_at_pct_mvc(trains, force, fs, mvc=100.0, levels=(40, 80))
+ assert rates[40] == pytest.approx((20 + 10) / 2, abs=1.0)
+ assert rates[80] == pytest.approx((20 + 30) / 2, abs=1.0)
+
+
+def test_high_force_rates_are_scored_against_de_luca_1982(opt):
+ assert opt.HIGH_FORCE_RATE_TARGETS == {40: 25.3, 80: 41.4}
+ assert opt.high_force_rate_penalty({40: 25.3, 80: 41.4}) == 0.0
+ expected = opt.HIGH_FORCE_RATE_WEIGHT * (5.06 / 25.3 + 8.28 / 41.4)
+ assert opt.high_force_rate_penalty({40: 20.24, 80: 33.12}) == pytest.approx(expected)
+
+
+def test_force_model_uses_the_selected_force_frequency_shape(opt):
+ from myogen.simulator.core.force.force_utils import FUGLEVAND_FORCE_FREQUENCY, HUMAN_UNIT_FORCE_FREQUENCY
+ thresholds = np.linspace(0.02, 1.0, 12)
+ assert opt.FORCE_FREQUENCY_SHAPE == FUGLEVAND_FORCE_FREQUENCY
+ old = opt.FORCE_FREQUENCY_SHAPE
+ try:
+ opt.FORCE_FREQUENCY_SHAPE = HUMAN_UNIT_FORCE_FREQUENCY
+ assert opt.make_force_model(thresholds).force_frequency_shape == HUMAN_UNIT_FORCE_FREQUENCY
+ finally:
+ opt.FORCE_FREQUENCY_SHAPE = old
diff --git a/tests/test_docs_gallery_execution.py b/tests/test_docs_gallery_execution.py
new file mode 100644
index 00000000..815a6550
--- /dev/null
+++ b/tests/test_docs_gallery_execution.py
@@ -0,0 +1,50 @@
+"""On pushes to main the docs build executes the examples that docs/gallery_conf.py selects.
+
+mkdocs-gallery runs them with exec(), where ``__file__`` is undefined, and without command-line arguments; a failed
+example fails the build (scripts/build_docs.py). Scripts that cannot run that way render source-only.
+"""
+import ast
+import re
+from pathlib import Path
+
+REPO = Path(__file__).resolve().parents[1]
+
+
+def _conf_value(key):
+ src = (REPO / "docs" / "gallery_conf.py").read_text(encoding="utf-8")
+ return ast.literal_eval(re.search(rf'"{key}":\s*(r"[^"]*")', src).group(1))
+
+
+def _executed_examples():
+ run, skip = re.compile(_conf_value("filename_pattern")), re.compile(_conf_value("ignore_pattern"))
+ for py in sorted((REPO / "examples").rglob("*.py")):
+ path = "/" + py.relative_to(REPO).as_posix()
+ if run.search(path) and not skip.search(path):
+ yield py
+
+
+def _unguarded_file_uses(py):
+ """Lines that use ``__file__`` outside a ``try`` whose handler catches NameError."""
+ tree = ast.parse(py.read_text(encoding="utf-8"))
+ guarded = set()
+ for node in ast.walk(tree):
+ if isinstance(node, ast.Try) and any(
+ h.type is not None and "NameError" in ast.unparse(h.type) for h in node.handlers
+ ):
+ for stmt in node.body:
+ guarded.update(id(n) for n in ast.walk(stmt))
+ return [n.lineno for n in ast.walk(tree)
+ if isinstance(n, ast.Name) and n.id == "__file__" and id(n) not in guarded]
+
+
+def test_executed_examples_do_not_need_file():
+ bad = {py.relative_to(REPO).as_posix(): lines for py in _executed_examples() if (lines := _unguarded_file_uses(py))}
+ assert bad == {}
+
+
+def test_source_only_lists_agree():
+ pattern = _conf_value("filename_pattern")
+ in_conf = set(re.findall(r"(\w+)\\.py\$", pattern))
+ src = (REPO / "scripts" / "build_docs.py").read_text(encoding="utf-8")
+ in_build = ast.literal_eval(re.search(r"^SOURCE_ONLY = (\{.*\})$", src, re.M).group(1))
+ assert in_conf == in_build
diff --git a/tests/test_force_model.py b/tests/test_force_model.py
index 9f35d662..22c3ffd4 100644
--- a/tests/test_force_model.py
+++ b/tests/test_force_model.py
@@ -6,7 +6,7 @@
from neo import Block, Segment, SpikeTrain
from myogen.simulator import ForceModel, ForceModelVectorized
-from myogen.simulator.core.force.force_utils import fuglevand_gain
+from myogen.simulator.core.force.force_utils import HUMAN_UNIT_FORCE_FREQUENCY, fuglevand_gain
MODELS = (ForceModel, ForceModelVectorized)
@@ -164,3 +164,47 @@ def test_force_single_spike_has_analytic_tail(model_class, unit):
lengths = np.ceil(10 * model.contraction_times__samples).astype(int)
assert model.twitch_mat.shape == (max(lengths), len(THRESHOLDS))
assert [len(twitch) for twitch in model.twitch_list] == list(lengths)
+
+
+def test_gain_saturation_shape_generalises_fuglevand():
+ r = np.array([0.1, 0.41, 0.6, 1, 2, 10])
+ np.testing.assert_array_equal(fuglevand_gain(r, a=2.0, c=3.0), fuglevand_gain(r))
+ a, c = HUMAN_UNIT_FORCE_FREQUENCY
+ s = lambda x: 1 - np.exp(-a * x**c) # noqa: E731
+ expected = np.where(r > 0.4, s(r) / r / (s(0.4) / 0.4), 1.0)
+ np.testing.assert_allclose(fuglevand_gain(r, a=a, c=c), expected, rtol=1e-14)
+
+
+@pytest.mark.parametrize("model_class", MODELS)
+def test_force_frequency_shape_default_is_fuglevand(model_class):
+ assert model_class(THRESHOLDS, 1000 * pq.Hz).force_frequency_shape == (2.0, 3.0)
+
+
+@pytest.mark.parametrize("model_class", MODELS)
+@pytest.mark.parametrize("rate", [5, 10, 20, 40, 120])
+def test_force_matches_closed_form_with_human_unit_shape(model_class, rate):
+ a, c = HUMAN_UNIT_FORCE_FREQUENCY
+ model = model_class(THRESHOLDS, 10000 * pq.Hz, force_frequency_shape=(a, c))
+ block = _spike_block([np.arange(0.01, 2, 1 / rate)] + [[] for _ in range(9)])
+ force = model.generate_force(block, verbose=False)
+ T = model.contraction_times__samples[0] / 10000
+ P = model.peak_twitch_forces__unitless[0]
+ expected = np.e * P * T * rate * fuglevand_gain(T * rate, a=a, c=c)
+ assert force.magnitude[10000:, 0].mean() == pytest.approx(expected, rel=0.005)
+
+
+@pytest.mark.parametrize("model_class", MODELS)
+def test_human_unit_shape_gives_measured_twitch_tetanus_ratio(model_class):
+ # Measured 0.22-0.36 (Macefield et al. 1996; Fuglevand et al. 1999; McNulty et al. 2000); Fuglevand 1993 gives 0.11.
+ model = model_class(THRESHOLDS, 10000 * pq.Hz, force_frequency_shape=HUMAN_UNIT_FORCE_FREQUENCY)
+ force = model.generate_force(_spike_block([np.arange(0.01, 2, 1 / 200)] + [[] for _ in range(9)]),
+ verbose=False)
+ ratio = model.peak_twitch_forces__unitless[0] / force.magnitude[10000:, 0].mean()
+ assert 0.22 <= ratio <= 0.25
+
+
+@pytest.mark.parametrize("model_class", MODELS)
+@pytest.mark.parametrize("shape", [(0, 3), (2, 0), (-1, 1), (1, -2)])
+def test_force_frequency_shape_invalid(model_class, shape):
+ with pytest.raises(ValueError, match="force_frequency_shape"):
+ model_class(THRESHOLDS, 1000 * pq.Hz, force_frequency_shape=shape)
diff --git a/tests/test_literature_register.py b/tests/test_literature_register.py
new file mode 100644
index 00000000..cdc9f52e
--- /dev/null
+++ b/tests/test_literature_register.py
@@ -0,0 +1,21 @@
+"""docs/literature.md is generated by scripts/build_literature.py; code cites sources as [key locator] comments."""
+import ast
+import subprocess
+import sys
+from pathlib import Path
+
+REPO = Path(__file__).resolve().parents[1]
+
+
+def test_register_is_current_and_code_cites_only_registered_keys():
+ r = subprocess.run([sys.executable, str(REPO / "scripts" / "build_literature.py"), "--check"],
+ capture_output=True, text=True, encoding="utf-8")
+ assert r.returncode == 0, r.stderr
+
+
+def test_every_scanned_directory_exists():
+ """A moved folder must not drop out of the citation check silently."""
+ src = (REPO / "scripts" / "build_literature.py").read_text(encoding="utf-8")
+ line = next(l for l in src.splitlines() if l.startswith("CITE_DIRS = "))
+ dirs = ast.literal_eval(line.split("=", 1)[1].strip())
+ assert [d for d in dirs if not (REPO / d).is_dir()] == []
diff --git a/tests/test_network_fan_in.py b/tests/test_network_fan_in.py
new file mode 100644
index 00000000..79bf5b22
--- /dev/null
+++ b/tests/test_network_fan_in.py
@@ -0,0 +1,65 @@
+"""Network.connect_fan_in: every target receives the same number of inputs, each with its target's weight."""
+import numpy as np
+import pytest
+import quantities as pq
+
+from myogen import get_random_generator, load_nmodl_mechanisms, set_random_seed
+from myogen.simulator.neuron import Network
+from myogen.simulator.neuron.populations import AlphaMN__Pool, DescendingDrive__Pool
+
+
+@pytest.fixture(scope="module", autouse=True)
+def mechanisms():
+ load_nmodl_mechanisms()
+
+
+def _network(n_dd=20):
+ set_random_seed(3)
+ mn = AlphaMN__Pool(recruitment_thresholds__array=np.linspace(0.05, 1.0, 6),
+ config_file="alpha_mn_default.yaml")
+ dd = DescendingDrive__Pool(n=n_dd, timestep__ms=0.1 * pq.ms, process_type="gamma", shape=3.0)
+ return Network({"DD": dd, "aMN": mn}), mn, dd
+
+
+def _target(mn, netcon):
+ sec = netcon.syn().get_segment().sec.hname()
+ return next(i for i, cell in enumerate(mn) if sec in {s.hname() for s in [cell.soma, *cell.dend]})
+
+
+def test_each_target_gets_n_excitatory_inputs_with_its_own_weight():
+ net, mn, dd = _network()
+ weights = np.linspace(0.05, 0.02, len(mn))
+ netcons = net.connect_fan_in("DD", "aMN", n_inputs=7, weights__uS=weights * pq.uS, delay__ms=1.0 * pq.ms)
+ targets = np.array([_target(mn, nc) for nc in netcons])
+ np.testing.assert_array_equal(np.bincount(targets, minlength=len(mn)), 7)
+ np.testing.assert_allclose([nc.weight[0] for nc in netcons], weights[targets])
+ assert all(nc.syn().e > -40 for nc in netcons)
+ assert all(nc.delay == 1.0 for nc in netcons)
+ assert net.get_netcons("DD", "aMN") == netcons
+
+
+def test_sources_are_drawn_without_replacement_per_target_in_target_order():
+ net, mn, dd = _network()
+ state = get_random_generator().bit_generator.state
+ netcons = net.connect_fan_in("DD", "aMN", n_inputs=5, weights__uS=0.05 * pq.uS)
+ rng = np.random.default_rng()
+ rng.bit_generator.state = state
+ expected = [list(rng.choice(len(dd), size=5, replace=False)) for _ in mn]
+ source = {cell.ns.hname(): j for j, cell in enumerate(dd)}
+ drawn = [source[nc.pre().hname()] for nc in netcons]
+ assert drawn == [j for row in expected for j in row]
+ assert [_target(mn, nc) for nc in netcons] == [i for i in range(len(mn)) for _ in range(5)]
+
+
+@pytest.mark.parametrize("n_inputs", [0, 21])
+def test_impossible_fan_in_is_refused(n_inputs):
+ net, _, _ = _network()
+ with pytest.raises(ValueError, match="n_inputs"):
+ net.connect_fan_in("DD", "aMN", n_inputs=n_inputs, weights__uS=0.05 * pq.uS)
+
+
+def test_print_network_lists_the_fan_in_connection(capsys):
+ net, _, _ = _network()
+ net.connect_fan_in("DD", "aMN", n_inputs=3, weights__uS=0.05 * pq.uS)
+ net.print_network()
+ assert "DD → aMN [fan-in 3]" in capsys.readouterr().out
diff --git a/tests/test_ramp_recruitment_force.py b/tests/test_ramp_recruitment_force.py
new file mode 100644
index 00000000..f37e798e
--- /dev/null
+++ b/tests/test_ramp_recruitment_force.py
@@ -0,0 +1,60 @@
+"""Recruitment force (% MVC) from an MVC ramp, and motor-unit CV mapped from it."""
+
+import inspect
+
+import numpy as np
+import pytest
+
+from myogen import set_random_seed
+
+from myogen.utils.ramp_analysis import cv_from_recruitment_force, recruitment_force_pct_mvc
+
+FS = 1000.0
+FORCE = np.linspace(0.0, 100.0, int(10 * FS) + 1) # 0 -> 100 % MVC over 10 s
+
+
+def _train(onset_s, rate_hz=10.0, end_s=10.0):
+ return np.arange(onset_s, end_s, 1.0 / rate_hz)
+
+
+def test_recruitment_force_is_the_force_at_the_first_sustained_firing():
+ trains = [_train(1.0), _train(4.25), np.array([]), np.r_[0.5, _train(6.0)]]
+ rt = recruitment_force_pct_mvc(trains, FORCE, FS)
+ np.testing.assert_allclose(rt[:2], [10.0, 42.5], atol=0.1)
+ assert np.isnan(rt[2])
+ assert rt[3] == pytest.approx(60.0, abs=0.1) # an isolated spike at 0.5 s does not count
+
+
+def test_cv_passes_through_the_del_vecchio_band_means():
+ rt = np.array([5.0, 10.0, 15.0, 40.0, 55.0, 65.0, np.nan])
+ cv, (a, b) = cv_from_recruitment_force(rt, residual_sd__m_s=0.0)
+ assert a + b * np.mean([5.0, 10.0, 15.0]) == pytest.approx(3.89)
+ assert a + b * np.mean([55.0, 65.0]) == pytest.approx(5.12)
+ np.testing.assert_allclose(cv[:-1], a + b * rt[:-1])
+ assert cv[-1] == pytest.approx(a + b * 100.0) # never recruited on the ramp: placed at 100 % MVC
+
+
+def test_cv_needs_units_in_both_bands():
+ with pytest.raises(ValueError):
+ cv_from_recruitment_force(np.array([5.0, 10.0, 40.0]))
+
+
+def test_cv_scatter_is_the_del_vecchio_residual():
+ """Units scatter about the line by sqrt(1 - R^2) * SD: R^2 0.70, SD 0.71 m/s [delvecchio2018 p. 2]."""
+ default = inspect.signature(cv_from_recruitment_force).parameters["residual_sd__m_s"].default
+ assert default == pytest.approx(0.39, abs=0.005)
+ rt = np.r_[np.linspace(0.0, 70.0, 4000), 5.0, 60.0]
+ cv, (a, b) = cv_from_recruitment_force(rt, rng=np.random.default_rng(1))
+ residual = cv - (a + b * rt)
+ assert np.std(residual) == pytest.approx(default, abs=0.02)
+ assert np.mean(residual) == pytest.approx(0.0, abs=0.03)
+
+
+def test_cv_scatter_is_seeded():
+ rt = np.array([5.0, 10.0, 15.0, 40.0, 55.0, 65.0])
+ first = cv_from_recruitment_force(rt, rng=np.random.default_rng(3))[0]
+ np.testing.assert_array_equal(first, cv_from_recruitment_force(rt, rng=np.random.default_rng(3))[0])
+ set_random_seed(5)
+ global_first = cv_from_recruitment_force(rt)[0]
+ set_random_seed(5)
+ np.testing.assert_array_equal(global_first, cv_from_recruitment_force(rt)[0])