From 21a412f27eff19bed8db37bb3e98a5f8beee4c28 Mon Sep 17 00:00:00 2001 From: Devon Rohlf Date: Fri, 30 Oct 2026 17:39:20 +0100 Subject: [PATCH 01/42] fix(semg): complex B(k_z), physical IAP scale, fibre timing and geometry - B(k_z) keeps its imaginary part in the Numba kernel, the NumPy fallback and the phi kernel; longitudinal bipolar outputs were previously lost. - Remove the ad-hoc `z /= 2` IAP stretch: the Rosenfalck source now uses the same mm grid as k_z (matches the time-domain hybrid path). - MUAP sampling rate follows the internal grid (muap_sampling_rate__Hz), shared by the default and execution routes. - Drop the per-fibre peak re-centring, which erased propagation delays and fibre-to-fibre timing within a motor unit. - L1/L2 are the +z/-z semi-lengths from the end-plate (fiber_semi_lengths). - Each electrode array keeps its own output shape. - Bump the execution-route cache key so stale templates are recomputed. --- CHANGELOG.md | 16 ++ myogen/simulator/core/emg/runtime/semg.py | 16 +- .../core/emg/surface/simulate_fiber.py | 111 ++++++---- .../simulator/core/emg/surface/surface_emg.py | 88 ++++---- tests/test_fiber_simulation.py | 6 +- tests/test_surface_model_consistency.py | 189 ++++++++++++++++++ 6 files changed, 325 insertions(+), 101 deletions(-) create mode 100644 tests/test_surface_model_consistency.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 96355418..931a6134 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,22 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +This release corrects the surface-EMG pipeline against Farina, Mesin, Martina & Merletti (2004). The problems were reported in an external tutorial notebook. Both the default route and `execution=` are affected. + +### Fixed + +- **Longitudinal bipolar channels.** The angular integral B(k_z) is now complex. Keeping only its real part made `bipolar_longitudinal` channels about 10⁻⁷ of their true value and distorted rotated arrays. A bipolar channel now equals the difference of the two monopolar electrodes. +- **MUAP time axis.** MUAP templates are labelled with the rate of the grid they were computed on (window `K / v`). The previous label was 2× too high, so MUAPs looked half as long, spectra sat an octave high, and conduction velocity estimated from the output was 2× the value set. +- **Rosenfalck IAP length scale.** The IAP is evaluated on the simulation's own z axis (λ = 1 mm, as published), not on `z / 2`. The stretched IAP was 2× too long; it now agrees with the independent time-domain path (`simulate_fiber_hybrid`). +- **Fibre semi-lengths.** `L1` (towards +z) and `L2` were swapped, which moved both tendons by 2·z_i with the end-plate. Fibres now span `±fiber_length / 2` with the end-plate at z_i (`fiber_semi_lengths`). +- **Fibre timing within a motor unit.** Each fibre's potential was circularly shifted onto its own peak before summation, misaligning fibres (up to ~60 % MUAP error on bipolar channels). All fibres now share the firing instant at the window centre. +- **Mixed electrode-array shapes.** On the default route each array's MUAPs have that array's own grid shape; previously every array was broadcast into the first array's shape. + +### Upgrade notes + +- 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. +- The `execution=` MUAP disk-cache key changes, so templates cached by earlier versions are recomputed. + ## [0.13.0] - 2026-10-28 Thanks to @joaohbbittar for the detailed report in #26, this release corrects the Fuglevand force model and the Hill muscle activation. diff --git a/myogen/simulator/core/emg/runtime/semg.py b/myogen/simulator/core/emg/runtime/semg.py index f88332be..980deff2 100644 --- a/myogen/simulator/core/emg/runtime/semg.py +++ b/myogen/simulator/core/emg/runtime/semg.py @@ -8,7 +8,11 @@ from myogen import get_random_generator from myogen.utils.neo import create_grid_signal from ..electrodes import SurfaceElectrodeArray -from ..surface.simulate_fiber import _simulate_fiber_v2_python +from ..surface.simulate_fiber import ( + _simulate_fiber_v2_python, + fiber_semi_lengths, + muap_sampling_rate__Hz, +) from .session import ExecutionSession from .surface_preparation import SurfacePreparation from .cache import DiskCache, fingerprint @@ -52,8 +56,8 @@ def simulate_muaps(owner, policy, n_jobs, verbose, use_gpu): rows, cols = next(iter(shapes)) shape = (rows, cols, N) final_samples = owner._sampling_points_in_t_and_z_domains - duration = (kernel_length / 2.0) / owner._mean_conduction_velocity__m_s / 1000 - rate = (final_samples - 1) / duration + rate = muap_sampling_rate__Hz(owner._mean_conduction_velocity__m_s, N, final_samples, + kernel_length) # Source-grid arithmetic retains the original precision and complete # N x N grids. Estimate all live complex/native temporaries generously. workspace = N * N * 16 * 24 + N * M * 8 * 180 + rows * cols * N * 16 * 4 @@ -67,7 +71,7 @@ def simulate_muaps(owner, policy, n_jobs, verbose, use_gpu): radius = float(electrode.electrode_radius__mm.rescale(pq.mm).magnitude) preparation = SurfacePreparation() dependencies = ( - "surface-muap-gitlab-fc077-v1", N, M, final_samples, kernel_length, rate, + "surface-muap-gitlab-fc077-v2", N, M, final_samples, kernel_length, rate, owner._internal_sampling_frequency__Hz, owner._mean_conduction_velocity__m_s, owner._radius_total, owner._radius_bone__mm, owner._fat_thickness__mm, owner._skin_thickness__mm, @@ -94,6 +98,7 @@ def simulate_muaps(owner, policy, n_jobs, verbose, use_gpu): def compute(start, stop): values = np.empty((stop - start, *shape)) for local, fibre in enumerate(range(start, stop)): + L1, L2 = fiber_semi_lengths(zone, lengths[fibre]) values[local], _, _ = _simulate_fiber_v2_python( Fs=float(owner._internal_sampling_frequency__Hz * 1e-3), v=float(owner._mean_conduction_velocity__m_s), @@ -101,8 +106,7 @@ def compute(start, stop): r_bone=float(owner._radius_bone__mm), th_fat=float(owner._fat_thickness__mm), th_skin=float(owner._skin_thickness__mm), - R=radial[fibre], L1=abs(zone + lengths[fibre] / 2), - L2=abs(zone - lengths[fibre] / 2), zi=zone, + R=radial[fibre], L1=L1, L2=L2, zi=zone, electrode_array=electrode, sig_muscle_rho=float(owner._muscle_conductivity_radial__S_m), sig_muscle_z=float(owner._muscle_conductivity_longitudinal__S_m), diff --git a/myogen/simulator/core/emg/surface/simulate_fiber.py b/myogen/simulator/core/emg/surface/simulate_fiber.py index 89a6f6b7..bc4efba6 100644 --- a/myogen/simulator/core/emg/surface/simulate_fiber.py +++ b/myogen/simulator/core/emg/surface/simulate_fiber.py @@ -101,6 +101,37 @@ def f_minus_t(y): return y[::-1].copy() +def fiber_semi_lengths(innervation_zone__mm: float, fiber_length__mm: float) -> tuple[float, float]: + """Return ``(L1, L2)`` for a fibre centred on z = 0 with its end-plate at ``zi``. + + In the source term (Farina et al. 2004, eq. 1) ``L1`` is the semi-length + towards +z and ``L2`` towards -z, so the fibre spans ``[zi - L2, zi + L1]``. + Keeping the tendons at ``±fiber_length / 2`` gives ``L1 = l/2 - zi`` and + ``L2 = l/2 + zi``. + """ + L1 = fiber_length__mm / 2 - innervation_zone__mm + L2 = fiber_length__mm / 2 + innervation_zone__mm + if L1 <= 0 or L2 <= 0: + raise ValueError( + f"End-plate at {innervation_zone__mm} mm lies outside a fibre of length {fiber_length__mm} mm" + ) + return L1, L2 + + +def muap_sampling_rate__Hz( + conduction_velocity__m_s: float, internal_points: int, output_points: int, kernel_length__mm: float +) -> float: + """Sampling rate of a fibre window resampled from ``internal_points`` to ``output_points``. + + The internal window has ``internal_points`` samples at + ``Fs_eff = v * (N - 1) / kernel_length`` (see ``_simulate_fiber_v2_python``), + i.e. a period of ``N / Fs_eff``. Periodic (Fourier) resampling keeps that + period, so the output rate is ``output_points * Fs_eff / N``. + """ + fs_effective__kHz = conduction_velocity__m_s * (internal_points - 1) / kernel_length__mm + return output_points * fs_effective__kHz / internal_points * 1e3 + + # Numba-optimized helper functions for simulate_fiber_v3 def _get_numba_functions(): """Get Numba-optimized functions if available, otherwise return None.""" @@ -119,24 +150,32 @@ def _compute_B_kz_fast( len_kz, len_ktheta, ): - """Numba-optimized computation of B_kz matrix.""" - B_kz = np.zeros((channels_0, channels_1, len_kz)) + """Numba-optimized computation of the complex B_kz matrix (Farina 2004, eq. 5). + + B_kz is complex in general: a longitudinal bipolar filter is purely + imaginary, so its whole signal lives in the imaginary part. + """ + B_kz = np.zeros((channels_0, channels_1, len_kz), dtype=np.complex128) for channel_z in prange(channels_0): for channel_theta in prange(channels_1): pos_theta_val = pos_theta[channel_z, channel_theta] for i in range(len_kz): sum_real = 0.0 + sum_imag = 0.0 for j in range(len_ktheta): phase = pos_theta_val * ktheta_mesh_kzktheta[i, j] cos_phase = math.cos(phase) sin_phase = math.sin(phase) # Complex multiplication: H_glo * exp(1j * phase) - result_real = H_glo_real[i, j] * cos_phase - H_glo_imag[i, j] * sin_phase - sum_real += result_real + sum_real += H_glo_real[i, j] * cos_phase - H_glo_imag[i, j] * sin_phase + sum_imag += H_glo_real[i, j] * sin_phase + H_glo_imag[i, j] * cos_phase - B_kz[channel_z, channel_theta, i] = sum_real * k_theta_diff / (2 * math.pi) + B_kz[channel_z, channel_theta, i] = complex( + sum_real * k_theta_diff / (2 * math.pi), + sum_imag * k_theta_diff / (2 * math.pi), + ) return B_kz @@ -144,7 +183,8 @@ def _compute_B_kz_fast( def _compute_phi_fast( I_kzkt_real, I_kzkt_imag, - B_kz, + B_kz_real, + B_kz_imag, pos_z, kz_mesh_kzkt, channels_0, @@ -166,9 +206,10 @@ def _compute_phi_fast( for i in range(len_kz): # I_kzkt[i, j] * B_kz[channel_z, channel_theta, i] - B_val = B_kz[channel_z, channel_theta, i] - arg_real = I_kzkt_real[i, j] * B_val - arg_imag = I_kzkt_imag[i, j] * B_val + B_re = B_kz_real[channel_z, channel_theta, i] + B_im = B_kz_imag[channel_z, channel_theta, i] + arg_real = I_kzkt_real[i, j] * B_re - I_kzkt_imag[i, j] * B_im + arg_imag = I_kzkt_real[i, j] * B_im + I_kzkt_imag[i, j] * B_re # arg * exp(1j * pos_z * kz_mesh_kzkt[i, j]) phase = pos_z_val * kz_mesh_kzkt[i, j] @@ -284,10 +325,11 @@ def _simulate_fiber_v2_python( **Recommended**: Set this explicitly (e.g., 50-100 mm for typical motor units) to ensure MUAP duration is physically accurate and independent of N. - The actual MUAP duration will be approximately: fiber_length__mm / (2 * v) ms - (factor of 2 due to internal scaling for Rosenfalck model). + The returned window spans fiber_length__mm / v ms (``N`` samples at + ``v * (N - 1) / fiber_length__mm`` kHz), with the fibre's firing + instant at sample ``N // 2``. See :func:`muap_sampling_rate__Hz`. - Example: fiber_length__mm=100 with v=4 m/s → ~12.5 ms MUAP duration. + Example: fiber_length__mm=100 with v=4 m/s → 25 ms window. Returns ------- @@ -327,11 +369,8 @@ def _simulate_fiber_v2_python( # Calculate effective Fs based on grid spacing # When fiber_length__mm is specified, we must use the corresponding Fs for FFT consistency if fiber_length__mm is not None: - # z-grid spans fiber_length__mm with N points (before z /= 2 scaling) - # Grid spacing (before scaling): dz = fiber_length__mm / (N-1) + # z-grid spans fiber_length__mm with N points: dz = fiber_length__mm / (N-1) # For FFT: dz = v / Fs_eff → Fs_eff = v * (N-1) / fiber_length__mm - # Note: The z /= 2 scaling affects the IAP kernel evaluation but the FFT - # operates on the array structure, so we use the pre-scaling grid spacing Fs_effective = v * (N - 1) / fiber_length__mm # kHz else: Fs_effective = Fs @@ -376,10 +415,9 @@ def _simulate_fiber_v2_python( # N controls only the resolution, not the physical extent z = np.linspace(-fiber_length__mm / 2, fiber_length__mm / 2, N, dtype=np.longdouble) - z /= 2 # mm - scaling for Rosenfalck model - - # Rosenfalck action potential: A * exp(-z) * (3z² - z³) for z ≥ 0 - # Uses raw z values in mm (NOT normalized) to match original implementation + # Rosenfalck action potential: A * exp(-z) * (3z² - z³) for z ≥ 0, with z in mm + # on the same grid as k_z (length scale 1 mm, as published). The source grid + # and the volume-conductor grid must share one z axis. aux = np.zeros_like(z) for i in range(len(z)): if z[i] >= 0: @@ -486,7 +524,7 @@ def _simulate_fiber_v2_python( channels[0], channels[1], len(k_z), len(k_theta), ) else: - B_kz = np.zeros((channels[0], channels[1], len(k_z))) + B_kz = np.zeros((channels[0], channels[1], len(k_z)), dtype=np.complex128) for channel_z in range(channels[0]): for channel_theta in range(channels[1]): arg = np.multiply( @@ -505,7 +543,8 @@ def _simulate_fiber_v2_python( if _prepared is not None else _numba_phi_func( np.ascontiguousarray(I_kzkt.real, dtype=np.float64), np.ascontiguousarray(I_kzkt.imag, dtype=np.float64), - np.ascontiguousarray(B_kz, dtype=np.float64), + np.ascontiguousarray(B_kz.real, dtype=np.float64), + np.ascontiguousarray(B_kz.imag, dtype=np.float64), np.ascontiguousarray(pos_z_mm, dtype=np.float64), np.ascontiguousarray(kz_mesh_kzkt, dtype=np.float64), channels[0], channels[1], len(k_z), len(k_t), @@ -549,24 +588,9 @@ def _simulate_fiber_v2_python( ) ) - # Center the MUAP signal in the time window by finding the peak and shifting - # For each electrode channel, find the peak and center it - N_time = phi.shape[2] - center_idx = N_time // 2 - - # Find the peak location (use center electrode as reference) - center_electrode_row = phi.shape[0] // 2 - center_electrode_col = phi.shape[1] // 2 - center_signal = phi[center_electrode_row, center_electrode_col, :] - - # Find peak (max absolute value) - peak_idx = np.argmax(np.abs(center_signal)) - - # Calculate shift needed to center the peak - shift = center_idx - peak_idx - - # Apply circular shift to all channels - phi = np.roll(phi, shift, axis=2) + # No per-fibre re-centring: the fibre fires at sample N // 2 (the IAP grid is + # centred on z = 0), and every fibre of a motor unit must keep that common + # instant so that the summed MUAP preserves the fibres' relative timing. # print(f"DEBUG: Final phi range = [{np.min(phi):.6e}, {np.max(phi):.6e}]") # print(f"DEBUG: Final phi peak-to-peak = {np.ptp(phi):.6e}") @@ -668,10 +692,11 @@ def simulate_fiber_v2( **Recommended**: Set this explicitly (e.g., 50-100 mm for typical motor units) to ensure MUAP duration is physically accurate and independent of N. - The actual MUAP duration will be approximately: fiber_length__mm / (2 * v) ms - (factor of 2 due to internal scaling for Rosenfalck model). + The returned window spans fiber_length__mm / v ms (``N`` samples at + ``v * (N - 1) / fiber_length__mm`` kHz), with the fibre's firing + instant at sample ``N // 2``. See :func:`muap_sampling_rate__Hz`. - Example: fiber_length__mm=100 with v=4 m/s → ~12.5 ms MUAP duration. + Example: fiber_length__mm=100 with v=4 m/s → 25 ms window. Returns ------- diff --git a/myogen/simulator/core/emg/surface/surface_emg.py b/myogen/simulator/core/emg/surface/surface_emg.py index 11153adc..b88aca86 100644 --- a/myogen/simulator/core/emg/surface/surface_emg.py +++ b/myogen/simulator/core/emg/surface/surface_emg.py @@ -22,7 +22,12 @@ from myogen import get_random_generator from myogen.simulator.core.emg.electrodes import SurfaceElectrodeArray -from myogen.simulator.core.emg.surface.simulate_fiber import simulate_fiber_v2, _simulate_fiber_v2_python +from myogen.simulator.core.emg.surface.simulate_fiber import ( + _simulate_fiber_v2_python, + fiber_semi_lengths, + muap_sampling_rate__Hz, + simulate_fiber_v2, +) from myogen.simulator.core.muscle import Muscle from myogen.utils.decorators import beartowertype from myogen.utils.types import ( @@ -72,18 +77,13 @@ class SurfaceEMG: If None, defaults to 10 kHz. Higher values provide better temporal resolution but increase computation time. Default is 10 kHz. iap_kernel_length__mm : float, optional - Physical spatial extent for intracellular action potential (IAP) kernel evaluation in mm. - If None (default), uses individual fiber lengths from muscle model, ensuring - MUAP duration is physiologically accurate and independent of sampling resolution. - - **Recommended**: Leave as None to use fiber-specific lengths for most realistic MUAPs. - - Alternatively, set to a fixed value (e.g., 80-100 mm) to use the same kernel - extent for all fibers, which can simplify analysis but may be less physiologically - accurate for muscles with variable fiber lengths. + Length ``K`` of the periodic z window (mm) on which each fibre is simulated. + If None (default), ``K = 2.5 × mean_fiber_length__mm``. - This parameter controls the spatial extent over which the IAP waveform is computed, - directly affecting MUAP duration: duration ≈ iap_kernel_length__mm / (2 * v) ms. + The MUAP template window spans ``K / v`` ms, with the firing instant at its + centre. ``K`` must hold the whole fibre plus the IAP, or the signal wraps + around; it does not change the physical MUAP, only the window and the + internal sampling step ``K / (v (N - 1))``. Attributes ---------- @@ -306,20 +306,16 @@ def simulate_muaps( for mu_idx in range(n_motor_units) ] - # Pre-allocate result shape at INTERNAL resolution (optimization: avoid repeated shape calculations) - # Will be downsampled to output resolution after simulation - internal_result_shape = ( - self._electrode_arrays[0].num_rows, - self._electrode_arrays[0].num_cols, - len(t_internal), - ) + # Result shapes follow each electrode array's own grid (arrays may differ). + def internal_result_shape(electrode_array: SurfaceElectrodeArray) -> tuple[int, int, int]: + return electrode_array.num_rows, electrode_array.num_cols, len(t_internal) - # Final output shape after downsampling - output_result_shape = ( - self._electrode_arrays[0].num_rows, - self._electrode_arrays[0].num_cols, - self._sampling_points_in_t_and_z_domains, - ) + def output_result_shape(electrode_array: SurfaceElectrodeArray) -> tuple[int, int, int]: + return ( + electrode_array.num_rows, + electrode_array.num_cols, + self._sampling_points_in_t_and_z_domains, + ) # Helper function to process a single motor unit def _process_single_mu( @@ -351,13 +347,17 @@ def _process_single_mu( electrode_array = electrode_array_original # Pre-allocated result array at INTERNAL resolution (optimization: use pre-computed shape) - array_result_internal = np.zeros(internal_result_shape, dtype=np.float64) + array_result_internal = np.zeros( + internal_result_shape(electrode_array), dtype=np.float64 + ) number_of_fibers = number_of_fibers_per_MUs[MU_index] if number_of_fibers == 0: # Return empty signal (downsampled to output resolution) - array_result_downsampled = np.zeros(output_result_shape, dtype=np.float64) + array_result_downsampled = np.zeros( + output_result_shape(electrode_array), dtype=np.float64 + ) return array_result_downsampled, f"MUAP_{MU_index}" # Get fiber positions @@ -428,9 +428,8 @@ def _process_single_mu( theta = theta_values[fiber_number] fiber_length__mm = fiber_lengths[fiber_number] - # Calculate fiber end positions - L1 = abs(innervation_zone + fiber_length__mm / 2) - L2 = abs(innervation_zone - fiber_length__mm / 2) + # Semi-lengths towards +z (L1) and -z (L2) from the end-plate + L1, L2 = fiber_semi_lengths(innervation_zone, fiber_length__mm) if self._use_unified: from myogen.simulator.core.emg.fiber_simulation import ( @@ -519,7 +518,7 @@ def _process_single_mu( f"Failed to process MU {MU_index} for electrode array {array_idx}: {e}" ) # Return empty signal with error marker at output resolution - empty_result = np.zeros(output_result_shape, dtype=np.float64) + empty_result = np.zeros(output_result_shape(electrode_array_original), dtype=np.float64) return empty_result, f"MUAP_{MU_index}_FAILED" block = Block() @@ -576,21 +575,13 @@ def _process_single_mu( IAP_SCALE_FACTOR = 2.5 kernel_length_mm = self._mean_fiber_length__mm * IAP_SCALE_FACTOR - # Physical duration based on IAP kernel length (after /=2 scaling in simulate_fiber) - # Duration = (kernel_length_mm / 2) / velocity - muap_duration__s = ( - kernel_length_mm / 2.0 - ) / self._mean_conduction_velocity__m_s / 1000.0 # Convert ms to s - - # Create custom time array for this duration - times__s = np.linspace(0, muap_duration__s, self._sampling_points_in_t_and_z_domains) - - # Calculate effective sampling rate for these times - effective_sampling_rate__Hz = ( - self._sampling_points_in_t_and_z_domains - 1 - ) / muap_duration__s - - use_custom_times = True + # The fibre window spans kernel_length / v (firing instant at its centre) + effective_sampling_rate__Hz = muap_sampling_rate__Hz( + self._mean_conduction_velocity__m_s, + self._internal_sampling_points, + self._sampling_points_in_t_and_z_domains, + kernel_length_mm, + ) # Create segments for ALL MUs (maintaining index order) # Non-computed MUs get empty signals at output resolution @@ -599,14 +590,13 @@ def _process_single_mu( array_result, segment_name = results[MU_index] else: # Create empty MUAP for non-computed MUs at output resolution - array_result = np.zeros(output_result_shape, dtype=np.float64) + array_result = np.zeros(output_result_shape(electrode_array), dtype=np.float64) segment_name = f"MUAP_{MU_index}" segment = Segment(name=segment_name) group.segments.append(segment) - # Use actual physical duration based on fiber lengths - grid_shape = (output_result_shape[0], output_result_shape[1]) + grid_shape = (electrode_array.num_rows, electrode_array.num_cols) segment.analogsignals.append( create_grid_signal( signal=np.transpose(array_result, (2, 0, 1)) * pq.mV, diff --git a/tests/test_fiber_simulation.py b/tests/test_fiber_simulation.py index f16cfc30..aa917acf 100644 --- a/tests/test_fiber_simulation.py +++ b/tests/test_fiber_simulation.py @@ -351,7 +351,7 @@ def test_muap_structural_properties(self, electrode_array): # different mathematical approaches: # - Old: frequency-domain convolution (FFT of source * H_vc, then IFFT) # - New: spatial kernel IFFT + time-domain traveling-wave convolution - # They also differ in Rosenfalck parameterization (D1=9.6 + z/2 vs D1=96). + # They also differ in the Rosenfalck amplitude constant (D1=9.6 vs D1=96). # Strict waveform correlation is not expected; instead we verify that both # produce physically plausible signals with consistent structural properties. @@ -392,12 +392,12 @@ def test_muap_structural_properties(self, electrode_array): ) # 4. The amplitude order of magnitude should be within a reasonable factor - # (allowing for the 10x D1 difference and z/2 scaling) + # (allowing for the 10x D1 difference) old_max = np.max(np.abs(phi_old)) new_max = np.max(np.abs(phi_new)) ratio = new_max / old_max print(f"Amplitude ratio (new/old): {ratio:.2f}") - # The D1 difference alone gives ~10x, plus z-scaling effects. + # The D1 difference alone gives ~10x. # Allow a wide range (0.1 to 10000) to accommodate the different formulations. assert 0.1 < ratio < 10000, ( f"Amplitude ratio {ratio:.1f} outside plausible range [0.1, 10000]" diff --git a/tests/test_surface_model_consistency.py b/tests/test_surface_model_consistency.py new file mode 100644 index 00000000..75a437ea --- /dev/null +++ b/tests/test_surface_model_consistency.py @@ -0,0 +1,189 @@ +"""Surface-EMG consistency with the Farina et al. (2004) model. + +Each test checks the analytical pipeline against an independent reference: +linearity of the detection system, the conduction velocity that was set, an +independent time-domain Rosenfalck path, or the documented source geometry. +""" +import numpy as np +import pytest +import quantities as pq +from scipy.signal import correlate, resample + +from myogen import set_random_seed +from myogen.simulator import ( + ExecutionPolicy, + Muscle, + RecruitmentThresholds, + SurfaceElectrodeArray, + SurfaceEMG, +) +from myogen.simulator.core.emg import fiber_simulation +from myogen.simulator.core.emg.surface import simulate_fiber +from myogen.utils.neo import signal_to_grid + +GEOMETRY = dict(r=8.5, r_bone=0.0, th_fat=0.3, th_skin=1.29, sig_muscle_rho=0.09, + sig_muscle_z=0.4, sig_fat=4.07e-2, sig_skin=4.88e-4) +V = 4.0 # mm/ms +N = 512 +KERNEL = 79.25 # mm +DT = KERNEL / (V * (N - 1)) # ms per internal sample + + +def _array(rows, cols, ied, z0, mode="monopolar"): + return SurfaceElectrodeArray(rows, cols, ied * pq.mm, 0.5 * pq.mm, + center_point__mm_deg=(z0 * pq.mm, 0 * pq.deg), + bending_radius__mm=GEOMETRY["r"] * pq.mm, + differentiation_mode=mode) + + +def _fibre(array, *, R=5.5, theta=0.15, L1=16.0, L2=16.0, zi=0.0, kernel=KERNEL, n=N): + phi, _, _ = simulate_fiber._simulate_fiber_v2_python( + Fs=10.0, v=V, N=n, M=32, R=R, L1=L1, L2=L2, zi=zi, electrode_array=array, + fiber_length__mm=kernel, use_gpu=False, theta_offset=-theta, **GEOMETRY) + return phi + + +def _upsampled_lag_ms(a, b, dt_ms, up=20): + a2, b2 = resample(a, len(a) * up), resample(b, len(b) * up) + return (np.argmax(correlate(b2, a2, "full")) - (len(a2) - 1)) * dt_ms / up + + +@pytest.fixture(params=["numba", "numpy"]) +def b_kz_backend(request, monkeypatch): + if request.param == "numpy": + monkeypatch.setattr(simulate_fiber, "_numba_B_kz_func", None) + monkeypatch.setattr(simulate_fiber, "_numba_phi_func", None) + return request.param + + +def test_longitudinal_bipolar_equals_difference_of_monopolar_pair(b_kz_backend): + bipolar = _fibre(_array(1, 1, 10, 8, "bipolar_longitudinal"))[0, 0] + pair = _fibre(_array(2, 1, 10, 8)) + expected = pair[1, 0] - pair[0, 0] + np.testing.assert_allclose(bipolar, expected, atol=1e-9 * np.ptp(expected)) + + +def test_transversal_bipolar_equals_difference_of_monopolar_pair(b_kz_backend): + bipolar = _fibre(_array(1, 1, 3, 8, "bipolar_transversal"))[0, 0] + pair = _fibre(_array(1, 2, 3, 8)) + expected = pair[0, 1] - pair[0, 0] + np.testing.assert_allclose(bipolar, expected, atol=1e-9 * np.ptp(expected)) + + +def test_separately_simulated_channels_keep_their_propagation_delay(): + # One call per electrode: a per-call re-centring would erase the delay. + near = _fibre(_array(1, 1, 5, 6), R=6.5, theta=0.0)[0, 0] + far = _fibre(_array(1, 1, 5, 11), R=6.5, theta=0.0)[0, 0] + assert _upsampled_lag_ms(near, far, DT) == pytest.approx(5.0 / V, rel=0.03) + + +def test_firing_instant_is_the_window_centre(): + # Monopolar negative peaks along a long fibre fall on t = t0 + z/v; the + # intercept t0 is the IAP depolarisation peak (z = 3 mm for Rosenfalck) + # leaving the end-plate, measured from the window centre. + n, kernel = 1024, 300.0 + dt = kernel / (V * (n - 1)) + z_values = np.array([10.0, 20.0, 30.0, 40.0]) + peaks = [np.argmin(_fibre(_array(1, 1, 5, z), R=6.7, theta=0.0, L1=60.0, L2=60.0, + kernel=kernel, n=n)[0, 0]) for z in z_values] + slope, intercept = np.polyfit(z_values, peaks, 1) + assert slope * dt == pytest.approx(1.0 / V, rel=0.01) + assert (intercept - n // 2) * dt * V == pytest.approx(3.0, abs=0.5) + + +def test_source_matches_independent_time_domain_rosenfalck(): + # simulate_fiber_hybrid builds the published Rosenfalck IAP (lambda = 1 mm) + # in the time domain; only the amplitude constant differs. + array = _array(1, 1, 5, 8) + kwargs = dict(Fs=10.0, v=V, N=N, M=32, R=6.5, L1=15.85, L2=15.85, zi=0.0, + electrode_array=array, fiber_length__mm=KERNEL, use_gpu=False, **GEOMETRY) + frequency = simulate_fiber._simulate_fiber_v2_python(**kwargs)[0][0, 0] + reference = fiber_simulation.simulate_fiber_hybrid(**kwargs, D1=96.0)[0][0, 0] + similarity = np.max(np.abs(correlate(frequency, reference, "full"))) + similarity /= np.linalg.norm(frequency) * np.linalg.norm(reference) + assert similarity > 0.99 + + +def test_l1_is_the_positive_z_semi_length(): + # Fibre spans [zi - L2, zi + L1]: pairs inside it see a travelling wave, + # pairs beyond its end see only a standing end-of-fibre component. + kernel = 100.0 + dt = kernel / (V * (N - 1)) + inside = _fibre(_array(2, 1, 5, 12.5), R=6.4, theta=0.0, L1=30.0, L2=8.0, kernel=kernel) + beyond = _fibre(_array(2, 1, 5, -12.5), R=6.4, theta=0.0, L1=30.0, L2=8.0, kernel=kernel) + assert _upsampled_lag_ms(inside[0, 0], inside[1, 0], dt) == pytest.approx(5.0 / V, rel=0.05) + assert abs(_upsampled_lag_ms(beyond[1, 0], beyond[0, 0], dt)) < 0.5 * 5.0 / V + + +def _muscle(**kwargs): + set_random_seed(42) + thresholds, _ = RecruitmentThresholds(N=20, recruitment_range__ratio=50, mode="konstantin") + muscle = Muscle(thresholds, fiber_density__fibers_per_mm2=10 * pq.mm**-2, **kwargs) + muscle.distribute_innervation_centers() + muscle.generate_muscle_fiber_centers(verbose=False) + muscle.assign_mfs2mns(verbose=False, n_jobs=1) + return muscle + + +@pytest.fixture(scope="module") +def muscle(): + return _muscle() + + +@pytest.fixture(scope="module") +def long_muscle(): + return _muscle(mean_fiber_length__mm=120 * pq.mm) + + +def _simulate(semg, route): + if route == "default": + return semg.simulate_muaps(n_jobs=1, verbose=False, use_gpu=False) + return semg.simulate_muaps(n_jobs=1, verbose=False, execution=ExecutionPolicy()) + + +@pytest.mark.parametrize("route", ["default", "execution"]) +@pytest.mark.usefixtures("ample_session_memory") +def test_muap_sources_span_the_fibre_symmetrically(route, muscle, monkeypatch): + from myogen.simulator.core.emg.runtime import semg as runtime_semg + from myogen.simulator.core.emg.surface import surface_emg + calls = [] + + def spy(**kwargs): + calls.append((kwargs["L1"], kwargs["L2"], kwargs["zi"])) + array = kwargs["electrode_array"] + return np.zeros((array.num_rows, array.num_cols, kwargs["N"])), None, None + + monkeypatch.setattr(surface_emg, "_simulate_fiber_v2_python", spy) + monkeypatch.setattr(runtime_semg, "_simulate_fiber_v2_python", spy) + _simulate(SurfaceEMG(muscle, [_array(1, 1, 5, 8)], MUs_to_simulate=[1, 2]), route) + L1, L2, zi = np.array(calls).T + assert np.all(zi != 0) + # Fibre of length L1 + L2 centred on z = 0 with its end-plate at zi. + np.testing.assert_allclose(zi + L1, (L1 + L2) / 2) + np.testing.assert_allclose(zi - L2, -(L1 + L2) / 2) + + +@pytest.mark.parametrize("route", ["default", "execution"]) +@pytest.mark.usefixtures("ample_session_memory") +def test_muap_time_axis_recovers_the_conduction_velocity(route, long_muscle): + muscle = long_muscle + mu = 1 + centre = muscle.innervation_center_positions__mm.magnitude[mu] + column = SurfaceElectrodeArray(5, 1, 5 * pq.mm, 0.5 * pq.mm, center_point__mm_deg=( + 30 * pq.mm, np.degrees(np.arctan2(centre[1], centre[0])) * pq.deg)) + signal = _simulate(SurfaceEMG(muscle, [column], MUs_to_simulate=[mu]), route + ).groups[0].segments[mu].analogsignals[0] + single_differential = np.diff(signal_to_grid(signal)[:, :, 0], axis=1) + dt_ms = 1e3 / float(signal.sampling_rate) + lags = [_upsampled_lag_ms(single_differential[:, i], single_differential[:, i + 1], dt_ms) + for i in range(single_differential.shape[1] - 1)] + velocity = float(muscle.mean_conduction_velocity__m_s.magnitude) + assert 5.0 / np.median(lags) == pytest.approx(velocity, rel=0.03) + + +def test_mixed_electrode_array_shapes_keep_their_own_shape(muscle): + semg = SurfaceEMG(muscle, [_array(2, 1, 10, 8), _array(1, 1, 10, 8, "bipolar_longitudinal")], + MUs_to_simulate=[1]) + block = semg.simulate_muaps(n_jobs=1, verbose=False, use_gpu=False) + shapes = [signal_to_grid(group.segments[1].analogsignals[0]).shape[1:] for group in block.groups] + assert shapes == [(2, 1), (1, 1)] From 4b71b4afed8db3c2f5979398a8601938d2175a23 Mon Sep 17 00:00:00 2001 From: Devon Rohlf Date: Fri, 30 Oct 2026 17:42:19 +0100 Subject: [PATCH 02/42] test(semg): cover the execution route's prepared phi contraction for bipolar --- tests/test_surface_model_consistency.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/tests/test_surface_model_consistency.py b/tests/test_surface_model_consistency.py index 75a437ea..de0e2342 100644 --- a/tests/test_surface_model_consistency.py +++ b/tests/test_surface_model_consistency.py @@ -18,6 +18,7 @@ SurfaceEMG, ) from myogen.simulator.core.emg import fiber_simulation +from myogen.simulator.core.emg.runtime.surface_preparation import SurfacePreparation from myogen.simulator.core.emg.surface import simulate_fiber from myogen.utils.neo import signal_to_grid @@ -36,10 +37,12 @@ def _array(rows, cols, ied, z0, mode="monopolar"): differentiation_mode=mode) -def _fibre(array, *, R=5.5, theta=0.15, L1=16.0, L2=16.0, zi=0.0, kernel=KERNEL, n=N): +def _fibre(array, *, R=5.5, theta=0.15, L1=16.0, L2=16.0, zi=0.0, kernel=KERNEL, n=N, + prepared=False): phi, _, _ = simulate_fiber._simulate_fiber_v2_python( Fs=10.0, v=V, N=n, M=32, R=R, L1=L1, L2=L2, zi=zi, electrode_array=array, - fiber_length__mm=kernel, use_gpu=False, theta_offset=-theta, **GEOMETRY) + fiber_length__mm=kernel, use_gpu=False, theta_offset=-theta, + _prepared=SurfacePreparation() if prepared else None, **GEOMETRY) return phi @@ -48,24 +51,25 @@ def _upsampled_lag_ms(a, b, dt_ms, up=20): return (np.argmax(correlate(b2, a2, "full")) - (len(a2) - 1)) * dt_ms / up -@pytest.fixture(params=["numba", "numpy"]) +@pytest.fixture(params=["numba", "numpy", "prepared"]) def b_kz_backend(request, monkeypatch): + # "prepared" is the execution route's phi contraction (SurfacePreparation). if request.param == "numpy": monkeypatch.setattr(simulate_fiber, "_numba_B_kz_func", None) monkeypatch.setattr(simulate_fiber, "_numba_phi_func", None) - return request.param + return request.param == "prepared" def test_longitudinal_bipolar_equals_difference_of_monopolar_pair(b_kz_backend): - bipolar = _fibre(_array(1, 1, 10, 8, "bipolar_longitudinal"))[0, 0] - pair = _fibre(_array(2, 1, 10, 8)) + bipolar = _fibre(_array(1, 1, 10, 8, "bipolar_longitudinal"), prepared=b_kz_backend)[0, 0] + pair = _fibre(_array(2, 1, 10, 8), prepared=b_kz_backend) expected = pair[1, 0] - pair[0, 0] np.testing.assert_allclose(bipolar, expected, atol=1e-9 * np.ptp(expected)) def test_transversal_bipolar_equals_difference_of_monopolar_pair(b_kz_backend): - bipolar = _fibre(_array(1, 1, 3, 8, "bipolar_transversal"))[0, 0] - pair = _fibre(_array(1, 2, 3, 8)) + bipolar = _fibre(_array(1, 1, 3, 8, "bipolar_transversal"), prepared=b_kz_backend)[0, 0] + pair = _fibre(_array(1, 2, 3, 8), prepared=b_kz_backend) expected = pair[0, 1] - pair[0, 0] np.testing.assert_allclose(bipolar, expected, atol=1e-9 * np.ptp(expected)) From abeec9751a1650c39c404dc7894f6bad96497b8f Mon Sep 17 00:00:00 2001 From: Devon Rohlf Date: Fri, 30 Oct 2026 18:41:04 +0100 Subject: [PATCH 03/42] fix(semg): continuous-transform scaling so amplitude no longer scales with 1/v The source FFT was normalised by 1/N and the inverse by N, dropping the dz (forward) and Fs/N (inverse) factors of the continuous transforms. Their product is dz * Fs = v, so every potential was divided by the conduction velocity. Use dz on the forward transform and Fs on the inverse. New test compares a single fibre against an exact, independently written Fourier-Bessel solution for a homogeneous insulated cylinder (both routes, two velocities, two geometries): now within 2 %, before it was off by 1/v. --- myogen/simulator/core/emg/runtime/semg.py | 2 +- .../core/emg/surface/simulate_fiber.py | 9 ++- tests/test_surface_amplitude.py | 74 +++++++++++++++++++ 3 files changed, 81 insertions(+), 4 deletions(-) create mode 100644 tests/test_surface_amplitude.py diff --git a/myogen/simulator/core/emg/runtime/semg.py b/myogen/simulator/core/emg/runtime/semg.py index 980deff2..05668fc5 100644 --- a/myogen/simulator/core/emg/runtime/semg.py +++ b/myogen/simulator/core/emg/runtime/semg.py @@ -71,7 +71,7 @@ def simulate_muaps(owner, policy, n_jobs, verbose, use_gpu): radius = float(electrode.electrode_radius__mm.rescale(pq.mm).magnitude) preparation = SurfacePreparation() dependencies = ( - "surface-muap-gitlab-fc077-v2", N, M, final_samples, kernel_length, rate, + "surface-muap-gitlab-fc077-v3", N, M, final_samples, kernel_length, rate, owner._internal_sampling_frequency__Hz, owner._mean_conduction_velocity__m_s, owner._radius_total, owner._radius_bone__mm, owner._fat_thickness__mm, owner._skin_thickness__mm, diff --git a/myogen/simulator/core/emg/surface/simulate_fiber.py b/myogen/simulator/core/emg/surface/simulate_fiber.py index bc4efba6..cdc5d862 100644 --- a/myogen/simulator/core/emg/surface/simulate_fiber.py +++ b/myogen/simulator/core/emg/surface/simulate_fiber.py @@ -429,7 +429,9 @@ def _simulate_fiber_v2_python( psi = -f_minus_t(aux) # print(f"DEBUG: psi range = [{np.min(psi):.6e}, {np.max(psi):.6e}]") - PSI = np.fft.fftshift(np.fft.fft(psi)) / len(psi) + # Continuous Fourier transform: integral psi(z) e^{-ikz} dz ~ dz * FFT. The + # inverse transforms in step 6 carry the matching dk_t / 2pi = Fs / N. + PSI = np.fft.fftshift(np.fft.fft(psi)) * float(z[1] - z[0]) PSI_conj = np.conj(PSI) PSI_conj = PSI_conj.reshape(-1, len(PSI_conj)) ones = np.ones((len(PSI_conj[0, :]), 1)) @@ -558,9 +560,10 @@ def _simulate_fiber_v2_python( # applied when building the centered spectrum. They are # equivalent for even N but differ by a one-sample phase ramp # for odd N — ifftshift is the correct inverse here. + # (1/2pi) integral dk_t e^{i k_t t} = (Fs / N) * sum = Fs * ifft phi[channel_z, channel_theta, :] = np.real( np.fft.ifft( - np.fft.ifftshift(PHI_complex[channel_z, channel_theta, :] * len(psi)) + np.fft.ifftshift(PHI_complex[channel_z, channel_theta, :] * Fs_effective) ) ) else: @@ -584,7 +587,7 @@ def _simulate_fiber_v2_python( # for odd N — ifftshift is the correct inverse here. phi[channel_z, channel_theta, :] = np.real( np.fft.ifft( - np.fft.ifftshift(PHI / 2 / math.pi * len(psi)) + np.fft.ifftshift(PHI / 2 / math.pi * Fs_effective) ) ) diff --git a/tests/test_surface_amplitude.py b/tests/test_surface_amplitude.py new file mode 100644 index 00000000..4f7e784c --- /dev/null +++ b/tests/test_surface_amplitude.py @@ -0,0 +1,74 @@ +"""Absolute sEMG amplitude against an exact volume-conductor solution. + +For a homogeneous, isotropic cylinder with an insulating surface the potential +of a point current has a closed Fourier-Bessel form (Jackson eq. 3.148 plus a +Neumann correction; Wronskian I_n K_n' - I_n' K_n = -1/x): + + phi(z, rho=r, theta) = sum_n int dk I_n(|k| rho0) / (4 pi^2 sigma |k| r I_n'(|k| r)) + * S(k) e^{ikz} e^{in theta} + +with S(k) the z-Fourier transform of the line current density. It is written +here independently of MyoGen's layered-conductor code. +""" +import math + +import numpy as np +import pytest +import quantities as pq +from scipy.special import ive + +from myogen.simulator.core.emg.electrodes import SurfaceElectrodeArray +from myogen.simulator.core.emg.runtime.surface_preparation import SurfacePreparation +from myogen.simulator.core.emg.surface import simulate_fiber + +SIGMA = 0.3 # S/m, every layer +SOURCE_SCALE = 96 / 10 # mV/mm^3, MyoGen's Rosenfalck dVm/dz amplitude +DU = 0.02 # mm + + +def exact_surface_profile(i_line, du, rho0, r, sigma, n_max=200): + """Surface potential along z directly above a line source sampled at spacing du.""" + n = len(i_line) + k = 2 * math.pi * np.fft.fftfreq(n, d=du) + S = np.fft.fft(i_line) * du + x0, x1 = np.abs(k[1:]) * rho0, np.abs(k[1:]) * r + total = np.zeros(n - 1) + for order in range(-n_max, n_max + 1): + m = abs(order) + dI = 0.5 * (ive(m - 1, x1) + ive(m + 1, x1)) # I_m'(x) e^{-x} + with np.errstate(divide="ignore", invalid="ignore"): + term = ive(m, x0) / dI * np.exp(x0 - x1) + total += np.where(np.isfinite(term), term, 0.0) # high orders underflow; they vanish + g = np.concatenate([[0.0], total / (4 * math.pi**2 * sigma * x1)]) + return np.real(np.fft.ifft(S * g)) * n * (2 * math.pi / (n * du)) + + +def rosenfalck_line_source(u): + """i(z) = d2Vm/dz2 * SOURCE_SCALE / A for a front at u = 0 travelling towards +z.""" + s = -u + dvm = np.where(s >= 0, SOURCE_SCALE * np.exp(-s) * (3 * s**2 - s**3), 0.0) + return -np.gradient(dvm, u[1] - u[0]) + + +def _myogen_passing_wave(v, r, depth, prepared, z0=30.0, n=1024, kernel=300.0): + array = SurfaceElectrodeArray(1, 1, 5 * pq.mm, 0.01 * pq.mm, center_point__mm_deg=(z0 * pq.mm, 0 * pq.deg), + bending_radius__mm=r * pq.mm, differentiation_mode="monopolar") + phi, _, _ = simulate_fiber._simulate_fiber_v2_python( + Fs=10.0, v=v, N=n, M=64, r=r, r_bone=0.0, th_fat=0.05, th_skin=0.05, R=r - depth, L1=140.0, L2=140.0, + zi=0.0, electrode_array=array, sig_muscle_rho=SIGMA, sig_muscle_z=SIGMA, sig_fat=SIGMA, sig_skin=SIGMA, + fiber_length__mm=kernel, use_gpu=False, theta_offset=0.0, + _prepared=SurfacePreparation() if prepared else None) + t = (np.arange(n) - n // 2) * kernel / (v * (n - 1)) # ms from the firing instant + window = (t > (z0 - 12) / v) & (t < (z0 + 12) / v) # the wave passing the electrode only + return phi[0, 0][window] + + +@pytest.mark.parametrize("prepared", [False, True], ids=["default", "prepared"]) +@pytest.mark.parametrize("v", [2.0, 4.0]) +@pytest.mark.parametrize("r, depth", [(10.0, 2.0), (20.0, 3.0)]) +def test_single_fibre_amplitude_matches_the_exact_cylinder_solution(r, depth, v, prepared): + u = np.arange(-60, 60, DU) + exact = exact_surface_profile(rosenfalck_line_source(u), DU, r - depth, r, SIGMA) + simulated = _myogen_passing_wave(v, r, depth, prepared) + # Amplitude must not depend on conduction velocity: the spatial wave is the same. + assert np.ptp(simulated) == pytest.approx(np.ptp(exact), rel=0.02) From 60df7c2e9df506c60ef54e6073418e7c7440a9ad Mon Sep 17 00:00:00 2001 From: Devon Rohlf Date: Fri, 30 Oct 2026 18:49:41 +0100 Subject: [PATCH 04/42] feat(semg)!: physical source current sigma_i * pi a^2 * d2Vm/dz2 The fibre source was the bare Rosenfalck slope with an ad-hoc A = 96/10 and no fibre properties, so outputs labelled mV were not potentials. Use the published A = 96 mV/mm^3 and the core-conductor line current density sigma_i * pi * (d/2)^2 * d2Vm/dz2; with S/m and mm the result is in mV. SurfaceEMG gains fiber_diameter__um (default 50) and intracellular_conductivity__S_m (default 1.01), passed to every fibre on both routes and included in the execution cache key. The experimental hybrid path (use_unified=True) keeps its own D1 source and is unchanged. Tests: exact-solution comparison now uses the physical source; amplitude scales with d^2 and sigma_i; both routes forward the parameters. --- myogen/simulator/core/emg/runtime/semg.py | 5 +- .../core/emg/surface/simulate_fiber.py | 35 +++++++++-- .../simulator/core/emg/surface/surface_emg.py | 14 +++++ tests/test_surface_amplitude.py | 58 ++++++++++++++++--- 4 files changed, 98 insertions(+), 14 deletions(-) diff --git a/myogen/simulator/core/emg/runtime/semg.py b/myogen/simulator/core/emg/runtime/semg.py index 05668fc5..e0a7aeb7 100644 --- a/myogen/simulator/core/emg/runtime/semg.py +++ b/myogen/simulator/core/emg/runtime/semg.py @@ -71,12 +71,13 @@ def simulate_muaps(owner, policy, n_jobs, verbose, use_gpu): radius = float(electrode.electrode_radius__mm.rescale(pq.mm).magnitude) preparation = SurfacePreparation() dependencies = ( - "surface-muap-gitlab-fc077-v3", N, M, final_samples, kernel_length, rate, + "surface-muap-gitlab-fc077-v4", N, M, final_samples, kernel_length, rate, owner._internal_sampling_frequency__Hz, owner._mean_conduction_velocity__m_s, owner._radius_total, owner._radius_bone__mm, owner._fat_thickness__mm, owner._skin_thickness__mm, owner._muscle_conductivity_radial__S_m, owner._muscle_conductivity_longitudinal__S_m, owner._fat_conductivity__S_m, owner._skin_conductivity__S_m, + owner._fiber_diameter__mm, owner._intracellular_conductivity__S_m, vars(electrode), base_z, base_theta, radius, ) results = {} @@ -116,6 +117,8 @@ def compute(start, stop): theta_offset=-angular[fibre], pos_z_precomputed=base_z, pos_theta_precomputed=base_theta, rele_precomputed=radius, _prepared=preparation, + sig_intracellular=owner._intracellular_conductivity__S_m, + fiber_diameter__mm=owner._fiber_diameter__mm, ) if not np.all(np.isfinite(values[local])): raise FloatingPointError(f"Nonfinite surface fibre in MU {mu_index}") diff --git a/myogen/simulator/core/emg/surface/simulate_fiber.py b/myogen/simulator/core/emg/surface/simulate_fiber.py index cdc5d862..5ccec8e9 100644 --- a/myogen/simulator/core/emg/surface/simulate_fiber.py +++ b/myogen/simulator/core/emg/surface/simulate_fiber.py @@ -266,6 +266,9 @@ def _simulate_fiber_v2_python( pos_theta_precomputed: np.ndarray | None = None, rele_precomputed: float | None = None, _prepared=None, + *, + sig_intracellular: float = 1.01, + fiber_diameter__mm: float = 0.05, ): """ Simulate a single fiber (Python implementation). @@ -330,6 +333,12 @@ def _simulate_fiber_v2_python( instant at sample ``N // 2``. See :func:`muap_sampling_rate__Hz`. Example: fiber_length__mm=100 with v=4 m/s → 25 ms window. + sig_intracellular : float, optional + Intracellular conductivity (S/m), default 1.01. + fiber_diameter__mm : float, optional + Fibre diameter (mm), default 0.05. Together with ``sig_intracellular`` it + sets the line current density ``sig_i * pi * (d/2)**2 * d2Vm/dz2`` of the + source, so the output is the potential in mV. Returns ------- @@ -359,7 +368,8 @@ def _simulate_fiber_v2_python( if _prepared is not None: _prepared.bind((Fs, v, N, M, r, r_bone, th_fat, th_skin, - sig_muscle_rho, sig_muscle_z, sig_fat, sig_skin, sig_bone, fiber_length__mm)) + sig_muscle_rho, sig_muscle_z, sig_fat, sig_skin, sig_bone, fiber_length__mm, + sig_intracellular, fiber_diameter__mm)) if _prepared is None or _prepared.grid is None: ## Model angular frequencies # FFT-consistent bin centers (matches np.fft.fft + fftshift): spacing = 1, bins at @@ -401,10 +411,12 @@ def _simulate_fiber_v2_python( ################################################################################################### ## 2. I(k_t, k_z) - # Rosenfalck (1969) intracellular action potential spatial distribution - # As used in Farina & Merletti (2001), equation 16 - # IMPORTANT: This matches the original Farina implementation exactly (no normalization) - A = 96 / 10 # mV/mm^3 (amplitude parameter from Farina 2001, eq 16) + # Rosenfalck (1969) intracellular action potential, Vm = 96 z^3 e^{-z} - 90 mV + # (Farina & Merletti 2001, eq. 16); aux below is dVm/dz in mV/mm. + A = 96 # mV/mm^3 + # Line current density i = sig_i * pi * (d/2)^2 * d2Vm/dz2 (core-conductor + # model). With conductivities in S/m and lengths in mm, phi comes out in mV. + source_scale = sig_intracellular * math.pi * (fiber_diameter__mm / 2) ** 2 # Define spatial grid for IAP kernel evaluation if fiber_length__mm is None: @@ -421,7 +433,7 @@ def _simulate_fiber_v2_python( aux = np.zeros_like(z) for i in range(len(z)): if z[i] >= 0: - aux[i] = A * np.exp(-z[i]) * (3 * z[i] ** 2 - z[i] ** 3) + aux[i] = source_scale * A * np.exp(-z[i]) * (3 * z[i] ** 2 - z[i] ** 3) # print(f"DEBUG: z range = [{np.min(z):.2f}, {np.max(z):.2f}] mm") # print(f"DEBUG: aux range = [{np.min(aux):.6e}, {np.max(aux):.6e}]") @@ -630,6 +642,9 @@ def simulate_fiber_v2( pos_z_precomputed: np.ndarray | None = None, pos_theta_precomputed: np.ndarray | None = None, rele_precomputed: float | None = None, + *, + sig_intracellular: float = 1.01, + fiber_diameter__mm: float = 0.05, ): """ Simulate a single fiber in the layered-cylinder volume conductor. @@ -700,6 +715,12 @@ def simulate_fiber_v2( instant at sample ``N // 2``. See :func:`muap_sampling_rate__Hz`. Example: fiber_length__mm=100 with v=4 m/s → 25 ms window. + sig_intracellular : float, optional + Intracellular conductivity (S/m), default 1.01. + fiber_diameter__mm : float, optional + Fibre diameter (mm), default 0.05. Together with ``sig_intracellular`` it + sets the line current density ``sig_i * pi * (d/2)**2 * d2Vm/dz2`` of the + source, so the output is the potential in mV. Returns ------- @@ -755,4 +776,6 @@ def simulate_fiber_v2( pos_z_precomputed, pos_theta_precomputed, rele_precomputed, + sig_intracellular=sig_intracellular, + fiber_diameter__mm=fiber_diameter__mm, ) diff --git a/myogen/simulator/core/emg/surface/surface_emg.py b/myogen/simulator/core/emg/surface/surface_emg.py index b88aca86..af036f1f 100644 --- a/myogen/simulator/core/emg/surface/surface_emg.py +++ b/myogen/simulator/core/emg/surface/surface_emg.py @@ -84,6 +84,12 @@ class SurfaceEMG: centre. ``K`` must hold the whole fibre plus the IAP, or the signal wraps around; it does not change the physical MUAP, only the window and the internal sampling step ``K / (v (N - 1))``. + fiber_diameter__um : float, optional + Muscle fibre diameter in µm, default 50. The source current scales with + the fibre cross-section (``d**2``). + intracellular_conductivity__S_m : float, optional + Intracellular conductivity in S/m, default 1.01. The source current + density is ``sigma_i * pi * (d/2)**2 * d2Vm/dz2`` (core-conductor model). Attributes ---------- @@ -112,6 +118,8 @@ def __init__( internal_sampling_frequency__Hz: Quantity__Hz | None = None, iap_kernel_length__mm: float | None = None, use_unified: bool = False, + fiber_diameter__um: float = 50.0, + intracellular_conductivity__S_m: float = 1.01, ): # Immutable public arguments - never modify these self.muscle_model = muscle_model @@ -122,6 +130,8 @@ def __init__( self.MUs_to_simulate = MUs_to_simulate self.iap_kernel_length__mm = iap_kernel_length__mm self._use_unified = use_unified + self.fiber_diameter__um = fiber_diameter__um + self.intracellular_conductivity__S_m = intracellular_conductivity__S_m # Internal sampling frequency for higher resolution MUAP computation # If not specified, defaults to 10 kHz for better MUAP resolution @@ -149,6 +159,8 @@ def __init__( self._sampling_points_in_theta_domain = sampling_points_in_theta_domain self._MUs_to_simulate = MUs_to_simulate self._iap_kernel_length__mm = iap_kernel_length__mm + self._fiber_diameter__mm = float(fiber_diameter__um) * 1e-3 + self._intracellular_conductivity__S_m = float(intracellular_conductivity__S_m) # Derived properties from muscle model - immutable public access self.mean_conduction_velocity__m_s = self._muscle_model.mean_conduction_velocity__m_s @@ -500,6 +512,8 @@ def _process_single_mu( pos_z_precomputed=base_pos_z, pos_theta_precomputed=base_pos_theta, rele_precomputed=base_rele, + sig_intracellular=self._intracellular_conductivity__S_m, + fiber_diameter__mm=self._fiber_diameter__mm, ) array_result_internal += phi_temp diff --git a/tests/test_surface_amplitude.py b/tests/test_surface_amplitude.py index 4f7e784c..b464d18b 100644 --- a/tests/test_surface_amplitude.py +++ b/tests/test_surface_amplitude.py @@ -22,7 +22,8 @@ from myogen.simulator.core.emg.surface import simulate_fiber SIGMA = 0.3 # S/m, every layer -SOURCE_SCALE = 96 / 10 # mV/mm^3, MyoGen's Rosenfalck dVm/dz amplitude +SIGMA_I = 1.01 # S/m, intracellular conductivity (default) +DIAMETER = 0.05 # mm, fibre diameter (default) DU = 0.02 # mm @@ -43,21 +44,24 @@ def exact_surface_profile(i_line, du, rho0, r, sigma, n_max=200): return np.real(np.fft.ifft(S * g)) * n * (2 * math.pi / (n * du)) -def rosenfalck_line_source(u): - """i(z) = d2Vm/dz2 * SOURCE_SCALE / A for a front at u = 0 travelling towards +z.""" +def rosenfalck_line_source(u, sigma_i=SIGMA_I, diameter=DIAMETER): + """Line current density sigma_i * pi a^2 * d2Vm/dz2 of a front at u = 0 travelling towards +z. + + Vm = 96 z^3 e^{-z} - 90 mV (Rosenfalck 1969), z in mm behind the front. + """ s = -u - dvm = np.where(s >= 0, SOURCE_SCALE * np.exp(-s) * (3 * s**2 - s**3), 0.0) - return -np.gradient(dvm, u[1] - u[0]) + dvm = np.where(s >= 0, 96 * np.exp(-s) * (3 * s**2 - s**3), 0.0) + return sigma_i * math.pi * (diameter / 2) ** 2 * -np.gradient(dvm, u[1] - u[0]) -def _myogen_passing_wave(v, r, depth, prepared, z0=30.0, n=1024, kernel=300.0): +def _myogen_passing_wave(v, r, depth, prepared, z0=30.0, n=1024, kernel=300.0, **fibre): array = SurfaceElectrodeArray(1, 1, 5 * pq.mm, 0.01 * pq.mm, center_point__mm_deg=(z0 * pq.mm, 0 * pq.deg), bending_radius__mm=r * pq.mm, differentiation_mode="monopolar") phi, _, _ = simulate_fiber._simulate_fiber_v2_python( Fs=10.0, v=v, N=n, M=64, r=r, r_bone=0.0, th_fat=0.05, th_skin=0.05, R=r - depth, L1=140.0, L2=140.0, zi=0.0, electrode_array=array, sig_muscle_rho=SIGMA, sig_muscle_z=SIGMA, sig_fat=SIGMA, sig_skin=SIGMA, fiber_length__mm=kernel, use_gpu=False, theta_offset=0.0, - _prepared=SurfacePreparation() if prepared else None) + _prepared=SurfacePreparation() if prepared else None, **fibre) t = (np.arange(n) - n // 2) * kernel / (v * (n - 1)) # ms from the firing instant window = (t > (z0 - 12) / v) & (t < (z0 + 12) / v) # the wave passing the electrode only return phi[0, 0][window] @@ -72,3 +76,43 @@ def test_single_fibre_amplitude_matches_the_exact_cylinder_solution(r, depth, v, simulated = _myogen_passing_wave(v, r, depth, prepared) # Amplitude must not depend on conduction velocity: the spatial wave is the same. assert np.ptp(simulated) == pytest.approx(np.ptp(exact), rel=0.02) + + +@pytest.mark.parametrize("fibre, scale", [(dict(fiber_diameter__mm=0.08), (0.08 / DIAMETER) ** 2), + (dict(sig_intracellular=0.5), 0.5 / SIGMA_I)]) +def test_source_scales_with_fibre_cross_section_and_intracellular_conductivity(fibre, scale): + reference = np.ptp(_myogen_passing_wave(4.0, 10.0, 2.0, False)) + assert np.ptp(_myogen_passing_wave(4.0, 10.0, 2.0, False, **fibre)) == pytest.approx(scale * reference, rel=1e-6) + + +@pytest.mark.parametrize("route", ["default", "execution"]) +@pytest.mark.usefixtures("ample_session_memory") +def test_surface_emg_passes_fibre_properties_to_every_fibre(route, monkeypatch): + from myogen import set_random_seed + from myogen.simulator import ExecutionPolicy, Muscle, RecruitmentThresholds, SurfaceEMG + from myogen.simulator.core.emg.runtime import semg as runtime_semg + from myogen.simulator.core.emg.surface import surface_emg + + set_random_seed(42) + thresholds, _ = RecruitmentThresholds(N=10, recruitment_range__ratio=20, mode="konstantin") + muscle = Muscle(thresholds, fiber_density__fibers_per_mm2=5 * pq.mm**-2) + muscle.distribute_innervation_centers() + muscle.generate_muscle_fiber_centers(verbose=False) + muscle.assign_mfs2mns(verbose=False, n_jobs=1) + calls = [] + + def spy(**kwargs): + calls.append((kwargs["fiber_diameter__mm"], kwargs["sig_intracellular"])) + array = kwargs["electrode_array"] + return np.zeros((array.num_rows, array.num_cols, kwargs["N"])), None, None + + monkeypatch.setattr(surface_emg, "_simulate_fiber_v2_python", spy) + monkeypatch.setattr(runtime_semg, "_simulate_fiber_v2_python", spy) + array = SurfaceElectrodeArray(1, 1, 5 * pq.mm, 0.5 * pq.mm, center_point__mm_deg=(8 * pq.mm, 0 * pq.deg)) + semg = SurfaceEMG(muscle, [array], MUs_to_simulate=[1, 2], fiber_diameter__um=60.0, + intracellular_conductivity__S_m=0.9) + if route == "default": + semg.simulate_muaps(n_jobs=1, verbose=False, use_gpu=False) + else: + semg.simulate_muaps(n_jobs=1, verbose=False, execution=ExecutionPolicy()) + assert calls and all(d == pytest.approx(0.06) and s == pytest.approx(0.9) for d, s in calls) From a8ff287284a5a47b729cb51029106537f882e556 Mon Sep 17 00:00:00 2001 From: Devon Rohlf Date: Fri, 30 Oct 2026 18:56:13 +0100 Subject: [PATCH 05/42] fix(semg)!: default skin conductivity 1 S/m (Farina 2004, from Roeleveld 1997) The previous default, 4.88e-4 S/m (measured composite skin, Botelho et al. 2019, tissue data at 150 Hz), raised surface amplitudes 3.5-6x. 1 S/m is the value used by Farina et al. (2004), taken from Roeleveld et al. (1997), who chose it so their layered cylinder matched how biceps MUAPs fall off with depth and around the arm. It is an effective parameter of the cylinder model, not a measured skin conductivity. Pass 4.88e-4 explicitly to reproduce the old default. The execution cache key already includes the skin conductivity. CHANGELOG: amplitude fixes, new SurfaceEMG fibre parameters, upgrade notes. --- CHANGELOG.md | 10 ++++++++++ myogen/simulator/core/muscle/muscle.py | 9 ++++++--- tests/test_surface_amplitude.py | 10 ++++++++++ 3 files changed, 26 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 931a6134..b9aba031 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,10 +17,20 @@ This release corrects the surface-EMG pipeline against Farina, Mesin, Martina & - **Fibre semi-lengths.** `L1` (towards +z) and `L2` were swapped, which moved both tendons by 2·z_i with the end-plate. Fibres now span `±fiber_length / 2` with the end-plate at z_i (`fiber_semi_lengths`). - **Fibre timing within a motor unit.** Each fibre's potential was circularly shifted onto its own peak before summation, misaligning fibres (up to ~60 % MUAP error on bipolar channels). All fibres now share the firing instant at the window centre. - **Mixed electrode-array shapes.** On the default route each array's MUAPs have that array's own grid shape; previously every array was broadcast into the first array's shape. +- **Absolute sEMG amplitude.** Three scaling errors are corrected; together they made MUAPs labelled mV roughly 50-100× too large. A single fibre now matches an exact solution for a homogeneous insulated cylinder within 2 %. + - The Fourier transforms dropped their `dz` and `Fs/N` factors, so every potential was divided by the conduction velocity. + - The source was the bare Rosenfalck slope with an ad-hoc `A = 96/10`. It is now the core-conductor line current `σᵢ·π·(d/2)²·∂²Vm/∂z²` with the published `A = 96 mV/mm³`. + - The default skin conductivity was 4.88·10⁻⁴ S/m (dry skin), which raised amplitudes 3.5-6×. It is now 1 S/m, as in Farina et al. (2004). + +### Added + +- `SurfaceEMG(fiber_diameter__um=50, intracellular_conductivity__S_m=1.01)` set the source strength; amplitude scales with the fibre cross-section. ### Upgrade notes - 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 dry-skin conductor, pass `skin_conductivity__S_per_m=4.88e-4 * pq.S / pq.m` to `Muscle`. +- 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. ## [0.13.0] - 2026-10-28 diff --git a/myogen/simulator/core/muscle/muscle.py b/myogen/simulator/core/muscle/muscle.py index f7c5f507..2e0bd91f 100644 --- a/myogen/simulator/core/muscle/muscle.py +++ b/myogen/simulator/core/muscle/muscle.py @@ -129,8 +129,11 @@ class Muscle: Muscle conductivity in longitudinal direction. Default is set to 0.4 S/m as determined by Botelho et al. 2019 [7] (Table 1). fat_conductivity__S_per_m : float, default=4.07E-2 Fat conductivity. Default is set to 4.07E-2 S/m as determined by Botelho et al. 2019 [7] (Table 1). - skin_conductivity__S_per_m : float, default=4.88E-4 - Skin conductivity. Default is set to 4.88E-4 S/m as determined by Botelho et al. 2019 [7] (Table 1). + skin_conductivity__S_per_m : float, default=1.0 + Skin conductivity. Default is 1 S/m as in Farina et al. 2004 [1]: gel electrodes bypass the + dry stratum corneum. Botelho et al. 2019 [7] used 4.88E-4 S/m (tissue data at 150 Hz, the + dry-skin range), which raises surface amplitudes 3.5-6x; pass it explicitly to model dry + electrodes. grid_resolution : int, default=256 Resolution of the computational grid used for innervation the muscle. Higher values provide more accurate spatial distribution but increase computational cost. @@ -192,7 +195,7 @@ def __init__( muscle_conductivity_radial__S_per_m: Quantity__S_per_m = 0.09 * pq.S / pq.m, muscle_conductivity_longitudinal__S_per_m: Quantity__S_per_m = 0.4 * pq.S / pq.m, fat_conductivity__S_per_m: Quantity__S_per_m = 4.07e-2 * pq.S / pq.m, - skin_conductivity__S_per_m: Quantity__S_per_m = 4.88e-4 * pq.S / pq.m, + skin_conductivity__S_per_m: Quantity__S_per_m = 1.0 * pq.S / pq.m, grid_resolution: int = 256, autorun: bool = False, ) -> None: diff --git a/tests/test_surface_amplitude.py b/tests/test_surface_amplitude.py index b464d18b..4b89be38 100644 --- a/tests/test_surface_amplitude.py +++ b/tests/test_surface_amplitude.py @@ -116,3 +116,13 @@ def spy(**kwargs): else: semg.simulate_muaps(n_jobs=1, verbose=False, execution=ExecutionPolicy()) assert calls and all(d == pytest.approx(0.06) and s == pytest.approx(0.9) for d, s in calls) + + +def test_default_skin_conductivity_is_hydrated_skin(): + # Gel electrodes bypass the dry stratum corneum; Farina et al. (2004) use 1 S/m. + # The previous default (4.88e-4 S/m, dry skin) raised surface amplitudes 3.5-6x. + from myogen.simulator import Muscle, RecruitmentThresholds + + thresholds, _ = RecruitmentThresholds(N=5, recruitment_range__ratio=10, mode="konstantin") + muscle = Muscle(thresholds) + assert float(muscle.skin_conductivity__S_m.rescale(pq.S / pq.m).magnitude) == pytest.approx(1.0) From 2380a068f0938446f6455eba2d268fd6b85b3288 Mon Sep 17 00:00:00 2001 From: Devon Rohlf Date: Fri, 30 Oct 2026 19:00:30 +0100 Subject: [PATCH 06/42] test(semg): equal-conductivity layers and bone reduce to the exact homogeneous solution --- tests/test_surface_amplitude.py | 22 +++++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/tests/test_surface_amplitude.py b/tests/test_surface_amplitude.py index 4b89be38..e9c0fe43 100644 --- a/tests/test_surface_amplitude.py +++ b/tests/test_surface_amplitude.py @@ -54,12 +54,14 @@ def rosenfalck_line_source(u, sigma_i=SIGMA_I, diameter=DIAMETER): return sigma_i * math.pi * (diameter / 2) ** 2 * -np.gradient(dvm, u[1] - u[0]) -def _myogen_passing_wave(v, r, depth, prepared, z0=30.0, n=1024, kernel=300.0, **fibre): +def _myogen_passing_wave(v, r, depth, prepared, z0=30.0, n=1024, kernel=300.0, th_fat=0.05, th_skin=0.05, + r_bone=0.0, **fibre): array = SurfaceElectrodeArray(1, 1, 5 * pq.mm, 0.01 * pq.mm, center_point__mm_deg=(z0 * pq.mm, 0 * pq.deg), bending_radius__mm=r * pq.mm, differentiation_mode="monopolar") phi, _, _ = simulate_fiber._simulate_fiber_v2_python( - Fs=10.0, v=v, N=n, M=64, r=r, r_bone=0.0, th_fat=0.05, th_skin=0.05, R=r - depth, L1=140.0, L2=140.0, - zi=0.0, electrode_array=array, sig_muscle_rho=SIGMA, sig_muscle_z=SIGMA, sig_fat=SIGMA, sig_skin=SIGMA, + Fs=10.0, v=v, N=n, M=64, r=r, r_bone=r_bone, th_fat=th_fat, th_skin=th_skin, R=r - depth, L1=140.0, + L2=140.0, zi=0.0, electrode_array=array, sig_muscle_rho=SIGMA, sig_muscle_z=SIGMA, sig_fat=SIGMA, + sig_skin=SIGMA, sig_bone=SIGMA, fiber_length__mm=kernel, use_gpu=False, theta_offset=0.0, _prepared=SurfacePreparation() if prepared else None, **fibre) t = (np.arange(n) - n // 2) * kernel / (v * (n - 1)) # ms from the firing instant @@ -78,6 +80,20 @@ def test_single_fibre_amplitude_matches_the_exact_cylinder_solution(r, depth, v, assert np.ptp(simulated) == pytest.approx(np.ptp(exact), rel=0.02) +@pytest.mark.parametrize("prepared", [False, True], ids=["default", "prepared"]) +@pytest.mark.parametrize("layers", [dict(th_fat=1.5, th_skin=1.0), dict(r_bone=1.0), + dict(th_fat=1.5, th_skin=1.0, r_bone=1.0)], + ids=["thick-fat-skin", "bone", "both"]) +def test_interfaces_between_equal_conductivities_are_transparent(layers, prepared): + # Production geometries use real layer thicknesses and a bone; with every + # conductivity equal they must reduce to the same exact homogeneous solution. + r, depth = 12.0, 5.0 + u = np.arange(-60, 60, DU) + exact = exact_surface_profile(rosenfalck_line_source(u), DU, r - depth, r, SIGMA) + simulated = _myogen_passing_wave(4.0, r, depth, prepared, **layers) + assert np.ptp(simulated) == pytest.approx(np.ptp(exact), rel=0.02) + + @pytest.mark.parametrize("fibre, scale", [(dict(fiber_diameter__mm=0.08), (0.08 / DIAMETER) ** 2), (dict(sig_intracellular=0.5), 0.5 / SIGMA_I)]) def test_source_scales_with_fibre_cross_section_and_intracellular_conductivity(fibre, scale): From 5ce91755538fb9d6eee2f746fa4c79a0d17995a7 Mon Sep 17 00:00:00 2001 From: Devon Rohlf Date: Sat, 31 Oct 2026 01:08:01 +0100 Subject: [PATCH 07/42] refactor(emg): move the two-layer nerve model into a shared module nmj_branch_schedule, two_layer_nmj and nerve_delays__ms hold the Konstantin et al. (2020) branch layout, end-plate positions and axon path lengths that MotorUnitSim used inline. Draw order is unchanged, so intramuscular EMG reproduces its previous end-plates and paths exactly. --- .../emg/intramuscular/intramuscular_emg.py | 17 ++--- .../core/emg/intramuscular/motor_unit_sim.py | 66 ++++------------- myogen/simulator/core/emg/nerve_branches.py | 73 +++++++++++++++++++ 3 files changed, 92 insertions(+), 64 deletions(-) create mode 100644 myogen/simulator/core/emg/nerve_branches.py diff --git a/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py b/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py index 48fc34d3..128bfd2d 100644 --- a/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py +++ b/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py @@ -38,6 +38,7 @@ except (ImportError, Exception): HAS_CUPY = False +from ..nerve_branches import nmj_branch_schedule from .motor_unit_sim import MotorUnitSim @@ -273,27 +274,19 @@ def _simulate_neuromuscular_junctions(self, verbose: bool = True) -> None: if not self._motor_units: raise ValueError("Must call _initialize_motor_units() first") - n_branches = 1 + np.round( - np.log( - self._muscle_model.recruitment_thresholds - / self._muscle_model.recruitment_thresholds[0] - ) + n_branches, branches_z_std, arborization_z_std = nmj_branch_schedule( + self._muscle_model.recruitment_thresholds ) for mu_idx, mu_sim in enumerate( tqdm(self._motor_units, desc="Setting up NMJ distributions", unit="Simulator", disable=not verbose) ): - spread_factor = np.sum(self._muscle_model.recruitment_thresholds[:mu_idx]) / np.sum( - self._muscle_model.recruitment_thresholds - ) - - # Create NMJ distribution # Branch spread increases with motor unit size mu_sim.sim_nmj_branches_two_layers( n_branches=int(n_branches[mu_idx]), endplate_center=self._endplate_center__mm, - branches_z_std=1.5 + spread_factor * 4.0, - arborization_z_std=0.5 + spread_factor * 1.5, + branches_z_std=float(branches_z_std[mu_idx]), + arborization_z_std=float(arborization_z_std[mu_idx]), ) def _calculate_muaps(self, n_jobs: int = -2, verbose: bool = True) -> INTRAMUSCULAR_MUAP__Block: diff --git a/myogen/simulator/core/emg/intramuscular/motor_unit_sim.py b/myogen/simulator/core/emg/intramuscular/motor_unit_sim.py index 694828ca..86ab3af8 100644 --- a/myogen/simulator/core/emg/intramuscular/motor_unit_sim.py +++ b/myogen/simulator/core/emg/intramuscular/motor_unit_sim.py @@ -13,7 +13,6 @@ from typing import Optional, List import numpy as np -from sklearn.cluster import KMeans from tqdm import tqdm try: @@ -25,6 +24,7 @@ from myogen import derive_subseed, get_random_generator from myogen.utils.decorators import beartowertype +from ..nerve_branches import two_layer_nmj from .bioelectric import ( get_current_density, get_elementary_current_response, @@ -139,60 +139,22 @@ def sim_nmj_branches_two_layers( arborization_z_std : float Standard deviation of secondary arborization in mm """ - rng = get_random_generator() - - self.nerve_paths = np.zeros( - (self._number_of_muscle_fibers, 2) - ) # Point coordinates - - kmeans = KMeans( - n_clusters=n_branches, init="k-means++", max_iter=100, random_state=derive_subseed(n_branches) - ) - idx = kmeans.fit_predict(self.muscle_fiber_centers__mm) - c = kmeans.cluster_centers_ - - self.branch_points_xy = c - self.branch_points_z = endplate_center + branches_z_std * rng.standard_normal( - size=(n_branches, 1) - ) - - self._neuromuscular_z_coordinates__mm = np.array( - [ - self.branch_points_z[idx[i]] - + arborization_z_std * rng.standard_normal() - for i in range(self._number_of_muscle_fibers) - ] + nmj_z, self.nerve_paths, branch_xy, branch_z, _ = two_layer_nmj( + self.muscle_fiber_centers__mm, + n_branches, + endplate_center, + branches_z_std, + arborization_z_std, + rng=get_random_generator(), + kmeans_seed=derive_subseed(n_branches), + return_branches=True, ) - + self.branch_points_xy = branch_xy + self.branch_points_z = branch_z + self._neuromuscular_z_coordinates__mm = nmj_z[:, None] self._actual_center = np.concatenate( - [ - np.mean(self.muscle_fiber_centers__mm, axis=0), - np.mean(self._neuromuscular_z_coordinates__mm, axis=0), - ] + [np.mean(self.muscle_fiber_centers__mm, axis=0), [np.mean(nmj_z)]] )[None] - for i in range(self._number_of_muscle_fibers): - cluster_center = np.concatenate( - [ - self.branch_points_xy[idx[i]], - self.branch_points_z[idx[i]], - ] - )[None] - nmj_coordinates = np.concatenate( - [ - self.muscle_fiber_centers__mm[i], - self._neuromuscular_z_coordinates__mm[i], - ] - )[None] - - # Extract scalar from the (1,)-shape norm result via .item() so the - # assignment to a 0-D slot doesn't trigger the numpy 1.25 ndim>0 - # deprecation (float() of a (1,)-shape array hits the same path). - self.nerve_paths[i, 0] = np.linalg.norm( - self._actual_center - cluster_center, axis=-1 - ).item() - self.nerve_paths[i, 1] = np.linalg.norm( - nmj_coordinates - cluster_center, axis=-1 - ).item() # Calculate delays # self._calculate_nmj_delays() diff --git a/myogen/simulator/core/emg/nerve_branches.py b/myogen/simulator/core/emg/nerve_branches.py new file mode 100644 index 00000000..e419269e --- /dev/null +++ b/myogen/simulator/core/emg/nerve_branches.py @@ -0,0 +1,73 @@ +"""Two-layer motor-nerve branching model (Konstantin et al. 2020). + +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. +Used by intramuscular EMG (``MotorUnitSim``) and surface EMG. + +Each motor unit's fibres are grouped into ``n_branches`` primary branches +(k-means on the fibre cross-section positions). A branch point sits at +``endplate_center + N(0, branches_z_std)``; each fibre's end-plate sits at its +branch point's z plus ``N(0, arborization_z_std)``. The nerve path to a fibre +is the distance from the motor unit's centre to its branch point (first layer) +plus the distance from the branch point to the end-plate (second layer). +""" +from __future__ import annotations + +import numpy as np +from sklearn.cluster import KMeans + + +def nmj_branch_schedule(recruitment_thresholds: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Branch count and z spreads per motor unit, growing with motor-unit size. + + Returns ``(n_branches, branches_z_std__mm, arborization_z_std__mm)``. The + spread factor is the share of total threshold below each unit. + """ + thresholds = np.asarray(recruitment_thresholds, dtype=float) + n_branches = (1 + np.round(np.log(thresholds / thresholds[0]))).astype(int) + spread = np.concatenate([[0.0], np.cumsum(thresholds)[:-1]]) / np.sum(thresholds) + return n_branches, 1.5 + spread * 4.0, 0.5 + spread * 1.5 + + +def two_layer_nmj( + fiber_xy__mm: np.ndarray, + n_branches: int, + endplate_center__mm: float, + branches_z_std__mm: float, + arborization_z_std__mm: float, + *, + rng: np.random.Generator, + kmeans_seed: int, + return_branches: bool = False, +): + """End-plate z per fibre and the two nerve-path lengths (mm) per fibre. + + Returns ``(nmj_z, nerve_paths)`` with ``nerve_paths[:, 0]`` the motor-unit + centre to branch point distance and ``nerve_paths[:, 1]`` the branch point + to end-plate distance. With ``return_branches`` the branch points + ``(xy, z)`` and each fibre's branch index are appended. + """ + fiber_xy__mm = np.asarray(fiber_xy__mm, dtype=float) + n_fibers = len(fiber_xy__mm) + kmeans = KMeans(n_clusters=n_branches, init="k-means++", max_iter=100, random_state=kmeans_seed) + idx = kmeans.fit_predict(fiber_xy__mm) + branch_xy = kmeans.cluster_centers_ + branch_z = endplate_center__mm + branches_z_std__mm * rng.standard_normal(size=(n_branches, 1)) + # One scalar draw per fibre, in fibre order (the published draw sequence). + nmj_z = np.array([branch_z[idx[i], 0] + arborization_z_std__mm * rng.standard_normal() + for i in range(n_fibers)]) + + centre = np.concatenate([np.mean(fiber_xy__mm, axis=0), [np.mean(nmj_z)]]) + branch = np.column_stack([branch_xy[idx], branch_z[idx, 0]]) + nmj = np.column_stack([fiber_xy__mm, nmj_z]) + paths = np.column_stack([np.linalg.norm(centre - branch, axis=1), + np.linalg.norm(nmj - branch, axis=1)]) + if return_branches: + return nmj_z, paths, branch_xy, branch_z, idx + return nmj_z, paths + + +def nerve_delays__ms(nerve_paths__mm: np.ndarray, branch_cvs__m_per_s) -> np.ndarray: + """Axon-branch conduction time per fibre (ms): sum of path / velocity per layer.""" + cvs = np.asarray(branch_cvs__m_per_s, dtype=float) # m/s == mm/ms + return np.sum(np.asarray(nerve_paths__mm, dtype=float) / cvs[None, :], axis=1) From 6a77faa935787d6590d2b842ea0c52d6b23ab12c Mon Sep 17 00:00:00 2001 From: Devon Rohlf Date: Sat, 31 Oct 2026 01:08:01 +0100 Subject: [PATCH 08/42] feat(semg): let a fibre propagate at its own velocity on a shared time grid v_grid fixes the (k_z, k_t) grid and window; the fibre velocity enters only the source term, whose IAP spectrum is read at k_t / v. Fibres at different velocities then sum on one time axis. delay_fibre_signal applies a fractional-sample delay. v_grid defaults to v, which leaves every existing call unchanged. --- .../core/emg/surface/simulate_fiber.py | 60 +++++++++++++++++-- 1 file changed, 54 insertions(+), 6 deletions(-) diff --git a/myogen/simulator/core/emg/surface/simulate_fiber.py b/myogen/simulator/core/emg/surface/simulate_fiber.py index 5ccec8e9..9293cafb 100644 --- a/myogen/simulator/core/emg/surface/simulate_fiber.py +++ b/myogen/simulator/core/emg/surface/simulate_fiber.py @@ -118,6 +118,37 @@ def fiber_semi_lengths(innervation_zone__mm: float, fiber_length__mm: float) -> return L1, L2 +def _source_spectrum(psi: np.ndarray, dz: float, stretch: float) -> np.ndarray: + """IAP spectrum for a fibre at ``v = v_grid / stretch`` on the ``v_grid`` time grid. + + ``psi`` is sampled on the symmetric grid ``z_n = (n - (N-1)/2) dz``. Column + ``m`` (wavenumber ``k_m`` of the centred FFT grid) needs the spectrum at + ``k_m * stretch``; the phase stays referred to the grid's first sample at + ``k_m`` so that, as on the fibre's own grid, the firing instant is the + window centre. ``stretch == 1`` gives ``fftshift(fft(psi)) * dz``. + """ + n = len(psi) + k = 2 * math.pi * np.fft.fftshift(np.fft.fftfreq(n, d=dz)) + support = np.nonzero(psi)[0] + z = (support - (n - 1) / 2) * dz + spectrum = (np.exp(-1j * np.outer(k * stretch, z)) @ psi[support]) * dz + return spectrum * np.exp(-1j * k * (n - 1) / 2 * dz) + + +def delay_fibre_signal(phi: np.ndarray, delay__ms: float, fs__kHz: float) -> np.ndarray: + """Delay ``phi`` along its last (time) axis by ``delay__ms``, allowing fractions of a sample. + + Applied as a linear phase in the frequency domain, i.e. a circular shift: + the window must leave room for the delayed signal. Zero delay returns + ``phi`` unchanged. + """ + if delay__ms == 0: + return phi + n = phi.shape[-1] + f__kHz = np.fft.rfftfreq(n, d=1.0 / fs__kHz) + return np.fft.irfft(np.fft.rfft(phi, axis=-1) * np.exp(-2j * np.pi * f__kHz * delay__ms), n=n, axis=-1) + + def muap_sampling_rate__Hz( conduction_velocity__m_s: float, internal_points: int, output_points: int, kernel_length__mm: float ) -> float: @@ -269,6 +300,7 @@ def _simulate_fiber_v2_python( *, sig_intracellular: float = 1.01, fiber_diameter__mm: float = 0.05, + v_grid: float | None = None, ): """ Simulate a single fiber (Python implementation). @@ -339,6 +371,11 @@ def _simulate_fiber_v2_python( Fibre diameter (mm), default 0.05. Together with ``sig_intracellular`` it sets the line current density ``sig_i * pi * (d/2)**2 * d2Vm/dz2`` of the source, so the output is the potential in mV. + v_grid : float, optional + Reference conduction velocity (m/s) that sets the time grid (``v`` above + in the window and sampling-rate formulas). Defaults to ``v``. Fibres + that propagate at different ``v`` but share ``v_grid`` return signals on + the same time axis and can be summed directly. Returns ------- @@ -349,6 +386,10 @@ def _simulate_fiber_v2_python( B_incomplete : np.ndarray B matrix for reuse in subsequent calls """ + # The grid velocity fixes the (k_z, k_t) grids; the fibre velocity enters + # only the source term below. + if v_grid is None: + v_grid = v # Get electrode configuration from the array channels = [electrode_array.num_rows, electrode_array.num_cols] @@ -367,7 +408,7 @@ def _simulate_fiber_v2_python( ## 1. Constants if _prepared is not None: - _prepared.bind((Fs, v, N, M, r, r_bone, th_fat, th_skin, + _prepared.bind((Fs, v_grid, N, M, r, r_bone, th_fat, th_skin, sig_muscle_rho, sig_muscle_z, sig_fat, sig_skin, sig_bone, fiber_length__mm, sig_intracellular, fiber_diameter__mm)) if _prepared is None or _prepared.grid is None: @@ -381,7 +422,7 @@ def _simulate_fiber_v2_python( if fiber_length__mm is not None: # z-grid spans fiber_length__mm with N points: dz = fiber_length__mm / (N-1) # For FFT: dz = v / Fs_eff → Fs_eff = v * (N-1) / fiber_length__mm - Fs_effective = v * (N - 1) / fiber_length__mm # kHz + Fs_effective = v_grid * (N - 1) / fiber_length__mm # kHz else: Fs_effective = Fs @@ -399,7 +440,7 @@ def _simulate_fiber_v2_python( if N > 1: dk_t = 2 * math.pi * Fs_effective / N k_t = np.where(k_t == 0.0, 0.5 * dk_t, k_t) - k_z = k_t / v + k_z = k_t / v_grid (kt_mesh_kzkt, kz_mesh_kzkt) = np.meshgrid(k_t, k_z) ## Model radii -> (Farina, 2004), Figure 1 - b) @@ -421,7 +462,7 @@ def _simulate_fiber_v2_python( # Define spatial grid for IAP kernel evaluation if fiber_length__mm is None: # Legacy behavior: z-range depends on N (INCORRECT - causes MUAP duration to scale with N) - z = np.linspace(-(N - 1) * (v / Fs) / 2, (N - 1) * (v / Fs) / 2, N, dtype=np.longdouble) + z = np.linspace(-(N - 1) * (v_grid / Fs) / 2, (N - 1) * (v_grid / Fs) / 2, N, dtype=np.longdouble) else: # Correct behavior: z-range defined by physical fiber length (independent of N) # N controls only the resolution, not the physical extent @@ -452,11 +493,18 @@ def _simulate_fiber_v2_python( _prepared.freeze_grid((k_theta, k_t, k_z, kt_mesh_kzkt, kz_mesh_kzkt, a, b, c, d, psi, PSI_mesh_conj, Fs_effective)) else: k_theta, k_t, k_z, kt_mesh_kzkt, kz_mesh_kzkt, a, b, c, d, psi, PSI_mesh_conj, Fs_effective = _prepared.grid + # The source term reads the IAP spectrum at k = k_t / v (one column per k_t). + # On the fibre's own grid that is the k_z grid; a fibre propagating at + # v != v_grid needs the spectrum at k_z * v_grid / v instead. + PSI_row_conj = PSI_mesh_conj[0] + if v != v_grid: + dz = (fiber_length__mm / (N - 1)) if fiber_length__mm is not None else v_grid / Fs + PSI_row_conj = np.conj(_source_spectrum(np.asarray(psi, dtype=np.float64), dz, v_grid / v)) if _prepared is not None: - I_kzkt = _prepared.source_current(k_z, k_t, PSI_mesh_conj[0], v, zi, L1, L2) + I_kzkt = _prepared.source_current(k_z, k_t, PSI_row_conj, v, zi, L1, L2) else: I_kzkt = 1j * np.multiply( - kz_mesh_kzkt / v, np.multiply(PSI_mesh_conj, np.exp(-1j * kz_mesh_kzkt * zi)) + kz_mesh_kzkt / v, np.multiply(PSI_row_conj[None, :], np.exp(-1j * kz_mesh_kzkt * zi)) ) k_eps = kz_mesh_kzkt + kt_mesh_kzkt / v k_beta = kz_mesh_kzkt - kt_mesh_kzkt / v From b72afc7a6d017672d43e247cde3ea92f6287fdac Mon Sep 17 00:00:00 2001 From: Devon Rohlf Date: Sat, 31 Oct 2026 01:08:02 +0100 Subject: [PATCH 09/42] feat(semg): optional per-fibre end-plates, nerve delays and fibre velocities SurfaceFiberDispersion switches each effect on independently: two-layer end-plate scatter (full or within each motor unit), axon-branch delays, size-ordered and within-unit fibre velocity spread, and diameter-scaled source strength. Draws use their own seeded generators; the default (off) is bit-identical to the undispersed model. The realised spread is reported per motor unit. The execution route and use_unified refuse an active configuration. --- myogen/simulator/core/emg/runtime/semg.py | 3 + .../core/emg/surface/fiber_dispersion.py | 183 +++++++++++++ .../simulator/core/emg/surface/surface_emg.py | 73 +++++- tests/test_surface_dispersion.py | 245 ++++++++++++++++++ 4 files changed, 501 insertions(+), 3 deletions(-) create mode 100644 myogen/simulator/core/emg/surface/fiber_dispersion.py create mode 100644 tests/test_surface_dispersion.py diff --git a/myogen/simulator/core/emg/runtime/semg.py b/myogen/simulator/core/emg/runtime/semg.py index e0a7aeb7..4e0dba2e 100644 --- a/myogen/simulator/core/emg/runtime/semg.py +++ b/myogen/simulator/core/emg/runtime/semg.py @@ -28,6 +28,9 @@ def simulate_muaps(owner, policy, n_jobs, verbose, use_gpu): raise ValueError("execution=ExecutionPolicy(backend='cpu') contradicts use_gpu=True") if owner._use_unified: raise ValueError("use_unified=True requires the existing route (execution=None)") + dispersion = getattr(owner, "_fiber_dispersion", None) + if dispersion is not None and dispersion.active: + raise ValueError("fiber_dispersion requires the existing route (execution=None)") if any(type(array) is not SurfaceElectrodeArray for array in owner._electrode_arrays): raise ValueError("Custom electrodes require the existing route (execution=None)") shapes = {(a.num_rows, a.num_cols) for a in owner._electrode_arrays} diff --git a/myogen/simulator/core/emg/surface/fiber_dispersion.py b/myogen/simulator/core/emg/surface/fiber_dispersion.py new file mode 100644 index 00000000..4646953c --- /dev/null +++ b/myogen/simulator/core/emg/surface/fiber_dispersion.py @@ -0,0 +1,183 @@ +"""Per-fibre dispersion for surface EMG: end-plates, axon-branch delays, fibre CV. + +Every option defaults to off, which reproduces the undispersed model (one +end-plate per motor unit, no nerve delays, every fibre at the mean conduction +velocity). The options are independent so that each effect can be measured +on its own. + +End-plates and nerve delays come from MyoGen's two-layer nerve model +(:mod:`myogen.simulator.core.emg.nerve_branches`), the same one intramuscular +EMG uses. Draws use generators seeded from ``seed`` only, never the global +MyoGen generator, so switching an option on leaves every other draw intact. +""" +from __future__ import annotations + +from dataclasses import asdict, dataclass, field +from typing import Literal, Optional, Sequence + +import numpy as np + +from ..nerve_branches import nerve_delays__ms, nmj_branch_schedule, two_layer_nmj + +NMJ_MODELS = ("uniform_mu", "two_layer", "two_layer_intra") + + +@dataclass(frozen=True) +class SurfaceFiberDispersion: + """Configuration of per-fibre dispersion (all off by default). + + Parameters + ---------- + nmj_model : {"uniform_mu", "two_layer", "two_layer_intra"} + ``"uniform_mu"``: every fibre of a motor unit has the unit's end-plate + (the undispersed model). ``"two_layer"``: the two-layer nerve model as + intramuscular EMG uses it (branch points scattered about the end-plate + band centre, fibres scattered about their branch point). + ``"two_layer_intra"``: the two-layer scatter re-centred on each unit's + undispersed end-plate, isolating the spread within a motor unit. + nmj_scatter_scale : float + Multiplies both two-layer z spreads (1 = the intramuscular default). + nerve_delays : bool + Delay each fibre by its axon-branch conduction time (two-layer nerve + paths over ``nerve_cvs__m_per_s``). Independent of ``nmj_model``: with + ``"uniform_mu"`` only the timing changes, not the end-plate positions. + nerve_cvs__m_per_s : (float, float) + Conduction velocity of the first (axon to branch point) and second + (branch point to end-plate) layer. Default (5, 2) m/s as in + intramuscular EMG (Konstantin et al. 2020). + fiber_cv_inter_mu_spread : float + Size-principle spread of motor-unit mean fibre CV: the first-recruited + unit propagates at ``mean * (1 - spread / 2)``, the last at + ``mean * (1 + spread / 2)``, linear in recruitment order. + fiber_cv_within_mu_cov : float + Coefficient of variation of fibre CV within a motor unit (normal, + clipped at three standard deviations). Stålberg & Karlsson (1973) + report 5-10 % in healthy muscle. + diameter_scales_source : bool + Scale each fibre's source by its cross-section, with diameter linear in + CV: ``d / d_ref = 1 + (cv - mean_cv) / (cv_per_diameter * d_ref)``. + cv_per_diameter__m_s_per_mm, reference_diameter__mm : float + Slope of CV against fibre diameter and the diameter at the mean CV. + seed : int + Seed of the dispersion draws. + """ + + nmj_model: Literal["uniform_mu", "two_layer", "two_layer_intra"] = "uniform_mu" + nmj_scatter_scale: float = 1.0 + nerve_delays: bool = False + nerve_cvs__m_per_s: tuple[float, float] = (5.0, 2.0) + fiber_cv_inter_mu_spread: float = 0.0 + fiber_cv_within_mu_cov: float = 0.0 + diameter_scales_source: bool = False + cv_per_diameter__m_s_per_mm: float = 77.0 + reference_diameter__mm: float = 0.05 + seed: int = 0 + + def __post_init__(self): + if self.nmj_model not in NMJ_MODELS: + raise ValueError(f"nmj_model must be one of {NMJ_MODELS}, got {self.nmj_model!r}") + if self.nmj_scatter_scale < 0: + raise ValueError("nmj_scatter_scale must be non-negative") + if len(self.nerve_cvs__m_per_s) != 2 or min(self.nerve_cvs__m_per_s) <= 0: + raise ValueError("nerve_cvs__m_per_s must be two positive velocities") + if not 0 <= self.fiber_cv_inter_mu_spread < 2: + raise ValueError("fiber_cv_inter_mu_spread must be in [0, 2)") + if not 0 <= self.fiber_cv_within_mu_cov < 1 / 3: + raise ValueError("fiber_cv_within_mu_cov must be in [0, 1/3)") + + @property + def active(self) -> bool: + return (self.nmj_model != "uniform_mu" or self.nerve_delays + or self.fiber_cv_inter_mu_spread > 0 or self.fiber_cv_within_mu_cov > 0 + or self.diameter_scales_source) + + +@dataclass +class FiberDispersionRealization: + """Per-motor-unit lists of per-fibre values, and a per-unit summary.""" + + zi__mm: list[np.ndarray] + cv__m_s: list[np.ndarray] + delay__ms: list[np.ndarray] + amplitude_scale: list[np.ndarray] + report: dict = field(default_factory=dict) + + +def realize_fiber_dispersion( + config: SurfaceFiberDispersion, + fiber_positions__mm: Sequence[np.ndarray], + recruitment_thresholds: np.ndarray, + mu_zones__mm: np.ndarray, + mean_cv__m_s: float, + fiber_lengths__mm: Optional[Sequence[np.ndarray]] = None, + end_margin__mm: float = 1.0, +) -> FiberDispersionRealization: + """Draw per-fibre end-plate, nerve delay, CV and source scale for every motor unit. + + ``mu_zones__mm`` holds each unit's undispersed end-plate (relative to the + fibre centre). With ``fiber_lengths__mm`` given, end-plates are clipped to + ``end_margin__mm`` inside each fibre and the clipped count is reported. + """ + n_mu = len(fiber_positions__mm) + counts = np.array([len(p) for p in fiber_positions__mm]) + nerve_rng = np.random.default_rng([config.seed, 0]) + cv_rng = np.random.default_rng([config.seed, 1]) + needs_nerve = config.nmj_model != "uniform_mu" or config.nerve_delays + if needs_nerve: + n_branches, branch_sd, arbor_sd = nmj_branch_schedule(recruitment_thresholds) + + report = {key: np.full(n_mu, np.nan) for key in ( + "zi_mean__mm", "zi_sd__mm", "cv_mean__m_s", "cv_cov", "delay_mean__ms", "delay_sd__ms", + "amplitude_scale_mean")} + report["n_fibers"] = counts + report["n_branches"] = np.zeros(n_mu, dtype=int) + report["zi_clipped"] = np.zeros(n_mu, dtype=int) + zi_all, cv_all, delay_all, scale_all = [], [], [], [] + rank = np.arange(n_mu) / max(n_mu - 1, 1) + for mu in range(n_mu): + n = int(counts[mu]) + zi = np.full(n, float(mu_zones__mm[mu])) + delay = np.zeros(n) + if needs_nerve and n: + branches = int(min(n_branches[mu], n)) + report["n_branches"][mu] = branches + nmj_z, paths = two_layer_nmj( + fiber_positions__mm[mu], branches, 0.0, + config.nmj_scatter_scale * branch_sd[mu], config.nmj_scatter_scale * arbor_sd[mu], + rng=nerve_rng, kmeans_seed=int(nerve_rng.integers(2**31))) + if config.nmj_model == "two_layer": + zi = nmj_z + elif config.nmj_model == "two_layer_intra": + zi = zi + (nmj_z - np.mean(nmj_z)) + if config.nerve_delays: + delay = nerve_delays__ms(paths, config.nerve_cvs__m_per_s) + if fiber_lengths__mm is not None and n: + limit = np.asarray(fiber_lengths__mm[mu]) / 2 - end_margin__mm + clipped = np.clip(zi, -limit, limit) + report["zi_clipped"][mu] = int(np.count_nonzero(clipped != zi)) + zi = clipped + + cv_mu = mean_cv__m_s * (1 - config.fiber_cv_inter_mu_spread / 2 + + config.fiber_cv_inter_mu_spread * rank[mu]) + cov = config.fiber_cv_within_mu_cov + if cov > 0: + cv = np.clip(cv_rng.normal(cv_mu, cov * cv_mu, size=n), cv_mu * (1 - 3 * cov), cv_mu * (1 + 3 * cov)) + else: + cv = np.full(n, cv_mu) + if config.diameter_scales_source: + ratio = 1 + (cv - mean_cv__m_s) / (config.cv_per_diameter__m_s_per_mm * config.reference_diameter__mm) + scale = ratio**2 + else: + scale = np.ones(n) + + zi_all.append(zi) + cv_all.append(cv) + delay_all.append(delay) + scale_all.append(scale) + if n: + report["zi_mean__mm"][mu], report["zi_sd__mm"][mu] = np.mean(zi), np.std(zi) + report["cv_mean__m_s"][mu], report["cv_cov"][mu] = np.mean(cv), np.std(cv) / np.mean(cv) + report["delay_mean__ms"][mu], report["delay_sd__ms"][mu] = np.mean(delay), np.std(delay) + report["amplitude_scale_mean"][mu] = np.mean(scale) + report["config"] = asdict(config) + return FiberDispersionRealization(zi_all, cv_all, delay_all, scale_all, report) diff --git a/myogen/simulator/core/emg/surface/surface_emg.py b/myogen/simulator/core/emg/surface/surface_emg.py index af036f1f..e36c532f 100644 --- a/myogen/simulator/core/emg/surface/surface_emg.py +++ b/myogen/simulator/core/emg/surface/surface_emg.py @@ -24,10 +24,15 @@ from myogen.simulator.core.emg.electrodes import SurfaceElectrodeArray from myogen.simulator.core.emg.surface.simulate_fiber import ( _simulate_fiber_v2_python, + delay_fibre_signal, fiber_semi_lengths, muap_sampling_rate__Hz, simulate_fiber_v2, ) +from myogen.simulator.core.emg.surface.fiber_dispersion import ( + SurfaceFiberDispersion, + realize_fiber_dispersion, +) from myogen.simulator.core.muscle import Muscle from myogen.utils.decorators import beartowertype from myogen.utils.types import ( @@ -90,6 +95,12 @@ class SurfaceEMG: intracellular_conductivity__S_m : float, optional Intracellular conductivity in S/m, default 1.01. The source current density is ``sigma_i * pi * (d/2)**2 * d2Vm/dz2`` (core-conductor model). + fiber_dispersion : SurfaceFiberDispersion, optional + Per-fibre end-plate scatter, axon-branch delays and fibre conduction + velocities (see :class:`SurfaceFiberDispersion`). None (default) or an + inactive configuration keeps one end-plate per motor unit, no nerve + delays and the mean velocity for every fibre. Supported on the default + route only (``execution=None``, ``use_unified=False``). Attributes ---------- @@ -101,6 +112,10 @@ class SurfaceEMG: Noisy surface EMG signals for each electrode array as a neo.Block. Available after add_noise(). spike_train__Block : SPIKE_TRAIN__Block Spike train block used for EMG generation signals. Available after simulate_surface_emg(). + fiber_dispersion_report : dict or None + Per-motor-unit summary of the realised dispersion (end-plate mean and SD, + CV mean and coefficient of variation, nerve delay mean and SD, source + scale), or None when no dispersion is active. Available after simulate_muaps(). References ---------- @@ -120,6 +135,7 @@ def __init__( use_unified: bool = False, fiber_diameter__um: float = 50.0, intracellular_conductivity__S_m: float = 1.01, + fiber_dispersion: SurfaceFiberDispersion | None = None, ): # Immutable public arguments - never modify these self.muscle_model = muscle_model @@ -132,6 +148,9 @@ def __init__( self._use_unified = use_unified self.fiber_diameter__um = fiber_diameter__um self.intracellular_conductivity__S_m = intracellular_conductivity__S_m + self.fiber_dispersion = fiber_dispersion + self._fiber_dispersion = fiber_dispersion + self.fiber_dispersion_report = None # Internal sampling frequency for higher resolution MUAP computation # If not specified, defaults to 10 kHz for better MUAP resolution @@ -318,6 +337,28 @@ def simulate_muaps( for mu_idx in range(n_motor_units) ] + # Per-fibre dispersion, drawn here from its own seeded generators (the + # global draws above are unchanged whether or not it is active). + dispersion = None + self.fiber_dispersion_report = None + if self._fiber_dispersion is not None and self._fiber_dispersion.active: + if self._use_unified: + raise ValueError("fiber_dispersion is not supported with use_unified=True") + fiber_positions = [] + for mu_idx in range(n_motor_units): + positions = self._muscle_model.resulting_fiber_assignment(mu_idx) + fiber_positions.append(np.asarray(getattr(positions, "magnitude", positions))[:, :2]) + dispersion = realize_fiber_dispersion( + self._fiber_dispersion, + fiber_positions, + np.asarray(self._muscle_model.recruitment_thresholds), + innervation_zones, + mean_cv__m_s=float(self._mean_conduction_velocity__m_s), + fiber_lengths__mm=[self._mean_fiber_length__mm + variation + for variation in fiber_length_variations_per_MU], + ) + self.fiber_dispersion_report = dispersion.report + # Result shapes follow each electrode array's own grid (arrays may differ). def internal_result_shape(electrode_array: SurfaceElectrodeArray) -> tuple[int, int, int]: return electrode_array.num_rows, electrode_array.num_cols, len(t_internal) @@ -426,6 +467,20 @@ def _process_single_mu( else: IAP_SCALE_FACTOR = 2.5 kernel_length = self._mean_fiber_length__mm * IAP_SCALE_FACTOR + Fs_grid = v_conduction * (N_internal - 1) / kernel_length # kHz + + if dispersion is not None: + # The fibre must finish inside the second half of the window + # (firing at its centre), or the circular grid wraps it. + zone = dispersion.zi__mm[MU_index] + reach = fiber_lengths / 2 + np.abs(zone) + 10.0 # mm, plus the IAP tail + finish = dispersion.delay__ms[MU_index] + reach / dispersion.cv__m_s[MU_index] + if np.max(finish) >= kernel_length / v_conduction / 2: + raise ValueError( + f"MU {MU_index}: dispersed fibres need {np.max(finish):.1f} ms after firing, " + f"the window holds {kernel_length / v_conduction / 2:.1f} ms; " + "increase iap_kernel_length__mm" + ) # Process each fiber (inner loop - must remain sequential) fiber_iter = tqdm( @@ -439,9 +494,14 @@ def _process_single_mu( R = R_values[fiber_number] theta = theta_values[fiber_number] fiber_length__mm = fiber_lengths[fiber_number] + if dispersion is not None: + fiber_zone = float(dispersion.zi__mm[MU_index][fiber_number]) + fiber_cv = float(dispersion.cv__m_s[MU_index][fiber_number]) + else: + fiber_zone, fiber_cv = innervation_zone, v_conduction # Semi-lengths towards +z (L1) and -z (L2) from the end-plate - L1, L2 = fiber_semi_lengths(innervation_zone, fiber_length__mm) + L1, L2 = fiber_semi_lengths(fiber_zone, fiber_length__mm) if self._use_unified: from myogen.simulator.core.emg.fiber_simulation import ( @@ -485,7 +545,8 @@ def _process_single_mu( # Existing frequency-domain path (unchanged) phi_temp, A_matrix, B_incomplete = _simulate_fiber_v2_python( Fs=Fs_internal, - v=v_conduction, + v=fiber_cv, + v_grid=v_conduction, N=N_internal, M=M_theta, r=r_total, @@ -495,7 +556,7 @@ def _process_single_mu( R=R, L1=L1, L2=L2, - zi=innervation_zone, + zi=fiber_zone, electrode_array=electrode_array, sig_muscle_rho=sig_rho, sig_muscle_z=sig_z, @@ -516,6 +577,12 @@ def _process_single_mu( fiber_diameter__mm=self._fiber_diameter__mm, ) + if dispersion is not None: + phi_temp = phi_temp * dispersion.amplitude_scale[MU_index][fiber_number] + phi_temp = delay_fibre_signal( + phi_temp, float(dispersion.delay__ms[MU_index][fiber_number]), Fs_grid + ) + array_result_internal += phi_temp # Downsample from internal resolution to output resolution diff --git a/tests/test_surface_dispersion.py b/tests/test_surface_dispersion.py new file mode 100644 index 00000000..f7d66a09 --- /dev/null +++ b/tests/test_surface_dispersion.py @@ -0,0 +1,245 @@ +"""Per-fibre dispersion in the surface-EMG model. + +Axon-branch delays and end-plate positions come from the two-layer nerve model +that basic MyoGen already uses for intramuscular EMG (Konstantin et al. 2020); +muscle-fibre conduction velocity may vary per fibre. Every option defaults to +off, and off must reproduce the undispersed model exactly. +""" +import numpy as np +import pytest +import quantities as pq +from scipy.signal import resample + +from myogen import derive_subseed, get_random_generator, set_random_seed +from myogen.simulator import ( + ExecutionPolicy, + Muscle, + RecruitmentThresholds, + SurfaceElectrodeArray, + SurfaceEMG, +) +from myogen.simulator.core.emg.intramuscular.motor_unit_sim import MotorUnitSim +from myogen.simulator.core.emg.nerve_branches import ( + nerve_delays__ms, + nmj_branch_schedule, + two_layer_nmj, +) +from myogen.simulator.core.emg.surface import simulate_fiber +from myogen.simulator.core.emg.surface.fiber_dispersion import ( + SurfaceFiberDispersion, + realize_fiber_dispersion, +) +from myogen.utils.neo import signal_to_grid + +GEOMETRY = dict(r=8.5, r_bone=0.0, th_fat=0.3, th_skin=1.29, sig_muscle_rho=0.09, + sig_muscle_z=0.4, sig_fat=4.07e-2, sig_skin=1.0) +N = 512 +KERNEL = 100.0 # mm + + +def _array(z0=10.0, mode="monopolar"): + return SurfaceElectrodeArray(1, 1, 10 * pq.mm, 0.5 * pq.mm, + center_point__mm_deg=(z0 * pq.mm, 0 * pq.deg), + bending_radius__mm=GEOMETRY["r"] * pq.mm, + differentiation_mode=mode) + + +def _fibre(v, v_grid=None, n=N, zi=0.0): + kwargs = {} if v_grid is None else {"v_grid": v_grid} + phi, _, _ = simulate_fiber._simulate_fiber_v2_python( + Fs=10.0, v=v, N=n, M=32, R=6.5, L1=20.0 - zi, L2=20.0 + zi, zi=zi, + electrode_array=_array(), fiber_length__mm=KERNEL, use_gpu=False, **kwargs, **GEOMETRY) + return phi[0, 0] + + +# ---- two-layer nerve model: one authority shared with intramuscular EMG ---------- + +def test_two_layer_function_reproduces_motor_unit_sim(): + set_random_seed(5) + xy = get_random_generator().uniform(-3, 3, size=(40, 2)) + set_random_seed(11) + sim = MotorUnitSim(muscle_fiber_centers__mm=xy, muscle_length__mm=80.0, + muscle_fiber_diameters__mm=np.full(40, 0.05), + muscle_fiber_conduction_velocity__mm_per_s=np.full(40, 4000.0), + nominal_center__mm=np.zeros(2)) + sim.sim_nmj_branches_two_layers(3, 40.0, 2.0, 0.8) + set_random_seed(11) + nmj_z, paths = two_layer_nmj(xy, 3, 40.0, 2.0, 0.8, rng=get_random_generator(), + kmeans_seed=derive_subseed(3)) + np.testing.assert_allclose(np.ravel(sim._neuromuscular_z_coordinates__mm), nmj_z, rtol=0, atol=1e-12) + np.testing.assert_allclose(sim.nerve_paths, paths, rtol=0, atol=1e-12) + + +def test_branch_schedule_grows_with_motor_unit_size(): + thresholds = np.array([1.0, np.e, np.e**2, np.e**3]) + n_branches, branch_sd, arbor_sd = nmj_branch_schedule(thresholds) + np.testing.assert_array_equal(n_branches, [1, 2, 3, 4]) + spread = np.concatenate([[0], np.cumsum(thresholds)[:-1]]) / thresholds.sum() + np.testing.assert_allclose(branch_sd, 1.5 + 4.0 * spread) + np.testing.assert_allclose(arbor_sd, 0.5 + 1.5 * spread) + + +def test_nerve_delay_is_path_over_branch_velocity(): + paths = np.array([[5.0, 2.0], [0.0, 4.0]]) # mm + np.testing.assert_allclose(nerve_delays__ms(paths, (5.0, 2.0)), [2.0, 2.0]) + + +# ---- fibre velocity is decoupled from the time grid ------------------------------- + +def test_reference_grid_velocity_defaults_to_the_fibre_velocity(): + np.testing.assert_array_equal(_fibre(4.2, v_grid=4.2), _fibre(4.2)) + + +@pytest.mark.parametrize("v", [3.4, 5.0]) +def test_fibre_on_a_shared_grid_matches_its_own_grid(v): + # Native run: dt = KERNEL / (v (N-1)). Shared run: dt = KERNEL / (4.2 (N-1)). + n = 1024 + native = _fibre(v, n=n) + shared = _fibre(v, v_grid=4.2, n=n) + t_native = (np.arange(n) - n // 2) * KERNEL / (v * (n - 1)) + t_shared = (np.arange(n) - n // 2) * KERNEL / (4.2 * (n - 1)) + up = 16 + fine = resample(native, n * up) + t_fine = (np.arange(n * up) - n // 2 * up) * KERNEL / (v * (n - 1)) / up + expected = np.interp(t_shared, t_fine, fine, left=0.0, right=0.0) + inside = np.abs(t_shared) < 0.45 * min(t_native.max(), t_shared.max()) + err = np.max(np.abs(shared[inside] - expected[inside])) / np.ptp(expected) + assert err < 0.01 + + +# ---- fractional delay --------------------------------------------------------------- + +def test_delay_moves_the_signal_by_the_requested_time(): + fs = 20.0 # kHz + t = np.arange(2048) / fs + x = np.exp(-((t - 40.0) / 0.8) ** 2)[None, None, :] + y = simulate_fiber.delay_fibre_signal(x, 1.37, fs) + up = resample(y[0, 0], 2048 * 50) + assert np.argmax(up) / 50 / fs == pytest.approx(41.37, abs=0.01) + np.testing.assert_array_equal(simulate_fiber.delay_fibre_signal(x, 0.0, fs), x) + + +# ---- realisation -------------------------------------------------------------------- + +def _muscle(): + set_random_seed(42) + thresholds, _ = RecruitmentThresholds(N=12, recruitment_range__ratio=30, mode="fuglevand") + muscle = Muscle(thresholds, radius__mm=6 * pq.mm, length__mm=60 * pq.mm, + fiber_density__fibers_per_mm2=6 * pq.mm**-2, mean_fiber_length__mm=50 * pq.mm, + var_fiber_length__mm=2 * pq.mm, radius_bone__mm=0 * pq.mm, + fat_thickness__mm=1 * pq.mm) + muscle.distribute_innervation_centers() + muscle.generate_muscle_fiber_centers(verbose=False) + muscle.assign_mfs2mns(verbose=False, n_jobs=1) + return muscle + + +@pytest.fixture(scope="module") +def muscle(): + return _muscle() + + +def _realize(muscle, config, zones=None): + counts = np.asarray(muscle.resulting_number_of_innervated_fibers) + zones = np.zeros(len(counts)) if zones is None else zones + positions = [np.asarray(muscle.resulting_fiber_assignment(i)) for i in range(len(counts))] + return realize_fiber_dispersion(config, positions, np.asarray(muscle.recruitment_thresholds), + zones, mean_cv__m_s=4.2) + + +def test_default_config_is_the_undispersed_model(muscle): + zones = np.linspace(-2, 2, 12) + real = _realize(muscle, SurfaceFiberDispersion(), zones) + for mu, z in enumerate(zones): + n = len(real.cv__m_s[mu]) + np.testing.assert_array_equal(real.zi__mm[mu], np.full(n, z)) + np.testing.assert_array_equal(real.cv__m_s[mu], np.full(n, 4.2)) + np.testing.assert_array_equal(real.delay__ms[mu], np.zeros(n)) + np.testing.assert_array_equal(real.amplitude_scale[mu], np.ones(n)) + assert not SurfaceFiberDispersion().active + + +def test_fibre_cv_spreads_follow_the_configuration(muscle): + real = _realize(muscle, SurfaceFiberDispersion(fiber_cv_inter_mu_spread=0.3, + fiber_cv_within_mu_cov=0.08)) + means = np.array([np.mean(c) for c in real.cv__m_s if len(c)]) + assert means[0] == pytest.approx(4.2 * 0.85, rel=0.05) + assert means[-1] == pytest.approx(4.2 * 1.15, rel=0.05) + big = max(range(12), key=lambda i: len(real.cv__m_s[i])) + c = real.cv__m_s[big] + assert np.std(c) / np.mean(c) == pytest.approx(0.08, rel=0.25) + assert np.all(real.amplitude_scale[big] == 1.0) # amplitude coupling is a separate switch + + +def test_diameter_coupling_scales_source_with_cross_section(muscle): + cfg = SurfaceFiberDispersion(fiber_cv_within_mu_cov=0.08, diameter_scales_source=True) + real = _realize(muscle, cfg) + big = max(range(12), key=lambda i: len(real.cv__m_s[i])) + d_ratio = 1 + (real.cv__m_s[big] - 4.2) / (cfg.cv_per_diameter__m_s_per_mm * cfg.reference_diameter__mm) + np.testing.assert_allclose(real.amplitude_scale[big], d_ratio**2) + + +def test_intra_mu_scatter_keeps_each_motor_unit_anchor(muscle): + zones = np.linspace(-2, 2, 12) + real = _realize(muscle, SurfaceFiberDispersion(nmj_model="two_layer_intra"), zones) + for mu in range(12): + if len(real.zi__mm[mu]) > 1: + assert np.mean(real.zi__mm[mu]) == pytest.approx(zones[mu], abs=1e-9) + assert np.std(real.zi__mm[mu]) > 0.3 + + +def test_nerve_delays_switch_independently_of_end_plates(muscle): + real = _realize(muscle, SurfaceFiberDispersion(nerve_delays=True)) + big = max(range(12), key=lambda i: len(real.cv__m_s[i])) + np.testing.assert_array_equal(real.zi__mm[big], np.zeros(len(real.zi__mm[big]))) + assert np.std(real.delay__ms[big]) > 0 + assert real.report["delay_sd__ms"][big] == pytest.approx(np.std(real.delay__ms[big])) + + +def test_realisation_leaves_the_global_generator_untouched(muscle): + set_random_seed(3) + expected = get_random_generator().uniform(size=4) + set_random_seed(3) + _realize(muscle, SurfaceFiberDispersion(nmj_model="two_layer", nerve_delays=True, + fiber_cv_within_mu_cov=0.05)) + np.testing.assert_array_equal(get_random_generator().uniform(size=4), expected) + + +# ---- SurfaceEMG integration --------------------------------------------------------- + +def _semg(muscle, dispersion=None, mus=(4, 11)): + array = SurfaceElectrodeArray(1, 1, 10 * pq.mm, 2.5 * pq.mm, + center_point__mm_deg=(12 * pq.mm, 0 * pq.deg), + bending_radius__mm=8.3 * pq.mm, + differentiation_mode="bipolar_longitudinal") + set_random_seed(7) + return SurfaceEMG(muscle, [array], sampling_frequency__Hz=2048 * pq.Hz, + internal_sampling_frequency__Hz=2048 * pq.Hz, + MUs_to_simulate=list(mus), fiber_dispersion=dispersion) + + +def _muaps(semg): + block = semg.simulate_muaps(n_jobs=1, verbose=False, use_gpu=False) + return np.array([np.asarray(signal_to_grid(s.analogsignals[0].rescale(pq.mV))) + for s in block.groups[0].segments]) + + +def test_inactive_dispersion_is_bit_identical(muscle): + np.testing.assert_array_equal(_muaps(_semg(muscle)), + _muaps(_semg(muscle, SurfaceFiberDispersion()))) + + +def test_dispersion_changes_the_muap_and_is_reported(muscle): + base = _muaps(_semg(muscle)) + semg = _semg(muscle, SurfaceFiberDispersion(nerve_delays=True, fiber_cv_within_mu_cov=0.08)) + dispersed = _muaps(semg) + assert not np.allclose(base[11], dispersed[11]) + report = semg.fiber_dispersion_report + assert report["delay_sd__ms"][11] > 0 and report["cv_cov"][11] > 0 + + +@pytest.mark.usefixtures("ample_session_memory") +def test_other_routes_refuse_active_dispersion(muscle): + semg = _semg(muscle, SurfaceFiberDispersion(nerve_delays=True)) + with pytest.raises(ValueError, match="dispersion"): + semg.simulate_muaps(n_jobs=1, verbose=False, execution=ExecutionPolicy()) From 4c80ca30eb61cb51e7beb3ab755558b2e29033e4 Mon Sep 17 00:00:00 2001 From: Devon Rohlf Date: Sat, 31 Oct 2026 11:51:15 +0100 Subject: [PATCH 10/42] feat(semg): fibre CV linear in recruitment threshold; literature CV-diameter line Add fiber_cv_threshold_slope (Del Vecchio et al. 2017: motor-unit CV rises linearly with recruitment threshold). Read fibre diameter from the Nandedkar & Stalberg (1983) line cv = 0.95 m/s + 0.05 m/s/um * d instead of an unsourced slope and reference diameter, and report the mean diameter. Drop an unverified citation for the within-unit CV spread. --- .../core/emg/surface/fiber_dispersion.py | 59 +++++++++++++------ tests/test_surface_dispersion.py | 16 ++++- 2 files changed, 56 insertions(+), 19 deletions(-) diff --git a/myogen/simulator/core/emg/surface/fiber_dispersion.py b/myogen/simulator/core/emg/surface/fiber_dispersion.py index 4646953c..043faa4e 100644 --- a/myogen/simulator/core/emg/surface/fiber_dispersion.py +++ b/myogen/simulator/core/emg/surface/fiber_dispersion.py @@ -46,18 +46,28 @@ class SurfaceFiberDispersion: (branch point to end-plate) layer. Default (5, 2) m/s as in intramuscular EMG (Konstantin et al. 2020). fiber_cv_inter_mu_spread : float - Size-principle spread of motor-unit mean fibre CV: the first-recruited - unit propagates at ``mean * (1 - spread / 2)``, the last at - ``mean * (1 + spread / 2)``, linear in recruitment order. + Spread of motor-unit mean fibre CV linear in recruitment order: the + first-recruited unit propagates at ``mean * (1 - spread / 2)``, the + last at ``mean * (1 + spread / 2)``. + fiber_cv_threshold_slope : float + Motor-unit mean fibre CV linear in recruitment threshold: + ``mean * (1 + slope * (rt - mean(rt)))`` with ``rt`` the threshold + relative to the largest. Del Vecchio et al. (2017) measured CV rising + linearly with recruitment threshold in tibialis anterior (406 units; + 3.89 m/s at 0-30 % MVC, 5.12 m/s at 50-70 % MVC, overall mean 4.32), + about 0.63 of the mean per unit of threshold. Cannot be combined with + ``fiber_cv_inter_mu_spread``. fiber_cv_within_mu_cov : float Coefficient of variation of fibre CV within a motor unit (normal, - clipped at three standard deviations). Stålberg & Karlsson (1973) - report 5-10 % in healthy muscle. + clipped at three standard deviations). No measurement of this spread + was found; Del Vecchio et al. (2017) call it unknown. diameter_scales_source : bool - Scale each fibre's source by its cross-section, with diameter linear in - CV: ``d / d_ref = 1 + (cv - mean_cv) / (cv_per_diameter * d_ref)``. - cv_per_diameter__m_s_per_mm, reference_diameter__mm : float - Slope of CV against fibre diameter and the diameter at the mean CV. + Scale each fibre's source by its cross-section, with the diameter read + from the CV line ``cv = cv_intercept + cv_per_diameter * d`` + (Nandedkar & Stålberg 1983: 0.95 m/s + 0.05 m/s per µm), relative to + the diameter at the mean CV. + cv_per_diameter__m_s_per_mm, cv_intercept__m_s : float + Slope and intercept of that line. seed : int Seed of the dispersion draws. """ @@ -67,10 +77,11 @@ class SurfaceFiberDispersion: nerve_delays: bool = False nerve_cvs__m_per_s: tuple[float, float] = (5.0, 2.0) fiber_cv_inter_mu_spread: float = 0.0 + fiber_cv_threshold_slope: float = 0.0 fiber_cv_within_mu_cov: float = 0.0 diameter_scales_source: bool = False - cv_per_diameter__m_s_per_mm: float = 77.0 - reference_diameter__mm: float = 0.05 + cv_per_diameter__m_s_per_mm: float = 50.0 + cv_intercept__m_s: float = 0.95 seed: int = 0 def __post_init__(self): @@ -82,13 +93,18 @@ def __post_init__(self): raise ValueError("nerve_cvs__m_per_s must be two positive velocities") if not 0 <= self.fiber_cv_inter_mu_spread < 2: raise ValueError("fiber_cv_inter_mu_spread must be in [0, 2)") + if self.fiber_cv_threshold_slope < 0: + raise ValueError("fiber_cv_threshold_slope must be non-negative") + if self.fiber_cv_threshold_slope and self.fiber_cv_inter_mu_spread: + raise ValueError("use either fiber_cv_threshold_slope or fiber_cv_inter_mu_spread, not both") if not 0 <= self.fiber_cv_within_mu_cov < 1 / 3: raise ValueError("fiber_cv_within_mu_cov must be in [0, 1/3)") @property def active(self) -> bool: return (self.nmj_model != "uniform_mu" or self.nerve_delays - or self.fiber_cv_inter_mu_spread > 0 or self.fiber_cv_within_mu_cov > 0 + or self.fiber_cv_inter_mu_spread > 0 or self.fiber_cv_threshold_slope > 0 + or self.fiber_cv_within_mu_cov > 0 or self.diameter_scales_source) @@ -128,12 +144,20 @@ def realize_fiber_dispersion( report = {key: np.full(n_mu, np.nan) for key in ( "zi_mean__mm", "zi_sd__mm", "cv_mean__m_s", "cv_cov", "delay_mean__ms", "delay_sd__ms", - "amplitude_scale_mean")} + "amplitude_scale_mean", "diameter_mean__um")} report["n_fibers"] = counts report["n_branches"] = np.zeros(n_mu, dtype=int) report["zi_clipped"] = np.zeros(n_mu, dtype=int) zi_all, cv_all, delay_all, scale_all = [], [], [], [] rank = np.arange(n_mu) / max(n_mu - 1, 1) + rt = np.asarray(recruitment_thresholds, dtype=float) + rt = rt / rt.max() + threshold_factor = 1 + config.fiber_cv_threshold_slope * (rt - rt.mean()) + if np.any(threshold_factor <= 0): + raise ValueError("fiber_cv_threshold_slope gives a non-positive CV for the first-recruited units") + d_mean__um = (mean_cv__m_s - config.cv_intercept__m_s) / config.cv_per_diameter__m_s_per_mm * 1e3 + if config.diameter_scales_source and d_mean__um <= 0: + raise ValueError("the mean CV lies below the CV-diameter intercept") for mu in range(n_mu): n = int(counts[mu]) zi = np.full(n, float(mu_zones__mm[mu])) @@ -157,16 +181,16 @@ def realize_fiber_dispersion( report["zi_clipped"][mu] = int(np.count_nonzero(clipped != zi)) zi = clipped - cv_mu = mean_cv__m_s * (1 - config.fiber_cv_inter_mu_spread / 2 - + config.fiber_cv_inter_mu_spread * rank[mu]) + cv_mu = mean_cv__m_s * threshold_factor[mu] * (1 - config.fiber_cv_inter_mu_spread / 2 + + config.fiber_cv_inter_mu_spread * rank[mu]) cov = config.fiber_cv_within_mu_cov if cov > 0: cv = np.clip(cv_rng.normal(cv_mu, cov * cv_mu, size=n), cv_mu * (1 - 3 * cov), cv_mu * (1 + 3 * cov)) else: cv = np.full(n, cv_mu) + diameter__um = (cv - config.cv_intercept__m_s) / config.cv_per_diameter__m_s_per_mm * 1e3 if config.diameter_scales_source: - ratio = 1 + (cv - mean_cv__m_s) / (config.cv_per_diameter__m_s_per_mm * config.reference_diameter__mm) - scale = ratio**2 + scale = (diameter__um / d_mean__um) ** 2 else: scale = np.ones(n) @@ -179,5 +203,6 @@ def realize_fiber_dispersion( report["cv_mean__m_s"][mu], report["cv_cov"][mu] = np.mean(cv), np.std(cv) / np.mean(cv) report["delay_mean__ms"][mu], report["delay_sd__ms"][mu] = np.mean(delay), np.std(delay) report["amplitude_scale_mean"][mu] = np.mean(scale) + report["diameter_mean__um"][mu] = np.mean(diameter__um) report["config"] = asdict(config) return FiberDispersionRealization(zi_all, cv_all, delay_all, scale_all, report) diff --git a/tests/test_surface_dispersion.py b/tests/test_surface_dispersion.py index f7d66a09..4c823993 100644 --- a/tests/test_surface_dispersion.py +++ b/tests/test_surface_dispersion.py @@ -171,12 +171,24 @@ def test_fibre_cv_spreads_follow_the_configuration(muscle): assert np.all(real.amplitude_scale[big] == 1.0) # amplitude coupling is a separate switch +def test_fibre_cv_rises_linearly_with_recruitment_threshold(muscle): + real = _realize(muscle, SurfaceFiberDispersion(fiber_cv_threshold_slope=0.63)) + rt = np.asarray(muscle.recruitment_thresholds, dtype=float) + rt = rt / rt.max() + for mu, cv in enumerate(real.cv__m_s): + np.testing.assert_allclose(cv, 4.2 * (1 + 0.63 * (rt[mu] - rt.mean()))) + with pytest.raises(ValueError): + SurfaceFiberDispersion(fiber_cv_threshold_slope=0.6, fiber_cv_inter_mu_spread=0.3) + + def test_diameter_coupling_scales_source_with_cross_section(muscle): cfg = SurfaceFiberDispersion(fiber_cv_within_mu_cov=0.08, diameter_scales_source=True) real = _realize(muscle, cfg) big = max(range(12), key=lambda i: len(real.cv__m_s[i])) - d_ratio = 1 + (real.cv__m_s[big] - 4.2) / (cfg.cv_per_diameter__m_s_per_mm * cfg.reference_diameter__mm) - np.testing.assert_allclose(real.amplitude_scale[big], d_ratio**2) + # Diameter from the CV line v = 0.95 m/s + 0.05 m/s/um * d (Nandedkar & Stalberg 1983) + diameter__um = (real.cv__m_s[big] - 0.95) / 0.05 + np.testing.assert_allclose(real.amplitude_scale[big], (diameter__um / ((4.2 - 0.95) / 0.05)) ** 2) + assert real.report["diameter_mean__um"][big] == pytest.approx(np.mean(diameter__um)) def test_intra_mu_scatter_keeps_each_motor_unit_anchor(muscle): From e805d299f3e327ce671a2d64d8717d6e8474272f Mon Sep 17 00:00:00 2001 From: Devon Rohlf Date: Sat, 31 Oct 2026 13:19:14 +0100 Subject: [PATCH 11/42] feat(emg): published two-layer end-plate coefficients as an option Konstantin et al. (2020, eq. 5-6) print a_mu = 1, b_mu = 2.5, a_sigma = 0.25, b_sigma = 1 mm; the reference simulator code (and MyoGen's default) uses 1.5, 4, 0.5, 1.5. Offer both in nmj_branch_schedule and SurfaceFiberDispersion, default unchanged. --- myogen/simulator/core/emg/nerve_branches.py | 21 ++++++++++++++++--- .../core/emg/surface/fiber_dispersion.py | 13 +++++++++--- tests/test_surface_dispersion.py | 17 +++++++++++++++ 3 files changed, 45 insertions(+), 6 deletions(-) diff --git a/myogen/simulator/core/emg/nerve_branches.py b/myogen/simulator/core/emg/nerve_branches.py index e419269e..2696e356 100644 --- a/myogen/simulator/core/emg/nerve_branches.py +++ b/myogen/simulator/core/emg/nerve_branches.py @@ -17,16 +17,31 @@ from sklearn.cluster import KMeans -def nmj_branch_schedule(recruitment_thresholds: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]: +# (branch intercept, branch slope, arborisation intercept, arborisation slope), mm. +# "reference_code": Konstantin's reference simulator (iemg_simulator, s08_cl_init_muaps.m). +# "published": Konstantin et al. 2020, eq. 5-6 (a_mu = 1, b_mu = 2.5, a_sigma = 0.25, b_sigma = 1). +NMJ_SCHEDULE_COEFFICIENTS = { + "reference_code": (1.5, 4.0, 0.5, 1.5), + "published": (1.0, 2.5, 0.25, 1.0), +} + + +def nmj_branch_schedule( + recruitment_thresholds: np.ndarray, coefficients: str = "reference_code" +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """Branch count and z spreads per motor unit, growing with motor-unit size. Returns ``(n_branches, branches_z_std__mm, arborization_z_std__mm)``. The - spread factor is the share of total threshold below each unit. + spread factor is the share of total threshold below each unit; + ``coefficients`` picks a set from ``NMJ_SCHEDULE_COEFFICIENTS``. """ + if coefficients not in NMJ_SCHEDULE_COEFFICIENTS: + raise ValueError(f"coefficients must be one of {tuple(NMJ_SCHEDULE_COEFFICIENTS)}, got {coefficients!r}") + branch_a, branch_b, arbor_a, arbor_b = NMJ_SCHEDULE_COEFFICIENTS[coefficients] thresholds = np.asarray(recruitment_thresholds, dtype=float) n_branches = (1 + np.round(np.log(thresholds / thresholds[0]))).astype(int) spread = np.concatenate([[0.0], np.cumsum(thresholds)[:-1]]) / np.sum(thresholds) - return n_branches, 1.5 + spread * 4.0, 0.5 + spread * 1.5 + return n_branches, branch_a + spread * branch_b, arbor_a + spread * arbor_b def two_layer_nmj( diff --git a/myogen/simulator/core/emg/surface/fiber_dispersion.py b/myogen/simulator/core/emg/surface/fiber_dispersion.py index 043faa4e..a9eb2034 100644 --- a/myogen/simulator/core/emg/surface/fiber_dispersion.py +++ b/myogen/simulator/core/emg/surface/fiber_dispersion.py @@ -17,7 +17,7 @@ import numpy as np -from ..nerve_branches import nerve_delays__ms, nmj_branch_schedule, two_layer_nmj +from ..nerve_branches import NMJ_SCHEDULE_COEFFICIENTS, nerve_delays__ms, nmj_branch_schedule, two_layer_nmj NMJ_MODELS = ("uniform_mu", "two_layer", "two_layer_intra") @@ -36,7 +36,11 @@ class SurfaceFiberDispersion: ``"two_layer_intra"``: the two-layer scatter re-centred on each unit's undispersed end-plate, isolating the spread within a motor unit. nmj_scatter_scale : float - Multiplies both two-layer z spreads (1 = the intramuscular default). + Multiplies both two-layer z spreads (1 = the coefficients as given). + nmj_coefficients : {"reference_code", "published"} + Two-layer z-spread coefficients: Konstantin's reference simulator code + (the intramuscular default) or the values printed in Konstantin et al. + (2020), eq. 5-6, which give about two-thirds of the spread. nerve_delays : bool Delay each fibre by its axon-branch conduction time (two-layer nerve paths over ``nerve_cvs__m_per_s``). Independent of ``nmj_model``: with @@ -74,6 +78,7 @@ class SurfaceFiberDispersion: nmj_model: Literal["uniform_mu", "two_layer", "two_layer_intra"] = "uniform_mu" nmj_scatter_scale: float = 1.0 + nmj_coefficients: str = "reference_code" nerve_delays: bool = False nerve_cvs__m_per_s: tuple[float, float] = (5.0, 2.0) fiber_cv_inter_mu_spread: float = 0.0 @@ -87,6 +92,8 @@ class SurfaceFiberDispersion: def __post_init__(self): if self.nmj_model not in NMJ_MODELS: raise ValueError(f"nmj_model must be one of {NMJ_MODELS}, got {self.nmj_model!r}") + if self.nmj_coefficients not in NMJ_SCHEDULE_COEFFICIENTS: + raise ValueError(f"nmj_coefficients must be one of {tuple(NMJ_SCHEDULE_COEFFICIENTS)}") if self.nmj_scatter_scale < 0: raise ValueError("nmj_scatter_scale must be non-negative") if len(self.nerve_cvs__m_per_s) != 2 or min(self.nerve_cvs__m_per_s) <= 0: @@ -140,7 +147,7 @@ def realize_fiber_dispersion( cv_rng = np.random.default_rng([config.seed, 1]) needs_nerve = config.nmj_model != "uniform_mu" or config.nerve_delays if needs_nerve: - n_branches, branch_sd, arbor_sd = nmj_branch_schedule(recruitment_thresholds) + n_branches, branch_sd, arbor_sd = nmj_branch_schedule(recruitment_thresholds, config.nmj_coefficients) report = {key: np.full(n_mu, np.nan) for key in ( "zi_mean__mm", "zi_sd__mm", "cv_mean__m_s", "cv_cov", "delay_mean__ms", "delay_sd__ms", diff --git a/tests/test_surface_dispersion.py b/tests/test_surface_dispersion.py index 4c823993..0fa4b812 100644 --- a/tests/test_surface_dispersion.py +++ b/tests/test_surface_dispersion.py @@ -77,6 +77,23 @@ def test_branch_schedule_grows_with_motor_unit_size(): spread = np.concatenate([[0], np.cumsum(thresholds)[:-1]]) / thresholds.sum() np.testing.assert_allclose(branch_sd, 1.5 + 4.0 * spread) np.testing.assert_allclose(arbor_sd, 0.5 + 1.5 * spread) + # Published coefficients (Konstantin et al. 2020, eq. 5-6): a_mu=1, b_mu=2.5, a_sigma=0.25, b_sigma=1 mm + _, branch_sd, arbor_sd = nmj_branch_schedule(thresholds, coefficients="published") + np.testing.assert_allclose(branch_sd, 1.0 + 2.5 * spread) + np.testing.assert_allclose(arbor_sd, 0.25 + 1.0 * spread) + with pytest.raises(ValueError): + nmj_branch_schedule(thresholds, coefficients="other") + + +def test_surface_dispersion_can_use_the_published_coefficients(muscle): + zones = np.zeros(12) + script = _realize(muscle, SurfaceFiberDispersion(nmj_model="two_layer_intra"), zones) + paper = _realize(muscle, SurfaceFiberDispersion(nmj_model="two_layer_intra", nmj_coefficients="published"), + zones) + big = max(range(12), key=lambda i: len(paper.zi__mm[i])) + assert np.std(paper.zi__mm[big]) < np.std(script.zi__mm[big]) + with pytest.raises(ValueError): + SurfaceFiberDispersion(nmj_coefficients="other") def test_nerve_delay_is_path_over_branch_velocity(): From 061ab2bde2620b5b1336e289f153ea9189f8d3e8 Mon Sep 17 00:00:00 2001 From: Devon Rohlf Date: Sat, 31 Oct 2026 13:54:53 +0100 Subject: [PATCH 12/42] feat(emg)!: published two-layer end-plate coefficients by default Intramuscular and surface EMG now use the z spreads printed in Konstantin et al. (2020, eq. 5-6). Within-unit end-plate SD falls from a median of about 1.3 mm to 0.8 mm (largest units about 1.9 mm instead of 2.9 mm), in line with simulator conventions of about 1 mm and the 2-3.5 mm upper bounds from single-fibre and surface-array studies. The reference-code values remain available as coefficients="reference_code". --- .../core/emg/intramuscular/intramuscular_emg.py | 4 ++-- myogen/simulator/core/emg/nerve_branches.py | 2 +- .../core/emg/surface/fiber_dispersion.py | 10 +++++----- tests/test_surface_dispersion.py | 17 +++++++++-------- 4 files changed, 17 insertions(+), 16 deletions(-) diff --git a/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py b/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py index 128bfd2d..f1a2104b 100644 --- a/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py +++ b/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py @@ -268,8 +268,8 @@ def _simulate_neuromuscular_junctions(self, verbose: bool = True) -> None: """ Simulate neuromuscular junction distributions for all motor units. - This implements the logic from s08_cl_init_muaps.m for generating - realistic NMJ branch patterns with size-dependent complexity. + Two-layer branch model of Konstantin et al. (2020) with the z spreads of + their eq. 5-6, growing with motor-unit size. """ if not self._motor_units: raise ValueError("Must call _initialize_motor_units() first") diff --git a/myogen/simulator/core/emg/nerve_branches.py b/myogen/simulator/core/emg/nerve_branches.py index 2696e356..6925b49e 100644 --- a/myogen/simulator/core/emg/nerve_branches.py +++ b/myogen/simulator/core/emg/nerve_branches.py @@ -27,7 +27,7 @@ def nmj_branch_schedule( - recruitment_thresholds: np.ndarray, coefficients: str = "reference_code" + recruitment_thresholds: np.ndarray, coefficients: str = "published" ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """Branch count and z spreads per motor unit, growing with motor-unit size. diff --git a/myogen/simulator/core/emg/surface/fiber_dispersion.py b/myogen/simulator/core/emg/surface/fiber_dispersion.py index a9eb2034..2724b0ab 100644 --- a/myogen/simulator/core/emg/surface/fiber_dispersion.py +++ b/myogen/simulator/core/emg/surface/fiber_dispersion.py @@ -37,10 +37,10 @@ class SurfaceFiberDispersion: undispersed end-plate, isolating the spread within a motor unit. nmj_scatter_scale : float Multiplies both two-layer z spreads (1 = the coefficients as given). - nmj_coefficients : {"reference_code", "published"} - Two-layer z-spread coefficients: Konstantin's reference simulator code - (the intramuscular default) or the values printed in Konstantin et al. - (2020), eq. 5-6, which give about two-thirds of the spread. + nmj_coefficients : {"published", "reference_code"} + Two-layer z-spread coefficients: the values printed in Konstantin et al. + (2020), eq. 5-6 (default, as in intramuscular EMG), or those of + Konstantin's reference simulator code, about 1.5 times wider. nerve_delays : bool Delay each fibre by its axon-branch conduction time (two-layer nerve paths over ``nerve_cvs__m_per_s``). Independent of ``nmj_model``: with @@ -78,7 +78,7 @@ class SurfaceFiberDispersion: nmj_model: Literal["uniform_mu", "two_layer", "two_layer_intra"] = "uniform_mu" nmj_scatter_scale: float = 1.0 - nmj_coefficients: str = "reference_code" + nmj_coefficients: str = "published" nerve_delays: bool = False nerve_cvs__m_per_s: tuple[float, float] = (5.0, 2.0) fiber_cv_inter_mu_spread: float = 0.0 diff --git a/tests/test_surface_dispersion.py b/tests/test_surface_dispersion.py index 0fa4b812..7a29c0c3 100644 --- a/tests/test_surface_dispersion.py +++ b/tests/test_surface_dispersion.py @@ -72,24 +72,25 @@ def test_two_layer_function_reproduces_motor_unit_sim(): def test_branch_schedule_grows_with_motor_unit_size(): thresholds = np.array([1.0, np.e, np.e**2, np.e**3]) + # Default: published coefficients (Konstantin et al. 2020, eq. 5-6): a_mu=1, b_mu=2.5, a_sigma=0.25, b_sigma=1 mm n_branches, branch_sd, arbor_sd = nmj_branch_schedule(thresholds) np.testing.assert_array_equal(n_branches, [1, 2, 3, 4]) spread = np.concatenate([[0], np.cumsum(thresholds)[:-1]]) / thresholds.sum() - np.testing.assert_allclose(branch_sd, 1.5 + 4.0 * spread) - np.testing.assert_allclose(arbor_sd, 0.5 + 1.5 * spread) - # Published coefficients (Konstantin et al. 2020, eq. 5-6): a_mu=1, b_mu=2.5, a_sigma=0.25, b_sigma=1 mm - _, branch_sd, arbor_sd = nmj_branch_schedule(thresholds, coefficients="published") np.testing.assert_allclose(branch_sd, 1.0 + 2.5 * spread) np.testing.assert_allclose(arbor_sd, 0.25 + 1.0 * spread) + # Konstantin's reference simulator code + _, branch_sd, arbor_sd = nmj_branch_schedule(thresholds, coefficients="reference_code") + np.testing.assert_allclose(branch_sd, 1.5 + 4.0 * spread) + np.testing.assert_allclose(arbor_sd, 0.5 + 1.5 * spread) with pytest.raises(ValueError): nmj_branch_schedule(thresholds, coefficients="other") def test_surface_dispersion_can_use_the_published_coefficients(muscle): zones = np.zeros(12) - script = _realize(muscle, SurfaceFiberDispersion(nmj_model="two_layer_intra"), zones) - paper = _realize(muscle, SurfaceFiberDispersion(nmj_model="two_layer_intra", nmj_coefficients="published"), - zones) + script = _realize(muscle, SurfaceFiberDispersion(nmj_model="two_layer_intra", nmj_coefficients="reference_code"), + zones) + paper = _realize(muscle, SurfaceFiberDispersion(nmj_model="two_layer_intra"), zones) big = max(range(12), key=lambda i: len(paper.zi__mm[i])) assert np.std(paper.zi__mm[big]) < np.std(script.zi__mm[big]) with pytest.raises(ValueError): @@ -214,7 +215,7 @@ def test_intra_mu_scatter_keeps_each_motor_unit_anchor(muscle): for mu in range(12): if len(real.zi__mm[mu]) > 1: assert np.mean(real.zi__mm[mu]) == pytest.approx(zones[mu], abs=1e-9) - assert np.std(real.zi__mm[mu]) > 0.3 + assert np.std(real.zi__mm[mu]) > 0.1 def test_nerve_delays_switch_independently_of_end_plates(muscle): From 4d57e08dff5623c9295200824dfc582619b7937f Mon Sep 17 00:00:00 2001 From: Devon Rohlf Date: Sat, 31 Oct 2026 14:12:55 +0100 Subject: [PATCH 13/42] feat(semg): optional tendon-end band with independent fibre ends tendon_end_band__mm places each fibre end independently, uniform over a band about +-mean_fiber_length/2, instead of the length variation that moves both ends together. Off by default; the dispersion realisation now carries each fibre's ends and the surface model reads L1/L2 from them. --- .../core/emg/surface/fiber_dispersion.py | 44 +++++++++++++++---- .../simulator/core/emg/surface/surface_emg.py | 10 ++++- tests/test_surface_dispersion.py | 35 +++++++++++++++ 3 files changed, 78 insertions(+), 11 deletions(-) diff --git a/myogen/simulator/core/emg/surface/fiber_dispersion.py b/myogen/simulator/core/emg/surface/fiber_dispersion.py index 2724b0ab..24cc7543 100644 --- a/myogen/simulator/core/emg/surface/fiber_dispersion.py +++ b/myogen/simulator/core/emg/surface/fiber_dispersion.py @@ -72,6 +72,10 @@ class SurfaceFiberDispersion: the diameter at the mean CV. cv_per_diameter__m_s_per_mm, cv_intercept__m_s : float Slope and intercept of that line. + tendon_end_band__mm : float + Width of the band over which each fibre end is placed, independently + at either end, uniform about ``±mean_fiber_length / 2``. Replaces the + fibre-length variation (which moves both ends together). 0 keeps it. seed : int Seed of the dispersion draws. """ @@ -87,6 +91,7 @@ class SurfaceFiberDispersion: diameter_scales_source: bool = False cv_per_diameter__m_s_per_mm: float = 50.0 cv_intercept__m_s: float = 0.95 + tendon_end_band__mm: float = 0.0 seed: int = 0 def __post_init__(self): @@ -100,6 +105,8 @@ def __post_init__(self): raise ValueError("nerve_cvs__m_per_s must be two positive velocities") if not 0 <= self.fiber_cv_inter_mu_spread < 2: raise ValueError("fiber_cv_inter_mu_spread must be in [0, 2)") + if self.tendon_end_band__mm < 0: + raise ValueError("tendon_end_band__mm must be non-negative") if self.fiber_cv_threshold_slope < 0: raise ValueError("fiber_cv_threshold_slope must be non-negative") if self.fiber_cv_threshold_slope and self.fiber_cv_inter_mu_spread: @@ -112,7 +119,7 @@ def active(self) -> bool: return (self.nmj_model != "uniform_mu" or self.nerve_delays or self.fiber_cv_inter_mu_spread > 0 or self.fiber_cv_threshold_slope > 0 or self.fiber_cv_within_mu_cov > 0 - or self.diameter_scales_source) + or self.diameter_scales_source or self.tendon_end_band__mm > 0) @dataclass @@ -123,6 +130,8 @@ class FiberDispersionRealization: cv__m_s: list[np.ndarray] delay__ms: list[np.ndarray] amplitude_scale: list[np.ndarray] + end_plus__mm: Optional[list[np.ndarray]] = None + end_minus__mm: Optional[list[np.ndarray]] = None report: dict = field(default_factory=dict) @@ -134,28 +143,35 @@ def realize_fiber_dispersion( mean_cv__m_s: float, fiber_lengths__mm: Optional[Sequence[np.ndarray]] = None, end_margin__mm: float = 1.0, + mean_fiber_length__mm: Optional[float] = None, ) -> FiberDispersionRealization: """Draw per-fibre end-plate, nerve delay, CV and source scale for every motor unit. ``mu_zones__mm`` holds each unit's undispersed end-plate (relative to the - fibre centre). With ``fiber_lengths__mm`` given, end-plates are clipped to - ``end_margin__mm`` inside each fibre and the clipped count is reported. + fibre centre). With ``fiber_lengths__mm`` given, each fibre's ends sit at + ``±length / 2`` (or are drawn in the tendon band, which needs + ``mean_fiber_length__mm``), and end-plates are clipped to + ``end_margin__mm`` inside the ends with the clipped count reported. """ n_mu = len(fiber_positions__mm) counts = np.array([len(p) for p in fiber_positions__mm]) nerve_rng = np.random.default_rng([config.seed, 0]) cv_rng = np.random.default_rng([config.seed, 1]) + end_rng = np.random.default_rng([config.seed, 2]) + band = config.tendon_end_band__mm + if band > 0 and (fiber_lengths__mm is None or mean_fiber_length__mm is None): + raise ValueError("tendon_end_band__mm needs fiber_lengths__mm and mean_fiber_length__mm") needs_nerve = config.nmj_model != "uniform_mu" or config.nerve_delays if needs_nerve: n_branches, branch_sd, arbor_sd = nmj_branch_schedule(recruitment_thresholds, config.nmj_coefficients) report = {key: np.full(n_mu, np.nan) for key in ( "zi_mean__mm", "zi_sd__mm", "cv_mean__m_s", "cv_cov", "delay_mean__ms", "delay_sd__ms", - "amplitude_scale_mean", "diameter_mean__um")} + "amplitude_scale_mean", "diameter_mean__um", "end_sd__mm")} report["n_fibers"] = counts report["n_branches"] = np.zeros(n_mu, dtype=int) report["zi_clipped"] = np.zeros(n_mu, dtype=int) - zi_all, cv_all, delay_all, scale_all = [], [], [], [] + zi_all, cv_all, delay_all, scale_all, plus_all, minus_all = [], [], [], [], [], [] rank = np.arange(n_mu) / max(n_mu - 1, 1) rt = np.asarray(recruitment_thresholds, dtype=float) rt = rt / rt.max() @@ -182,11 +198,20 @@ def realize_fiber_dispersion( zi = zi + (nmj_z - np.mean(nmj_z)) if config.nerve_delays: delay = nerve_delays__ms(paths, config.nerve_cvs__m_per_s) - if fiber_lengths__mm is not None and n: - limit = np.asarray(fiber_lengths__mm[mu]) / 2 - end_margin__mm - clipped = np.clip(zi, -limit, limit) + if fiber_lengths__mm is not None: + if band > 0: + end_plus = mean_fiber_length__mm / 2 + end_rng.uniform(-band / 2, band / 2, size=n) + end_minus = -mean_fiber_length__mm / 2 + end_rng.uniform(-band / 2, band / 2, size=n) + else: + end_plus = np.asarray(fiber_lengths__mm[mu]) / 2 + end_minus = -np.asarray(fiber_lengths__mm[mu]) / 2 + clipped = np.clip(zi, end_minus + end_margin__mm, end_plus - end_margin__mm) report["zi_clipped"][mu] = int(np.count_nonzero(clipped != zi)) zi = clipped + plus_all.append(end_plus) + minus_all.append(end_minus) + if n: + report["end_sd__mm"][mu] = np.std(end_plus) cv_mu = mean_cv__m_s * threshold_factor[mu] * (1 - config.fiber_cv_inter_mu_spread / 2 + config.fiber_cv_inter_mu_spread * rank[mu]) @@ -212,4 +237,5 @@ def realize_fiber_dispersion( report["amplitude_scale_mean"][mu] = np.mean(scale) report["diameter_mean__um"][mu] = np.mean(diameter__um) report["config"] = asdict(config) - return FiberDispersionRealization(zi_all, cv_all, delay_all, scale_all, report) + ends = (plus_all, minus_all) if fiber_lengths__mm is not None else (None, None) + return FiberDispersionRealization(zi_all, cv_all, delay_all, scale_all, *ends, report=report) diff --git a/myogen/simulator/core/emg/surface/surface_emg.py b/myogen/simulator/core/emg/surface/surface_emg.py index e36c532f..e8f8a8c2 100644 --- a/myogen/simulator/core/emg/surface/surface_emg.py +++ b/myogen/simulator/core/emg/surface/surface_emg.py @@ -356,6 +356,7 @@ def simulate_muaps( mean_cv__m_s=float(self._mean_conduction_velocity__m_s), fiber_lengths__mm=[self._mean_fiber_length__mm + variation for variation in fiber_length_variations_per_MU], + mean_fiber_length__mm=float(self._mean_fiber_length__mm), ) self.fiber_dispersion_report = dispersion.report @@ -473,7 +474,8 @@ def _process_single_mu( # The fibre must finish inside the second half of the window # (firing at its centre), or the circular grid wraps it. zone = dispersion.zi__mm[MU_index] - reach = fiber_lengths / 2 + np.abs(zone) + 10.0 # mm, plus the IAP tail + reach = np.maximum(dispersion.end_plus__mm[MU_index] - zone, + zone - dispersion.end_minus__mm[MU_index]) + 10.0 # mm, plus the IAP tail finish = dispersion.delay__ms[MU_index] + reach / dispersion.cv__m_s[MU_index] if np.max(finish) >= kernel_length / v_conduction / 2: raise ValueError( @@ -501,7 +503,11 @@ def _process_single_mu( fiber_zone, fiber_cv = innervation_zone, v_conduction # Semi-lengths towards +z (L1) and -z (L2) from the end-plate - L1, L2 = fiber_semi_lengths(fiber_zone, fiber_length__mm) + if dispersion is not None: + L1 = float(dispersion.end_plus__mm[MU_index][fiber_number]) - fiber_zone + L2 = fiber_zone - float(dispersion.end_minus__mm[MU_index][fiber_number]) + else: + L1, L2 = fiber_semi_lengths(fiber_zone, fiber_length__mm) if self._use_unified: from myogen.simulator.core.emg.fiber_simulation import ( diff --git a/tests/test_surface_dispersion.py b/tests/test_surface_dispersion.py index 7a29c0c3..3aa69474 100644 --- a/tests/test_surface_dispersion.py +++ b/tests/test_surface_dispersion.py @@ -226,6 +226,35 @@ def test_nerve_delays_switch_independently_of_end_plates(muscle): assert real.report["delay_sd__ms"][big] == pytest.approx(np.std(real.delay__ms[big])) +def _realize_with_lengths(muscle, config, mean_length=70.0, var=2.8): + counts = np.asarray(muscle.resulting_number_of_innervated_fibers) + positions = [np.asarray(muscle.resulting_fiber_assignment(i)) for i in range(len(counts))] + lengths = [mean_length + np.random.default_rng(i).uniform(-var, var, size=n) for i, n in enumerate(counts)] + real = realize_fiber_dispersion(config, positions, np.asarray(muscle.recruitment_thresholds), + np.zeros(len(counts)), mean_cv__m_s=4.2, fiber_lengths__mm=lengths, + mean_fiber_length__mm=mean_length) + return real, lengths + + +def test_fibre_ends_follow_the_fibre_length_by_default(muscle): + real, lengths = _realize_with_lengths(muscle, SurfaceFiberDispersion(nerve_delays=True)) + for mu, length in enumerate(lengths): + np.testing.assert_array_equal(real.end_plus__mm[mu], length / 2) + np.testing.assert_array_equal(real.end_minus__mm[mu], -length / 2) + + +def test_tendon_band_staggers_each_end_independently(muscle): + real, _ = _realize_with_lengths(muscle, SurfaceFiberDispersion(tendon_end_band__mm=8.0)) + assert SurfaceFiberDispersion(tendon_end_band__mm=8.0).active + plus = np.concatenate(real.end_plus__mm) + minus = np.concatenate(real.end_minus__mm) + assert np.all(np.abs(plus - 35.0) <= 4.0) and np.all(np.abs(minus + 35.0) <= 4.0) + assert np.std(plus) == pytest.approx(8.0 / np.sqrt(12), rel=0.1) + assert abs(np.corrcoef(plus, minus)[0, 1]) < 4 / np.sqrt(len(plus)) # independent within sampling error + big = max(range(12), key=lambda i: len(real.end_plus__mm[i])) + assert real.report["end_sd__mm"][big] == pytest.approx(np.std(real.end_plus__mm[big])) + + def test_realisation_leaves_the_global_generator_untouched(muscle): set_random_seed(3) expected = get_random_generator().uniform(size=4) @@ -268,6 +297,12 @@ def test_dispersion_changes_the_muap_and_is_reported(muscle): assert report["delay_sd__ms"][11] > 0 and report["cv_cov"][11] > 0 +def test_tendon_band_changes_the_muap(muscle): + base = _muaps(_semg(muscle)) + staggered = _muaps(_semg(muscle, SurfaceFiberDispersion(tendon_end_band__mm=10.0))) + assert not np.allclose(base[11], staggered[11]) + + @pytest.mark.usefixtures("ample_session_memory") def test_other_routes_refuse_active_dispersion(muscle): semg = _semg(muscle, SurfaceFiberDispersion(nerve_delays=True)) From 1d95f27dfdb3c91f4cfafab514108ebe086b9d12 Mon Sep 17 00:00:00 2001 From: Devon Rohlf Date: Sat, 31 Oct 2026 14:51:47 +0100 Subject: [PATCH 14/42] feat(semg): spatially correlated tendon ends tendon_end_correlation__mm makes each fibre end follow a smooth Gaussian random field over the cross-section (one per tendon, squared-exponential covariance via random Fourier features) plus a per-fibre jitter, keeping the total SD of the uniform band. Neighbouring fibres then end close to each other. 0 keeps the independent band. Off by default. --- .../core/emg/surface/fiber_dispersion.py | 40 +++++++++++++++++-- tests/test_surface_dispersion.py | 32 +++++++++++++++ 2 files changed, 69 insertions(+), 3 deletions(-) diff --git a/myogen/simulator/core/emg/surface/fiber_dispersion.py b/myogen/simulator/core/emg/surface/fiber_dispersion.py index 24cc7543..f04b42f6 100644 --- a/myogen/simulator/core/emg/surface/fiber_dispersion.py +++ b/myogen/simulator/core/emg/surface/fiber_dispersion.py @@ -76,6 +76,14 @@ class SurfaceFiberDispersion: Width of the band over which each fibre end is placed, independently at either end, uniform about ``±mean_fiber_length / 2``. Replaces the fibre-length variation (which moves both ends together). 0 keeps it. + tendon_end_correlation__mm : float + Correlation length of the fibre-end position across the cross-section. + 0 draws every end independently over the band. Above 0, each end + follows a smooth Gaussian random field (squared-exponential covariance, + one field per tendon) plus ``tendon_end_jitter__mm`` per fibre, with the + total SD of the uniform band, ``tendon_end_band__mm / sqrt(12)``. + tendon_end_jitter__mm : float + SD of each fibre's own offset from the smooth field. seed : int Seed of the dispersion draws. """ @@ -92,6 +100,8 @@ class SurfaceFiberDispersion: cv_per_diameter__m_s_per_mm: float = 50.0 cv_intercept__m_s: float = 0.95 tendon_end_band__mm: float = 0.0 + tendon_end_correlation__mm: float = 0.0 + tendon_end_jitter__mm: float = 0.0 seed: int = 0 def __post_init__(self): @@ -105,8 +115,14 @@ def __post_init__(self): raise ValueError("nerve_cvs__m_per_s must be two positive velocities") if not 0 <= self.fiber_cv_inter_mu_spread < 2: raise ValueError("fiber_cv_inter_mu_spread must be in [0, 2)") - if self.tendon_end_band__mm < 0: - raise ValueError("tendon_end_band__mm must be non-negative") + if min(self.tendon_end_band__mm, self.tendon_end_correlation__mm, self.tendon_end_jitter__mm) < 0: + raise ValueError("tendon_end_band__mm, tendon_end_correlation__mm and tendon_end_jitter__mm " + "must be non-negative") + if self.tendon_end_correlation__mm > 0: + if self.tendon_end_band__mm == 0: + raise ValueError("tendon_end_correlation__mm needs tendon_end_band__mm") + if self.tendon_end_jitter__mm >= self.tendon_end_band__mm / np.sqrt(12): + raise ValueError("tendon_end_jitter__mm must be below the band SD, tendon_end_band__mm / sqrt(12)") if self.fiber_cv_threshold_slope < 0: raise ValueError("fiber_cv_threshold_slope must be non-negative") if self.fiber_cv_threshold_slope and self.fiber_cv_inter_mu_spread: @@ -135,6 +151,15 @@ class FiberDispersionRealization: report: dict = field(default_factory=dict) +def _smooth_field(rng: np.random.Generator, correlation__mm: float, sd__mm: float, n_features: int = 1000): + """Gaussian random field over the cross-section with covariance + ``sd² exp(-r² / (2 correlation²))``, by random Fourier features.""" + k = rng.normal(0.0, 1.0 / correlation__mm, size=(n_features, 2)) + phase = rng.uniform(0.0, 2 * np.pi, size=n_features) + weight = sd__mm * np.sqrt(2.0 / n_features) + return lambda xy: weight * np.cos(xy @ k.T + phase).sum(axis=1) + + def realize_fiber_dispersion( config: SurfaceFiberDispersion, fiber_positions__mm: Sequence[np.ndarray], @@ -161,6 +186,10 @@ def realize_fiber_dispersion( band = config.tendon_end_band__mm if band > 0 and (fiber_lengths__mm is None or mean_fiber_length__mm is None): raise ValueError("tendon_end_band__mm needs fiber_lengths__mm and mean_fiber_length__mm") + end_fields = None + if band > 0 and config.tendon_end_correlation__mm > 0: + field_sd = np.sqrt(band**2 / 12 - config.tendon_end_jitter__mm**2) + end_fields = [_smooth_field(end_rng, config.tendon_end_correlation__mm, field_sd) for _ in range(2)] needs_nerve = config.nmj_model != "uniform_mu" or config.nerve_delays if needs_nerve: n_branches, branch_sd, arbor_sd = nmj_branch_schedule(recruitment_thresholds, config.nmj_coefficients) @@ -199,7 +228,12 @@ def realize_fiber_dispersion( if config.nerve_delays: delay = nerve_delays__ms(paths, config.nerve_cvs__m_per_s) if fiber_lengths__mm is not None: - if band > 0: + if end_fields is not None: + xy = np.asarray(fiber_positions__mm[mu], dtype=float)[:, :2] + jitter = config.tendon_end_jitter__mm + end_plus = mean_fiber_length__mm / 2 + end_fields[0](xy) + end_rng.normal(0, jitter, size=n) + end_minus = -mean_fiber_length__mm / 2 + end_fields[1](xy) + end_rng.normal(0, jitter, size=n) + elif band > 0: end_plus = mean_fiber_length__mm / 2 + end_rng.uniform(-band / 2, band / 2, size=n) end_minus = -mean_fiber_length__mm / 2 + end_rng.uniform(-band / 2, band / 2, size=n) else: diff --git a/tests/test_surface_dispersion.py b/tests/test_surface_dispersion.py index 3aa69474..595a1bd3 100644 --- a/tests/test_surface_dispersion.py +++ b/tests/test_surface_dispersion.py @@ -255,6 +255,38 @@ def test_tendon_band_staggers_each_end_independently(muscle): assert real.report["end_sd__mm"][big] == pytest.approx(np.std(real.end_plus__mm[big])) +def _end_variogram(muscle, real, near__mm, far__mm): + xy = np.concatenate([np.asarray(muscle.resulting_fiber_assignment(i))[:, :2] for i in range(12)]) + plus = np.concatenate(real.end_plus__mm) + i, j = np.triu_indices(len(xy), k=1) + d = np.hypot(*(xy[i] - xy[j]).T) + sq = (plus[i] - plus[j]) ** 2 + return np.mean(sq[d < near__mm]), np.mean(sq[d > far__mm]) + + +def test_correlated_tendon_ends_vary_smoothly_between_neighbours(muscle): + config = SurfaceFiberDispersion(tendon_end_band__mm=8.0, tendon_end_correlation__mm=2.0, + tendon_end_jitter__mm=0.3) + real, _ = _realize_with_lengths(muscle, config) + near, far = _end_variogram(muscle, real, 0.5, 8.0) + assert near < 0.2 * far + independent, _ = _realize_with_lengths(muscle, SurfaceFiberDispersion(tendon_end_band__mm=8.0)) + near_independent, _ = _end_variogram(muscle, independent, 0.5, 8.0) + assert near < 0.2 * near_independent + sd = real.report["end_sd__mm"] + assert np.nanmedian(sd) < np.nanmedian(independent.report["end_sd__mm"]) + + +def test_tendon_end_correlation_needs_a_band_and_room_for_the_jitter(): + with pytest.raises(ValueError, match="tendon_end_band"): + SurfaceFiberDispersion(tendon_end_correlation__mm=2.0) + with pytest.raises(ValueError, match="jitter"): + SurfaceFiberDispersion(tendon_end_band__mm=1.0, tendon_end_correlation__mm=2.0, + tendon_end_jitter__mm=0.5) + with pytest.raises(ValueError, match="non-negative"): + SurfaceFiberDispersion(tendon_end_band__mm=5.0, tendon_end_correlation__mm=-1.0) + + def test_realisation_leaves_the_global_generator_untouched(muscle): set_random_seed(3) expected = get_random_generator().uniform(size=4) From 8d3e811900aaee3de8025dba0b4c7ccaa3dc1a4f Mon Sep 17 00:00:00 2001 From: Devon Rohlf Date: Sat, 31 Oct 2026 15:37:59 +0100 Subject: [PATCH 15/42] docs(muscle): fibre-length defaults are FDI values; var is a uniform half-width --- myogen/simulator/core/muscle/muscle.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/myogen/simulator/core/muscle/muscle.py b/myogen/simulator/core/muscle/muscle.py index 2e0bd91f..38897a5f 100644 --- a/myogen/simulator/core/muscle/muscle.py +++ b/myogen/simulator/core/muscle/muscle.py @@ -114,9 +114,14 @@ class Muscle: Mean conduction velocity in m/s. Default is set to 4.2 m/s as determined by Nishizono et al. 1990 [4]. Experimental range determined by Nishizono et al. 1990 [4] is between 3.2 and 5.0 m/s. mean_fiber_length__mm : float, default=31.7 - Mean fiber length in mm. Default is set to 31.7 mm as determined by Jacobson et al. 1992 [3] (Table 1). + Mean fiber length in mm. Default is the first dorsal interosseous (FDI) value of Jacobson et al. 1992 [3] + (Table 1). var_fiber_length__mm : float, default=2.8 - Fiber length variance in mm. Default is set to 2.8 mm as determined by Jacobson et al. 1992 [3] (Table 1). + Half-width in mm of the uniform fibre-length variation: each fibre is ``mean_fiber_length__mm`` plus a + uniform draw in ``±var_fiber_length__mm``, with its ends at ``±length / 2``, so both ends move together + by up to half this value. Default 2.8 mm, the FDI fibre-length spread of Jacobson et al. 1992 [3] + (Table 1), which most likely describes variation between specimens rather than between fibres of one + muscle (not verified against the table). radius_bone__mm : float, default=1 Bone radius in mm. Default is set to 1 mm. fat_thickness__mm : float, default=0.3 From 0dd6fb363adb425b13c781711ec15a69f760aca3 Mon Sep 17 00:00:00 2001 From: Devon Rohlf Date: Sat, 31 Oct 2026 17:06:39 +0100 Subject: [PATCH 16/42] feat(semg): literature_candidate dispersion preset SurfaceFiberDispersion.literature_candidate() collects the values chosen from the literature: published two-layer end-plates, 5/2 m/s nerve branches, threshold-linear unit CV (Del Vecchio 2017), 5 % within-unit CV and correlated tendon ends (the last two are assumptions). The diameter-scaled source stays off. Opt-in. --- .../core/emg/surface/fiber_dispersion.py | 24 +++++++++++++++++++ tests/test_surface_dispersion.py | 13 ++++++++++ 2 files changed, 37 insertions(+) diff --git a/myogen/simulator/core/emg/surface/fiber_dispersion.py b/myogen/simulator/core/emg/surface/fiber_dispersion.py index f04b42f6..dadd31a2 100644 --- a/myogen/simulator/core/emg/surface/fiber_dispersion.py +++ b/myogen/simulator/core/emg/surface/fiber_dispersion.py @@ -130,6 +130,30 @@ def __post_init__(self): if not 0 <= self.fiber_cv_within_mu_cov < 1 / 3: raise ValueError("fiber_cv_within_mu_cov must be in [0, 1/3)") + @classmethod + def literature_candidate(cls, seed: int = 0) -> "SurfaceFiberDispersion": + """Dispersion with each value taken from the literature or marked as an assumption. + + - End-plates: two-layer within each unit, coefficients of Konstantin et al. + (2020), eq. 5-6. + - Nerve branches: 5 / 2 m/s, the values of Konstantin's reference code. They + keep the per-fibre delay SD (about 0.2 ms) inside the 0-1 ms range of + Stalberg & Karlsson (2001); the paper's 10 / 1 m/s gives about 0.4 ms. + - Unit conduction velocity: linear in recruitment threshold, slope 0.63 + (Del Vecchio et al. 2017). + - Fibre velocity within a unit: CoV 5 %. A modelling assumption; no + measurement was found. + - Tendon ends: 5 mm band (Keenan et al. 2005), correlated over 5 mm with + 0.3 mm per-fibre jitter (assumptions; no measurement of fibre-end stagger + was found). + - Diameter-scaled source: off, pending an amplitude check against published + ranges. + """ + return cls(nmj_model="two_layer_intra", nmj_scatter_scale=1.0, nmj_coefficients="published", + nerve_delays=True, nerve_cvs__m_per_s=(5.0, 2.0), fiber_cv_threshold_slope=0.63, + fiber_cv_within_mu_cov=0.05, tendon_end_band__mm=5.0, tendon_end_correlation__mm=5.0, + tendon_end_jitter__mm=0.3, seed=seed) + @property def active(self) -> bool: return (self.nmj_model != "uniform_mu" or self.nerve_delays diff --git a/tests/test_surface_dispersion.py b/tests/test_surface_dispersion.py index 595a1bd3..50f51601 100644 --- a/tests/test_surface_dispersion.py +++ b/tests/test_surface_dispersion.py @@ -340,3 +340,16 @@ def test_other_routes_refuse_active_dispersion(muscle): semg = _semg(muscle, SurfaceFiberDispersion(nerve_delays=True)) with pytest.raises(ValueError, match="dispersion"): semg.simulate_muaps(n_jobs=1, verbose=False, execution=ExecutionPolicy()) + + +def test_literature_candidate_preset(): + preset = SurfaceFiberDispersion.literature_candidate() + assert preset.active + assert (preset.nmj_model, preset.nmj_scatter_scale, preset.nmj_coefficients) == ("two_layer_intra", 1.0, "published") + assert preset.nerve_delays and preset.nerve_cvs__m_per_s == (5.0, 2.0) + assert (preset.fiber_cv_threshold_slope, preset.fiber_cv_within_mu_cov, preset.fiber_cv_inter_mu_spread) == ( + 0.63, 0.05, 0.0) + assert not preset.diameter_scales_source + assert (preset.tendon_end_band__mm, preset.tendon_end_correlation__mm, preset.tendon_end_jitter__mm) == ( + 5.0, 5.0, 0.3) + assert SurfaceFiberDispersion.literature_candidate(seed=4).seed == 4 From 2ea926dd7e3ad12d9c4fead6de2b9240fccaeab1 Mon Sep 17 00:00:00 2001 From: Devon Rohlf Date: Sat, 31 Oct 2026 17:13:10 +0100 Subject: [PATCH 17/42] feat(iemg): nmj_coefficients selects the end-plate spreads IntramuscularEMG(nmj_coefficients="reference_code") restores the end-plate spreads used up to 0.13.0; the default is the published Konstantin et al. (2020) values. --- .../emg/intramuscular/intramuscular_emg.py | 12 +++++-- tests/test_surface_dispersion.py | 31 +++++++++++++++++++ 2 files changed, 41 insertions(+), 2 deletions(-) diff --git a/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py b/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py index f1a2104b..25be8e24 100644 --- a/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py +++ b/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py @@ -38,7 +38,7 @@ except (ImportError, Exception): HAS_CUPY = False -from ..nerve_branches import nmj_branch_schedule +from ..nerve_branches import NMJ_SCHEDULE_COEFFICIENTS, nmj_branch_schedule from .motor_unit_sim import MotorUnitSim @@ -74,6 +74,10 @@ class IntramuscularEMG: 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 Konstantin et al. 2020 [1]. + nmj_coefficients : {"published", "reference_code"}, default="published" + z spreads of the two-layer end-plate model: the values printed in + Konstantin et al. 2020 [1], eq. 5-6, or those of their reference simulator + code (about 1.5 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. Default is None. For computational efficiency, consider @@ -109,7 +113,10 @@ def __init__( 2.0 * pq.m / pq.s, ), MUs_to_simulate: list[int] | None = None, + nmj_coefficients: str = "published", ): + if nmj_coefficients not in NMJ_SCHEDULE_COEFFICIENTS: + raise ValueError(f"nmj_coefficients must be one of {tuple(NMJ_SCHEDULE_COEFFICIENTS)}") # Immutable public arguments - never modify these self.muscle_model = muscle_model self.electrode_array = electrode_array @@ -119,6 +126,7 @@ def __init__( self.nmj_jitter__s = nmj_jitter__s self.branch_cvs__m_per_s = branch_cvs__m_per_s self.MUs_to_simulate = MUs_to_simulate + self.nmj_coefficients = nmj_coefficients # Private copies for internal modifications (extract magnitudes) self._muscle_model = muscle_model @@ -275,7 +283,7 @@ def _simulate_neuromuscular_junctions(self, verbose: bool = True) -> None: raise ValueError("Must call _initialize_motor_units() first") n_branches, branches_z_std, arborization_z_std = nmj_branch_schedule( - self._muscle_model.recruitment_thresholds + self._muscle_model.recruitment_thresholds, self.nmj_coefficients ) for mu_idx, mu_sim in enumerate( diff --git a/tests/test_surface_dispersion.py b/tests/test_surface_dispersion.py index 50f51601..f2227841 100644 --- a/tests/test_surface_dispersion.py +++ b/tests/test_surface_dispersion.py @@ -353,3 +353,34 @@ def test_literature_candidate_preset(): assert (preset.tendon_end_band__mm, preset.tendon_end_correlation__mm, preset.tendon_end_jitter__mm) == ( 5.0, 5.0, 0.3) assert SurfaceFiberDispersion.literature_candidate(seed=4).seed == 4 + + +def test_intramuscular_emg_selects_the_end_plate_coefficients(monkeypatch): + from myogen.simulator import IntramuscularElectrodeArray, IntramuscularEMG + from myogen.simulator.core.emg.intramuscular.motor_unit_sim import MotorUnitSim + + set_random_seed(5) + thresholds, _ = RecruitmentThresholds(N=4, recruitment_range__ratio=30) + muscle = Muscle(recruitment_thresholds=thresholds, radius__mm=2.0 * pq.mm, radius_bone__mm=0.5 * pq.mm, + fiber_density__fibers_per_mm2=50 * pq.mm**-2, grid_resolution=32, autorun=True) + electrode = IntramuscularElectrodeArray(num_electrodes=2, inter_electrode_distance__mm=0.5 * pq.mm, + position__mm=(0.0 * pq.mm, 0.0 * pq.mm, 15.0 * pq.mm)) + seen = [] + original = MotorUnitSim.sim_nmj_branches_two_layers + + def capture(self, **kwargs): + seen.append((kwargs["branches_z_std"], kwargs["arborization_z_std"])) + return original(self, **kwargs) + + monkeypatch.setattr(MotorUnitSim, "sim_nmj_branches_two_layers", capture) + for coefficients in ("published", "reference_code"): + seen.clear() + kwargs = {} if coefficients == "published" else {"nmj_coefficients": coefficients} + iemg = IntramuscularEMG(muscle_model=muscle, electrode_array=electrode, **kwargs) + iemg._initialize_motor_units(verbose=False) + iemg._simulate_neuromuscular_junctions(verbose=False) + _, branch_sd, arbor_sd = nmj_branch_schedule(np.asarray(muscle.recruitment_thresholds), coefficients) + np.testing.assert_allclose([s[0] for s in seen], branch_sd) + np.testing.assert_allclose([s[1] for s in seen], arbor_sd) + with pytest.raises(ValueError, match="nmj_coefficients"): + IntramuscularEMG(muscle_model=muscle, electrode_array=electrode, nmj_coefficients="paper") From b269c24013020bd77a40903618cea8c99d6fbd4e Mon Sep 17 00:00:00 2001 From: Devon Rohlf Date: Sat, 31 Oct 2026 17:25:07 +0100 Subject: [PATCH 18/42] docs(semg): dispersion and version-comparison harnesses, changelog benchmarks/semg_dispersion: paired sweeps of the SurfaceFiberDispersion options on a limb-like cylinder and the default (FDI) muscle; arms run before the published end-plate default pin reference_code coefficients and reproduce the stored outputs bit for bit. benchmarks/semg_versions: the before/after comparison of the surface-EMG fixes. CHANGELOG: fibre dispersion, the literature preset, the published end-plate spreads as iEMG default (nmj_coefficients restores the old ones) and the fibre-length docstring. Coefficient wording corrected to 1.5-2x; preset docstring drops an unverified delay comparison. --- CHANGELOG.md | 16 + benchmarks/semg_dispersion/README.md | 57 ++++ benchmarks/semg_dispersion/arms.json | 316 ++++++++++++++++++ benchmarks/semg_dispersion/gen_arm.py | 127 +++++++ benchmarks/semg_dispersion/plot_amplitude.py | 134 ++++++++ benchmarks/semg_dispersion/plot_signals.py | 142 ++++++++ benchmarks/semg_dispersion/run_arms.sh | 26 ++ benchmarks/semg_dispersion/score.py | 209 ++++++++++++ benchmarks/semg_versions/README.md | 15 + benchmarks/semg_versions/gen_real.py | 107 ++++++ benchmarks/semg_versions/plot_real.py | 173 ++++++++++ .../emg/intramuscular/intramuscular_emg.py | 2 +- .../core/emg/surface/fiber_dispersion.py | 8 +- 13 files changed, 1327 insertions(+), 5 deletions(-) create mode 100644 benchmarks/semg_dispersion/README.md create mode 100644 benchmarks/semg_dispersion/arms.json create mode 100644 benchmarks/semg_dispersion/gen_arm.py create mode 100644 benchmarks/semg_dispersion/plot_amplitude.py create mode 100644 benchmarks/semg_dispersion/plot_signals.py create mode 100644 benchmarks/semg_dispersion/run_arms.sh create mode 100644 benchmarks/semg_dispersion/score.py create mode 100644 benchmarks/semg_versions/README.md create mode 100644 benchmarks/semg_versions/gen_real.py create mode 100644 benchmarks/semg_versions/plot_real.py diff --git a/CHANGELOG.md b/CHANGELOG.md index b9aba031..b825b144 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -25,6 +25,21 @@ This release corrects the surface-EMG pipeline against Farina, Mesin, Martina & ### Added - `SurfaceEMG(fiber_diameter__um=50, intracellular_conductivity__S_m=1.01)` set the source strength; amplitude scales with the fibre cross-section. +- **Fibre dispersion in surface EMG.** `SurfaceEMG(fiber_dispersion=SurfaceFiberDispersion(...))` adds, each switched on independently: + - per-fibre end-plates from the two-layer nerve model of Konstantin et al. (2020), shared with intramuscular EMG; + - axon-branch conduction delays; + - motor-unit conduction velocity linear in recruitment threshold (Del Vecchio et al. 2017), and a spread of fibre velocities within each unit, each fibre propagating at its own velocity; + - a source scaled by fibre diameter (from the Nandedkar & Stålberg 1983 velocity–diameter line); + - tendon-end stagger: each fibre end in a band, independently or following a smooth random field across the cross-section so that neighbouring fibres end close together. + + The realised spreads are reported per motor unit in `SurfaceEMG.fiber_dispersion_report`. The default is off and leaves outputs bit-identical. `SurfaceFiberDispersion.literature_candidate()` gives a preset whose values come from the literature or are marked as assumptions. The `execution=` route and `use_unified=True` refuse an active configuration. +- `IntramuscularEMG(nmj_coefficients=...)` selects the end-plate spreads (see Changed). +- `benchmarks/semg_dispersion/`: paired sweeps of these options on a limb-like cylinder and on the default (FDI) muscle. + +### Changed + +- **End-plate spreads in intramuscular EMG.** The two-layer end-plate model now uses the spreads printed in Konstantin et al. (2020, eq. 5-6) instead of the 1.5–2× wider values of their reference code: branch-point SD 1.0 + 2.5·s mm (was 1.5 + 4.0·s) and arborisation SD 0.25 + 1.0·s mm (was 0.5 + 1.5·s), with s the share of total threshold below the unit. In the sEMG limb sweep the realised within-unit end-plate SD fell from a median of 1.68 mm to 1.07 mm. iEMG MUAPs and signals change. +- `Muscle.var_fiber_length__mm` is documented as the half-width of a uniform fibre-length variation (the FDI value of Jacobson et al. 1992), not a variance. Behaviour is unchanged. ### Upgrade notes @@ -32,6 +47,7 @@ This release corrects the surface-EMG pipeline against Farina, Mesin, Martina & - 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 dry-skin conductor, pass `skin_conductivity__S_per_m=4.88e-4 * pq.S / pq.m` to `Muscle`. - Pickled `SurfaceEMG` objects from earlier versions lack the new fibre properties; create them again. - The `execution=` MUAP disk-cache key changes, so templates cached by earlier versions are recomputed. +- iEMG outputs change through the narrower end-plate spreads. Pass `IntramuscularEMG(..., nmj_coefficients="reference_code")` to reproduce earlier end-plates. ## [0.13.0] - 2026-10-28 diff --git a/benchmarks/semg_dispersion/README.md b/benchmarks/semg_dispersion/README.md new file mode 100644 index 00000000..5fdb1a77 --- /dev/null +++ b/benchmarks/semg_dispersion/README.md @@ -0,0 +1,57 @@ +# Surface-EMG fibre-dispersion sweep + +Paired simulations behind the `SurfaceFiberDispersion` options and the +`SurfaceFiberDispersion.literature_candidate()` preset. Each arm reuses the same +muscle, fibres, spike trains and electrode as its scenario's baseline and +changes only the dispersion configuration, so MUAPs are compared unit by unit. +Fibre and spike hashes are stored in every output and checked when scoring. + +## Scenarios + +| Scenario | Muscle | Electrode | +|---|---|---| +| `limb` | cylinder r 15 mm, bone 3 mm, fat 4 mm, 70 ± 5 mm fibres, 350 fibres/mm², 100 MUs | bipolar, 5 mm contacts, IED 10 mm, 20 mm from the end-plate band | +| `fdi` | default `Muscle` (tuned for the first dorsal interosseous), 120 MUs | bipolar, 2 mm contacts, IED 5 mm, 8 mm from the end-plate band | + +Both: 30 % excitation trapezoid, Fuglevand thresholds (50-fold range) recruited up +to 70 % excitation, 8–35 Hz rate coding, ISI CV 0.2, 2048 Hz. Scoring applies a +20–450 Hz band-pass and no noise. + +## Arms + +`arms.json` holds every arm per scenario (`baseline` is the undispersed model). +Arms run before the published end-plate coefficients became the default pin +`"nmj_coefficients": "reference_code"`, so they reproduce. `_skin_thickness__mm` +is a harness key for the skin layer. + +## Running + +```bash +# from the repository root; PY selects the interpreter. Check that the +# recorded myogen path (key "myogen" in each output) is this checkout; set +# PYTHONPATH to the checkout if another one is installed. +PY=python bash benchmarks/semg_dispersion/run_arms.sh fdi out/fdi # all FDI arms +PY=python bash benchmarks/semg_dispersion/run_arms.sh limb out/limb baseline candidate_published +python benchmarks/semg_dispersion/score.py fdi out/fdi # scorecard.tsv, semg_dispersion_sweep.pdf +python benchmarks/semg_dispersion/plot_signals.py fdi out/fdi # semg_dispersion_signals.pdf +python benchmarks/semg_dispersion/plot_amplitude.py out/limb out/fdi out/amplitude.pdf +``` + +A limb arm takes about 20–25 min on CPU and an FDI arm about 5 min; the runner +does three at a time and skips arms already present. + +## Findings (2026-10-31) + +- Within-unit fibre CV spread is the strongest lever on MUAP peak-to-peak; a + CV difference between units only rescales time. +- In the limb scenario, dispersion mostly blunts the end-of-fibre phase of small + and mid units; plateau RMS falls by about 6 %. In the FDI scenario the short, + sharp MUAPs of large superficial units lose more, and plateau RMS falls by + about 17 %. +- Tendon ends: a 5 mm band of independent ends lowers median MUAP + peak-to-peak by about 17 % (limb); ends correlated over 5 mm across the + cross-section, so that neighbours end close together, by about 11 %. +- FDI skin thickness 1.29 → 3 mm lowers RMS by about 60 %; it is the strongest + amplitude lever found so far. + +`benchmarks/semg_versions` compares outputs before and after the unreleased surface-EMG fixes. diff --git a/benchmarks/semg_dispersion/arms.json b/benchmarks/semg_dispersion/arms.json new file mode 100644 index 00000000..856d0a23 --- /dev/null +++ b/benchmarks/semg_dispersion/arms.json @@ -0,0 +1,316 @@ +{ + "limb": { + "baseline": {}, + "nerve_cv_10_4": { + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 10.0, + 4.0 + ], + "nmj_coefficients": "reference_code" + }, + "nerve_cv_5_2": { + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 5.0, + 2.0 + ], + "nmj_coefficients": "reference_code" + }, + "nerve_cv_2.5_1": { + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 2.5, + 1.0 + ], + "nmj_coefficients": "reference_code" + }, + "nmj_intra_x0.5": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 0.5, + "nmj_coefficients": "reference_code" + }, + "nmj_intra_x1": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 1.0, + "nmj_coefficients": "reference_code" + }, + "nmj_intra_x2": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 2.0, + "nmj_coefficients": "reference_code" + }, + "nmj_two_layer_x1": { + "nmj_model": "two_layer", + "nmj_scatter_scale": 1.0, + "nmj_coefficients": "reference_code" + }, + "fcv_within_5": { + "fiber_cv_within_mu_cov": 0.05 + }, + "fcv_within_10": { + "fiber_cv_within_mu_cov": 0.1 + }, + "fcv_inter_30": { + "fiber_cv_inter_mu_spread": 0.3 + }, + "fcv_both": { + "fiber_cv_inter_mu_spread": 0.3, + "fiber_cv_within_mu_cov": 0.05 + }, + "fcv_both_diameter": { + "fiber_cv_inter_mu_spread": 0.3, + "fiber_cv_within_mu_cov": 0.05, + "diameter_scales_source": true + }, + "combined": { + "nmj_model": "two_layer", + "nerve_delays": true, + "fiber_cv_inter_mu_spread": 0.3, + "fiber_cv_within_mu_cov": 0.05, + "nmj_coefficients": "reference_code" + }, + "combined_diameter": { + "nmj_model": "two_layer", + "nerve_delays": true, + "fiber_cv_inter_mu_spread": 0.3, + "fiber_cv_within_mu_cov": 0.05, + "diameter_scales_source": true, + "nmj_coefficients": "reference_code" + }, + "candidate": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 1.0, + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 5.0, + 2.0 + ], + "fiber_cv_threshold_slope": 0.63, + "nmj_coefficients": "reference_code" + }, + "candidate_within5": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 1.0, + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 5.0, + 2.0 + ], + "fiber_cv_threshold_slope": 0.63, + "fiber_cv_within_mu_cov": 0.05, + "nmj_coefficients": "reference_code" + }, + "candidate_published": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 1.0, + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 5.0, + 2.0 + ], + "fiber_cv_threshold_slope": 0.63, + "fiber_cv_within_mu_cov": 0.05, + "nmj_coefficients": "published" + }, + "tendon_w5": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 1.0, + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 5.0, + 2.0 + ], + "fiber_cv_threshold_slope": 0.63, + "fiber_cv_within_mu_cov": 0.05, + "nmj_coefficients": "published", + "tendon_end_band__mm": 5.0 + }, + "tendon_w10": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 1.0, + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 5.0, + 2.0 + ], + "fiber_cv_threshold_slope": 0.63, + "fiber_cv_within_mu_cov": 0.05, + "nmj_coefficients": "published", + "tendon_end_band__mm": 10.0 + }, + "tendon_w5_corr2": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 1.0, + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 5.0, + 2.0 + ], + "fiber_cv_threshold_slope": 0.63, + "fiber_cv_within_mu_cov": 0.05, + "nmj_coefficients": "published", + "tendon_end_band__mm": 5.0, + "tendon_end_correlation__mm": 2.0, + "tendon_end_jitter__mm": 0.3 + }, + "tendon_w5_corr5": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 1.0, + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 5.0, + 2.0 + ], + "fiber_cv_threshold_slope": 0.63, + "fiber_cv_within_mu_cov": 0.05, + "nmj_coefficients": "published", + "tendon_end_band__mm": 5.0, + "tendon_end_correlation__mm": 5.0, + "tendon_end_jitter__mm": 0.3 + }, + "tendon_w5_corr10": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 1.0, + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 5.0, + 2.0 + ], + "fiber_cv_threshold_slope": 0.63, + "fiber_cv_within_mu_cov": 0.05, + "nmj_coefficients": "published", + "tendon_end_band__mm": 5.0, + "tendon_end_correlation__mm": 10.0, + "tendon_end_jitter__mm": 0.3 + }, + "full_diameter": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 1.0, + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 5.0, + 2.0 + ], + "fiber_cv_threshold_slope": 0.63, + "fiber_cv_within_mu_cov": 0.05, + "nmj_coefficients": "published", + "tendon_end_band__mm": 5.0, + "tendon_end_correlation__mm": 5.0, + "tendon_end_jitter__mm": 0.3, + "diameter_scales_source": true + } + }, + "fdi": { + "baseline": {}, + "candidate_published": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 1.0, + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 5.0, + 2.0 + ], + "fiber_cv_threshold_slope": 0.63, + "fiber_cv_within_mu_cov": 0.05, + "nmj_coefficients": "published" + }, + "tendon_w5": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 1.0, + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 5.0, + 2.0 + ], + "fiber_cv_threshold_slope": 0.63, + "fiber_cv_within_mu_cov": 0.05, + "nmj_coefficients": "published", + "tendon_end_band__mm": 5.0 + }, + "tendon_w5_corr5": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 1.0, + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 5.0, + 2.0 + ], + "fiber_cv_threshold_slope": 0.63, + "fiber_cv_within_mu_cov": 0.05, + "nmj_coefficients": "published", + "tendon_end_band__mm": 5.0, + "tendon_end_correlation__mm": 5.0, + "tendon_end_jitter__mm": 0.3 + }, + "tendon_w10": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 1.0, + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 5.0, + 2.0 + ], + "fiber_cv_threshold_slope": 0.63, + "fiber_cv_within_mu_cov": 0.05, + "nmj_coefficients": "published", + "tendon_end_band__mm": 10.0 + }, + "full_diameter": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 1.0, + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 5.0, + 2.0 + ], + "fiber_cv_threshold_slope": 0.63, + "fiber_cv_within_mu_cov": 0.05, + "nmj_coefficients": "published", + "tendon_end_band__mm": 5.0, + "tendon_end_correlation__mm": 5.0, + "tendon_end_jitter__mm": 0.3, + "diameter_scales_source": true + }, + "baseline_skin2": { + "_skin_thickness__mm": 2.0 + }, + "full_skin2": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 1.0, + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 5.0, + 2.0 + ], + "fiber_cv_threshold_slope": 0.63, + "fiber_cv_within_mu_cov": 0.05, + "nmj_coefficients": "published", + "tendon_end_band__mm": 5.0, + "tendon_end_correlation__mm": 5.0, + "tendon_end_jitter__mm": 0.3, + "diameter_scales_source": true, + "_skin_thickness__mm": 2.0 + }, + "baseline_skin3": { + "_skin_thickness__mm": 3.0 + }, + "full_skin3": { + "nmj_model": "two_layer_intra", + "nmj_scatter_scale": 1.0, + "nerve_delays": true, + "nerve_cvs__m_per_s": [ + 5.0, + 2.0 + ], + "fiber_cv_threshold_slope": 0.63, + "fiber_cv_within_mu_cov": 0.05, + "nmj_coefficients": "published", + "tendon_end_band__mm": 5.0, + "tendon_end_correlation__mm": 5.0, + "tendon_end_jitter__mm": 0.3, + "diameter_scales_source": true, + "_skin_thickness__mm": 3.0 + } + } +} \ No newline at end of file diff --git a/benchmarks/semg_dispersion/gen_arm.py b/benchmarks/semg_dispersion/gen_arm.py new file mode 100644 index 00000000..5e6c71f9 --- /dev/null +++ b/benchmarks/semg_dispersion/gen_arm.py @@ -0,0 +1,127 @@ +"""Simulate one sweep arm: a frozen sEMG scenario with a fibre-dispersion configuration. + +Usage: python gen_arm.py