Ported fixes: iEMG SFAP units and execution route, fibre CV units, Gfluctdv noise

CHANGELOG.md

@@ -9,14 +9,19 @@ 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.
+
 ### Breaking changes
 
 This release changes default outputs; regenerate anything calibrated against earlier versions.
 
 - sEMG outputs change in amplitude, shape, duration and spectrum. Bipolar outputs from 0.5.0–0.13.0 are invalid. Durations, mean/median frequencies and conduction-velocity estimates from those versions are off by about 2× and must be regenerated.
 - Absolute sEMG amplitudes change by one to two orders of magnitude (amplitude no longer depends on conduction velocity). Anything calibrated in µV or mV against earlier outputs must be regenerated. To keep the earlier resistive (measured) skin layer, pass `skin_conductivity__S_per_m=4.88e-4 * pq.S / pq.m` to `Muscle`. To keep the earlier layers, pass `skin_thickness__mm=1.29 * pq.mm, fat_thickness__mm=0.3 * pq.mm`.
 - Pickled `SurfaceEMG` objects from earlier versions lack the new fibre properties; create them again.
-- The `execution=` MUAP disk-cache key changes, so templates cached by earlier versions are recomputed.
+- The `execution=` MUAP disk-cache key changes, so templates cached by earlier versions are recomputed (sEMG and iEMG).
+- Pools that call `insert_Gfluctdv()` now produce noisy discharge: spike trains change, and descending-drive parameters fitted while the noise was silent should be refitted.
+- iEMG outputs change through the fibre conduction velocities (mean ~2.0 → ~3.45 m/s at the default diameter). Muscles with more than ~313 000 fibres get a different, now circular, fibre layout.
+- iEMG amplitudes, shapes and spectra change through the SFAP unit fixes (typically 20–50× larger at a fine wire). Anything tuned against earlier iEMG amplitudes, including `add_noise` floors given in absolute units, must be redone. Values are in mV.
 - iEMG outputs change through the narrower end-plate spreads. Pass `IntramuscularEMG(..., nmj_coefficients="reference_code")` to reproduce earlier end-plates.
 
 ### Fixed
@@ -34,6 +39,15 @@ This release changes default outputs; regenerate anything calibrated against ear
 - **Default skin and fat thickness (FDI).** Skin 1.29 → 1.2 mm (hand dorsum median, Li et al. 2018) and fat 0.3 → 1.3 mm, so the total depth to muscle is 2.5 mm, the low end of the FDI depth measured by Nishikawa et al. (2024); the fat value is derived from that total, not measured. Surface amplitudes from the default `Muscle` fall by about 40 %.
 - **Failed motor units raise.** On the default `SurfaceEMG` route an exception while simulating one motor unit was logged and replaced by an all-zero MUAP, which passes for a valid, silent unit. The exception is now re-raised with the unit and array index; the `ExecutionPolicy` route already raised.
 - **Angular sampling overflow raises.** With `sampling_points_in_theta_domain` ≥ 256 the scaled modified Bessel functions leave the float64 range at the highest angular orders (I_n underflows to 0, K_n overflows). The volume conductor replaced the overflowing entries with zero, which gave a singular system on the CPU and NaN MUAPs on the GPU. It now raises a `ValueError` that names the setting.
+- **Gfluctdv membrane noise was silent.** `AlphaMN.insert_Gfluctdv()` inserted the Destexhe (2001) fluctuating-conductance mechanism but left its GLOBAL multipliers `multex_Gfluctdv` / `multin_Gfluctdv` at the `.mod` default of 0, so it added no noise. They are now set to 1 on insertion. Any pool that calls `insert_Gfluctdv()` now has conductance noise and more variable discharge; this also makes the `gfluctdv` search in `examples/02_finetune/01_optimize_dd_for_target_firing_rate.py` effective.
+- **Fibre conduction velocities.** `Muscle` mixed units in the velocity–diameter line (5.5 (m/s)/mm with diameters in mm), giving ~1.1 m/s, which the 2 m/s clip then raised: every fibre had exactly 2.0 m/s. Fibre velocities now follow 0.95 + 0.05·d[µm] m/s (Nandedkar & Stålberg 1983; ~3.45 m/s at 50 µm), the same line `SurfaceFiberDispersion` uses; noise (SD 0.2 m/s) and clip (2–6 m/s) are unchanged. Intramuscular EMG uses these velocities; surface EMG uses the mean velocity and is unaffected.
+- **Intramuscular SFAP units and polarity.** Four errors in the iEMG source and volume conductor (`bioelectric.py`, `MotorUnitSim.calc_sfaps`) are corrected on both the default and `execution=` routes (`runtime/density.py`, `runtime/iemg.py`); a single-fibre SFAP now matches an independent SI line-source solution within 3 % at 0.1–1 mm (`tests/test_iemg_sfap_units.py`).
+  - Muscle conductivities defaulted to 63 / 330 "S/m" (1000× the tissue values of 0.063 / 0.33 S/m), making potentials 1000× too small. Defaults are now 0.063 / 0.33 S/m.
+  - The Rosenfalck IAP was evaluated at twice the distance and differentiated with the chain rule applied twice: an IAP half the published length and 8× the current. It now uses λ = 1 mm, as the sEMG model does.
+  - The wave travelling towards the shorter fibre half was subtracted instead of added, so it had inverted polarity; the two halves cancelled near the end-plate and in the far field.
+  - The sum over the fibre omitted the spatial step, so amplitudes scaled with `1/spatial_resolution__mm`.
+  Together the bugs nearly cancelled at the default 0.01 mm step for a nearby electrode. On an FDI cylinder (bipolar fine wire, 2 mm IED, 15 units) MUAPs go from 0.03–1.8 mV to 0.9–37 mV peak-to-peak. These are point electrodes in an infinite homogeneous medium, which overestimates the near field (no electrode area, no near-field correction); see Breaking changes.
+- **Large circular fibre layouts.** Fibre centres came from a 399 499-point table covering a square, of which only ~313 000 points lie inside the circle. Beyond that the outer ring was filled only at the diagonals, and beyond 399 499 the fibre count was silently capped. Layouts the table covers are unchanged (FDI defaults ≈ 54 000 fibres); larger counts now use a scrambled Halton sequence in the disc, seeded from MyoGen's generator, with a warning.
 
 ### Added
 

myogen/simulator/core/emg/fiber_simulation.py

@@ -55,8 +55,8 @@ def compute_intramuscular_kernel(
     z_grid: np.ndarray,
     electrode_z: np.ndarray,
     r: float,
-    sigma_r: float = 63.0,
-    sigma_z: float = 330.0,
+    sigma_r: float = 0.063,
+    sigma_z: float = 0.33,
 ) -> np.ndarray:
     """
     Compute the intramuscular volume conductor kernel using the Green's function.
@@ -72,9 +72,9 @@ def compute_intramuscular_kernel(
         Electrode z-positions in mm, shape (n_electrodes,).
     r : float
         Radial distance from fiber to electrode in mm.
-    sigma_r : float, default=63.0
+    sigma_r : float, default=0.063
         Radial conductivity in S/m.
-    sigma_z : float, default=330.0
+    sigma_z : float, default=0.33
         Longitudinal conductivity in S/m.
 
     Returns

myogen/simulator/core/emg/intramuscular/bioelectric.py

@@ -4,14 +4,14 @@
 This module contains the core bioelectric modeling functions for simulating
 single fiber action potentials (SFAPs) and motor unit action potentials (MUAPs)
 in intramuscular EMG. The functions implement the volume conductor models
-from Farina et al. 2004 and the transmembrane current models from
-Rosenfalck 1969.
+from Farina et al. 2004 [farina2004] and the transmembrane current models from
+Rosenfalck 1969 [rosenfalck1969].
 
 References
 ----------
-[1] Farina, D., Merletti, R., 2001. A novel approach for precise simulation of the EMG signal detected by surface electrodes. IEEE Transactions on Biomedical Engineering 48, 637–646. <br>
-[2] Rosenfalck, P., 1969. Intra- and extracellular potential fields of active nerve and muscle fibres. Acta Physiologica Scandinavica Supplementum 321, 1–168. <br>
-[3] Nandedkar, S.D., Stålberg, E., 1983. Simulation of single muscle fibre action potentials. Medical & Biological Engineering & Computing 21, 158–165.
+[1] Farina, D., Merletti, R., 2001. A novel approach for precise simulation of the EMG signal detected by surface electrodes. IEEE Transactions on Biomedical Engineering 48, 637–646. [farinamerletti2001] <br>
+[2] Rosenfalck, P., 1969. Intra- and extracellular potential fields of active nerve and muscle fibres. Acta Physiologica Scandinavica Supplementum 321, 1–168. [rosenfalck1969] <br>
+[3] Nandedkar, S.D., Stålberg, E., 1983. Simulation of single muscle fibre action potentials. Medical & Biological Engineering & Computing 21, 158–165. [nandedkar1983]
 """
 
 import numpy as np
@@ -22,7 +22,7 @@ def get_tm_current(z: np.ndarray, D1: float = 96.0, D2: float = -90.0) -> np.nda
     Calculate transmembrane current using Rosenfalck's model.
 
     This function implements the transmembrane current model from:
-    P. Rosenfalck "Intra and extracellular fields of active nerve and muscle fibers" (1969)
+    P. Rosenfalck "Intra and extracellular fields of active nerve and muscle fibers" (1969) [rosenfalck1969]
 
     Parameters
     ----------
@@ -45,7 +45,7 @@ def get_tm_current(z: np.ndarray, D1: float = 96.0, D2: float = -90.0) -> np.nda
 
 def get_tm_current_dz(z: np.ndarray, D1: float = 96.0, xp=np) -> np.ndarray:
     """
-    Calculate first derivative of transmembrane current (Rosenfalck model).
+    Calculate first derivative of transmembrane current (Rosenfalck model) [rosenfalck1969].
 
     This is the spatial derivative of the transmembrane current model used
     for action potential propagation simulation.
@@ -73,7 +73,7 @@ def get_tm_current_dz(z: np.ndarray, D1: float = 96.0, xp=np) -> np.ndarray:
 
 def get_tm_current_ddz(z: np.ndarray, D1: float = 96.0) -> np.ndarray:
     """
-    Calculate second derivative of transmembrane current (Rosenfalck model).
+    Calculate second derivative of transmembrane current (Rosenfalck model) [rosenfalck1969].
 
     Parameters
     ----------
@@ -98,16 +98,16 @@ def get_elementary_current_response(
     z: np.ndarray,
     z_electrode: float,
     r: np.ndarray,
-    sigma_r: float = 63.0,  # S/m
-    sigma_z: float = 330.0,  # S/m
+    sigma_r: float = 0.063,  # S/m
+    sigma_z: float = 0.33,  # S/m
     xp=np,
 ) -> np.ndarray:
     """
     Calculate elementary current response for volume conductor.
 
     This function calculates the potential response at electrode location
     due to a unit current source at different positions along the muscle fiber.
-    Based on Nandedkar & Stålberg 1983.
+    Based on Nandedkar & Stålberg 1983 [nandedkar1983].
 
     Parameters
     ----------
@@ -117,22 +117,23 @@ def get_elementary_current_response(
         Electrode position along z-axis in mm
     r : np.ndarray
         Radial distance from fiber to electrode in mm
-    sigma_r : float, default=63.0
-        Radial conductivity in S/m (from Andreassen & Rosenfalck 1980)
-    sigma_z : float, default=330.0
-        Longitudinal conductivity in S/m (from Andreassen & Rosenfalck 1980)
+    sigma_r : float, default=0.063
+        Radial conductivity in S/m (from Andreassen & Rosenfalck 1980) [andreassenrosenfalck1980]
+    sigma_z : float, default=0.33
+        Longitudinal conductivity in S/m (from Andreassen & Rosenfalck 1980) [andreassenrosenfalck1980]
     xp : module, default=np
         Array backend (numpy or cupy)
 
     Returns
     -------
     np.ndarray
-        Elementary current response (transfer function)
+        Transfer function in Ω: potential (mV) per mA of point current. Multiply the current per unit length from ``get_current_density``
+        (mA/mm) by this and by the spatial step dz (mm), and sum over z, to get the SFAP in mV.
     """
-    # ---- FIXED UNIT CONVERSIONS ----
-    # Convert conductivities from S/m to S/mm to match spatial units (mm)
-    sigma_r_S_per_mm = sigma_r / 1000.0  # CORRECTED: convert S/m → S/mm
-    sigma_z_S_per_mm = sigma_z / 1000.0  # CORRECTED: convert S/m → S/mm
+    # Unit system of the iEMG kernel: mm, mV, mA. 1 S/m = 1e-3 S/mm. The defaults were
+    # 63 / 330 "S/m" until 0.13.0 (1000x the tissue values), which made SFAPs 1000x too small.
+    sigma_r_S_per_mm = sigma_r / 1000.0
+    sigma_z_S_per_mm = sigma_z / 1000.0
 
     # Normalize inputs to computation backend (prevents numpy/cupy mixing)
     z = xp.asarray(z)
@@ -228,7 +229,7 @@ def get_current_density(
 ):
     """
     Model the individual action potential (IAP) or single fiber action potential (SFAP) in space and time.
-    Translated from Farina & Merletti (2001) and Nandedkar & Stålberg (1998).
+    Translated from Farina & Merletti (2001) [farinamerletti2001] and Nandedkar & Stålberg (1998).
 
     Parameters
     ----------
@@ -276,14 +277,17 @@ def tendon_terminator(z_inline, L_inline):
 
     # Compute psi (transmembrane current derivative)
     if L1 >= L2:
-        psi = -4 * get_tm_current_dz(-2 * (Z - zi - v * T), xp=xp)
+        # psi = dV_m/dz for the wave travelling towards +z, zeta = v t - (z - zi) in mm (Rosenfalck
+        # length scale 1 mm, as in the sEMG model). Until 0.13.0 this read -4 * V'(-2 * arg): the chain
+        # rule counted twice and the IAP compressed to half its length (8x the current).
+        psi = -get_tm_current_dz(-(Z - zi - v * T), xp=xp)
         longest_wave = xp.diff(psi, axis=0) / dz
         longest_wave *= tendon_terminator(Z[:-1, :] - zi - L1 / 2, L1)
         # Explicit bool→float64 cast required: CuPy does not support
         # implicit multiplication of bool arrays with float arrays.
         longest_wave *= ((Z[:-1, :] - zi) / v > 0).astype(xp.float64)
     else:
-        psi = 4 * get_tm_current_dz(-2 * (-Z + zi - v * T), xp=xp)
+        psi = get_tm_current_dz(-(-Z + zi - v * T), xp=xp)
         longest_wave = xp.diff(psi, axis=0) / dz
         longest_wave *= tendon_terminator(Z[:-1, :] - zi + L2 / 2, L2)
         longest_wave *= ((-Z[:-1, :] + zi) / v > 0).astype(xp.float64)  # bool→float64
@@ -296,10 +300,13 @@ def tendon_terminator(z_inline, L_inline):
 
     if L1 >= L2:
         shortest_wave *= tendon_terminator(Z[:-1, :] - zi + L2 / 2, L2)
-        iap = longest_wave - shortest_wave
+        # The mirrored array is already d2V/dz2 of the opposite wave (a second derivative is
+        # unchanged by reflection), so the two waves add. Until 0.13.0 they were subtracted,
+        # which inverted the wave travelling away from the longer half.
+        iap = longest_wave + shortest_wave
     else:
         shortest_wave *= tendon_terminator(Z[:-1, :] - zi - L1 / 2, L1)
-        iap = shortest_wave - longest_wave
+        iap = shortest_wave + longest_wave
 
     # Suppress endplate density if required
     if suppress_endplate_density:
@@ -315,7 +322,8 @@ def tendon_terminator(z_inline, L_inline):
     # Compute cross-sectional area in mm²
     area_mm2 = xp.pi * (d / 2) ** 2  # CORRECTED: removed extra /4
 
-    # Scale current density by intracellular conductivity and fiber cross-section area
+    # Scale current density by intracellular conductivity and fiber cross-section area:
+    # transmembrane current per unit length, S/mm * mm^2 * mV/mm^2 = mA/mm
     iap *= sigma_i * area_mm2
 
     return iap
@@ -329,7 +337,7 @@ def get_current_density_fast(
     L1: float,
     L2: float,
     v: float,
-    d: float = 55e-6,
+    d: float = 55e-3,
     suppress_endplate_density: bool = True,
     endplate_width: float = 0.5,
 ) -> np.ndarray:
@@ -355,7 +363,7 @@ def get_current_density_fast(
         Length from endplate to negative tendon in mm
     v : float
         Conduction velocity in mm/s
-    d : float, default=55e-6
+    d : float, default=55e-3
         Fiber diameter in mm
     suppress_endplate_density : bool, default=True
         Whether to suppress endplate region

myogen/simulator/core/emg/intramuscular/intramuscular_emg.py

@@ -314,7 +314,8 @@ def _calculate_muaps(self, n_jobs: int = -2, verbose: bool = True) -> INTRAMUSCU
         Returns
         -------
         INTRAMUSCULAR_MUAP__Block
-            Intramuscular MUAP shapes for all motor units stored in a neo.Block.
+            Intramuscular MUAP shapes for all motor units stored in a neo.Block, in mV
+            (tagged ``pq.dimensionless`` for backward compatibility).
         """
         if not self._motor_units:
             raise ValueError("Must call _initialize_motor_units() first")
@@ -516,8 +517,8 @@ def simulate_intramuscular_emg(
         Returns
         -------
         INTRAMUSCULAR_EMG__Block
-            Intramuscular EMG signals for the electrode array stored in a neo.Block.
-            Results are stored in the `intramuscular_emg__Block` property after execution.
+            Intramuscular EMG signals for the electrode array stored in a neo.Block, in mV
+            (tagged ``pq.dimensionless`` for backward compatibility). Results are stored in the `intramuscular_emg__Block` property after execution.
 
         Raises
         ------

myogen/simulator/core/emg/intramuscular/motor_unit_sim.py

@@ -321,7 +321,8 @@ def calc_sfaps(
                         ))
 
                     ecr_all = cp.column_stack(ecr_list)
-                    sfaps_dev[:, :, fiber_idx] = current_density.T @ ecr_all
+                    # mA/mm x Ω x dz (mm) summed over z: SFAP in mV
+                    sfaps_dev[:, :, fiber_idx] = (current_density.T @ ecr_all) * dz
 
             cp.cuda.Device(0).synchronize()
             self.sfaps = cp.asnumpy(sfaps_dev)
@@ -378,7 +379,8 @@ def calc_sfaps(
                     ))
 
                 ecr_all = np.column_stack(ecr_list)
-                self.sfaps[:, :, fiber_idx] = current_density.T @ ecr_all
+                # mA/mm x Ω x dz (mm) summed over z: SFAP in mV
+                self.sfaps[:, :, fiber_idx] = (current_density.T @ ecr_all) * dz
 
         self.shift_sfaps(dt)
 

myogen/simulator/core/emg/runtime/density.py

@@ -22,13 +22,17 @@ def _density(t, z, active, mirror, short_mask, end_mask, zi, v, dz, positive, sc
         for col in range(columns):
             value = 0.0
             if needed:
+                # zeta = v t - |z - zi| in mm, Rosenfalck length scale 1 mm, as
+                # in bioelectric.get_current_density (until 0.13.0 both read
+                # -2 * (...) and scaled by 4: half the IAP length, 8x the current).
                 if positive:
-                    x = -2 * (z[row] - zi - v * t[col])
+                    x = -(z[row] - zi - v * t[col])
                 else:
-                    x = -2 * (-z[row] + zi - v * t[col])
+                    x = -(-z[row] + zi - v * t[col])
                 if x > 0:
                     value = 96.0 * (3 * x**2 - x**3) * math.exp(-x)
-                value *= -4 if positive else 4
+                if positive:
+                    value = -value
             wave[row, col] = value
     # Ascending rows retain the unmodified next endpoint until it is consumed.
     for row in range(rows):
@@ -42,7 +46,9 @@ def _density(t, z, active, mirror, short_mask, end_mask, zi, v, dz, positive, sc
         source = mirror[row]
         for col in range(columns):
             short = wave[source, col] if source >= 0 and short_mask[row] else 0.0
-            value = wave[row, col] - short if positive else short - wave[row, col]
+            # The mirrored wave is already d2V/dz2 of the opposite wave, so the
+            # two add (they were subtracted until 0.13.0).
+            value = wave[row, col] + short
             result[row, col] = (value if end_mask[row] else 0.0) * scale
     return result