Skip to content

Simulator

Recruitment

RecruitmentThresholds

RecruitmentThresholds(N: int, recruitment_range__ratio: float, deluca__slope: float | None = None, konstantin__max_threshold__ratio: float = 1.0, mode: Literal['fuglevand', 'deluca', 'konstantin', 'combined'] = 'konstantin')

Motor unit recruitment threshold generator using physiological models.

This class computes recruitment thresholds for motor unit pools using different established models from the literature.

Generate recruitment thresholds for a pool of motor units using different models.

This class computes the recruitment thresholds (and zero-based thresholds) for a pool of N motor units according to one of several models from the literature. The distribution of thresholds is controlled by the recruitment range (RR) and, for some models, additional parameters.

Following models are available:
- Fuglevand et al. (1993) [1] - De Luca & Contessa (2012) [2] - Konstantin et al. (2020) [3] - Combined model

Parameters:

Name Type Description Default
N int

Number of motor units in the pool.

required
recruitment_range__ratio float

Recruitment range (dimensionless ratio), defined as the ratio of the largest to smallest threshold \((rt(N)/rt(1))\).

required
deluca__slope float

Dimensionless slope parameter for the 'deluca' mode. Required if mode='deluca'. Controls the curvature of the threshold distribution. Typical values range from 0.001-100.

None
konstantin__max_threshold__ratio float

Maximum recruitment threshold (dimensionless ratio) for the 'konstantin' mode. Required if mode='konstantin'. Sets the absolute scale of all thresholds. Default is 1.0.

1.0
mode RecruitmentMode

Model to use for threshold generation. One of 'fuglevand', 'deluca', 'konstantin', or 'combined'. Default is 'konstantin'.

'konstantin'

Attributes:

Name Type Description
rt RECRUITMENT_THRESHOLDS__ARRAY

Recruitment thresholds for each motor unit (shape: (N,)). Values are monotonically increasing from rt[0] to rt[N-1].

rtz RECRUITMENT_THRESHOLDS__ARRAY

Zero-based recruitment thresholds where \(rtz[0] = 0\) (shape: (N,)). Computed as \(rtz = rt - rt[0]\), convenient for simulation.

Raises:

Type Description
ValueError

If a required mode-specific parameter is not provided or if an unknown mode is specified.

References

[1] Fuglevand, A.J., Winter, D.A., Patla, A.E., 1993. Models of recruitment and rate coding organization in motor-unit pools. Journal of Neurophysiology 70, 2470-2488. https://doi.org/10.1152/jn.1993.70.6.2470
[2] De Luca, C.J., Contessa, P., 2012. Hierarchical control of motor units in voluntary contractions. Journal of Neurophysiology 107, 178-195. https://doi.org/10.1152/jn.00961.2010
[3] Konstantin, A., Yu, T., Le Carpentier, E., Aoustin, Y., Farina, D., 2020. Simulation of Motor Unit Action Potential Recordings From Intramuscular Multichannel Scanning Electrodes. IEEE Transactions on Biomedical Engineering 67, 2005-2014. https://doi.org/10.1109/TBME.2019.2953680

Notes

fuglevand : Fuglevand et al. (1993) [1] exponential model

\[rt(i) = \exp\left( \frac{i \cdot \ln(RR)}{N} \right) / 100\]

where \(i = 1, 2, \ldots, N\)

deluca : De Luca & Contessa (2012) [2] model with slope correction

\[rt(i) = \frac{b \cdot i}{N} \cdot \exp\left(\frac{i \cdot \ln(RR / b)}{N}\right) / 100\]

where \(b\) = deluca__slope, \(i = 1, 2, \ldots, N\)

konstantin : Konstantin et al. (2020) [3] model allowing explicit maximum threshold control

\[ \begin{aligned} rt(i) &= \frac{RT_{max}}{RR} \cdot \exp\left(\frac{(i - 1) \cdot \ln(RR)}{N - 1}\right) \\ rtz(i) &= \frac{RT_{max}}{RR} \cdot \left(\exp\left(\frac{(i - 1) \cdot \ln(RR + 1)}{N}\right) - 1\right) \end{aligned} \]

where \(RT_{max}\) = konstantin__max_threshold__ratio, \(i = 1, 2, \ldots, N\)

combined : A corrected De Luca model that uses the slope parameter for shape control but properly respects the RR constraint and maximum threshold like the Konstantin model

\[rt(i) = \frac{RT_{max}}{RR} + \left(\frac{b \cdot i}{N} \cdot \exp\left(\frac{i \cdot \ln(RR / b)}{N}\right) - \frac{RT_{max}}{RR}\right) \cdot \left(\frac{RT_{max} - RT_{max}/RR}{b \cdot N \cdot \exp\left(\frac{i \cdot \ln(RR / b)}{N}\right) - \frac{RT_{max}}{RR}}\right)\]

where \(b\) = deluca__slope, \(RT_{max}\) = konstantin__max_threshold__ratio, \(i = 1, 2, \ldots, N\)

Examples:

>>> # Generate thresholds using Fuglevand model
>>> thresholds = RecruitmentThresholds(
...     N=100, recruitment_range__ratio=50.0, mode='fuglevand'
... )
>>> rt, rtz = thresholds  # Tuple unpacking works
>>> # Or access directly
>>> rt = thresholds.rt
>>> rtz = thresholds.rtz
>>>
>>> # Generate thresholds using Konstantin model with explicit max threshold
>>> thresholds = RecruitmentThresholds(
...     N=100, recruitment_range__ratio=50.0, konstantin__max_threshold__ratio=1.0, mode='konstantin'
... )
>>> rt, rtz = thresholds
Source code in myogen/simulator/core/physiological_distribution.py
def __init__(
    self,
    N: int,
    recruitment_range__ratio: float,
    deluca__slope: float | None = None,
    konstantin__max_threshold__ratio: float = 1.0,
    mode: Literal["fuglevand", "deluca", "konstantin", "combined"] = "konstantin",
) -> None:
    r"""
    Generate recruitment thresholds for a pool of motor units using different models.

    This class computes the recruitment thresholds (and zero-based thresholds) for a pool of N motor units
    according to one of several models from the literature. The distribution of thresholds is controlled by the
    recruitment range (RR) and, for some models, additional parameters.

    Following models are available:  
        - Fuglevand et al. (1993) [1]
        - De Luca & Contessa (2012) [2]
        - Konstantin et al. (2020) [3]
        - Combined model

    Parameters
    ----------
    N : int
        Number of motor units in the pool.
    recruitment_range__ratio : float
        Recruitment range (dimensionless ratio), defined as the ratio of the largest to smallest threshold
        $(rt(N)/rt(1))$.
    deluca__slope : float, optional
        Dimensionless slope parameter for the ``'deluca'`` mode. Required if ``mode='deluca'``.
        Controls the curvature of the threshold distribution. Typical values range from 0.001-100.
    konstantin__max_threshold__ratio : float, optional
        Maximum recruitment threshold (dimensionless ratio) for the ``'konstantin'`` mode. Required if ``mode='konstantin'``.
        Sets the absolute scale of all thresholds. Default is 1.0.
    mode : RecruitmentMode, optional
        Model to use for threshold generation. One of ``'fuglevand'``, ``'deluca'``, ``'konstantin'``, or ``'combined'``.
        Default is ``'konstantin'``.

    Attributes
    ----------
    rt : RECRUITMENT_THRESHOLDS__ARRAY
        Recruitment thresholds for each motor unit (shape: (N,)).
        Values are monotonically increasing from ``rt[0]`` to ``rt[N-1]``.
    rtz : RECRUITMENT_THRESHOLDS__ARRAY
        Zero-based recruitment thresholds where $rtz[0] = 0$ (shape: (N,)).
        Computed as $rtz = rt - rt[0]$, convenient for simulation.

    Raises
    ------
    ValueError
        If a required mode-specific parameter is not provided or if an unknown mode is specified.

    References
    ----------
    [1] Fuglevand, A.J., Winter, D.A., Patla, A.E., 1993. Models of recruitment and rate coding organization in motor-unit pools. Journal of Neurophysiology 70, 2470-2488. https://doi.org/10.1152/jn.1993.70.6.2470 <br>
    [2] De Luca, C.J., Contessa, P., 2012. Hierarchical control of motor units in voluntary contractions. Journal of Neurophysiology 107, 178-195. https://doi.org/10.1152/jn.00961.2010 <br>
    [3] Konstantin, A., Yu, T., Le Carpentier, E., Aoustin, Y., Farina, D., 2020. Simulation of Motor Unit Action Potential Recordings From Intramuscular Multichannel Scanning Electrodes. IEEE Transactions on Biomedical Engineering 67, 2005-2014. https://doi.org/10.1109/TBME.2019.2953680

    Notes
    -----
    **fuglevand** : Fuglevand et al. (1993) [1] exponential model

    $$rt(i) = \exp\left( \frac{i \cdot \ln(RR)}{N} \right) / 100$$

    where $i = 1, 2, \ldots, N$

    **deluca** : De Luca & Contessa (2012) [2] model with slope correction

    $$rt(i) = \frac{b \cdot i}{N} \cdot \exp\left(\frac{i \cdot \ln(RR / b)}{N}\right) / 100$$

    where $b$ = ``deluca__slope``, $i = 1, 2, \ldots, N$

    **konstantin** : Konstantin et al. (2020) [3] model allowing explicit maximum threshold control

    $$
    \begin{aligned}
    rt(i) &= \frac{RT_{max}}{RR} \cdot \exp\left(\frac{(i - 1) \cdot \ln(RR)}{N - 1}\right) \\
    rtz(i) &= \frac{RT_{max}}{RR} \cdot \left(\exp\left(\frac{(i - 1) \cdot \ln(RR + 1)}{N}\right) - 1\right)
    \end{aligned}
    $$

    where $RT_{max}$ = ``konstantin__max_threshold__ratio``, $i = 1, 2, \ldots, N$

    **combined** : A corrected De Luca model that uses the slope parameter for shape control but properly respects the RR constraint and maximum threshold like the Konstantin model

    $$rt(i) = \frac{RT_{max}}{RR} + \left(\frac{b \cdot i}{N} \cdot \exp\left(\frac{i \cdot \ln(RR / b)}{N}\right) - \frac{RT_{max}}{RR}\right) \cdot \left(\frac{RT_{max} - RT_{max}/RR}{b \cdot N \cdot \exp\left(\frac{i \cdot \ln(RR / b)}{N}\right) - \frac{RT_{max}}{RR}}\right)$$

    where $b$ = ``deluca__slope``, $RT_{max}$ = ``konstantin__max_threshold__ratio``, $i = 1, 2, \ldots, N$

    Examples
    --------
    >>> # Generate thresholds using Fuglevand model
    >>> thresholds = RecruitmentThresholds(
    ...     N=100, recruitment_range__ratio=50.0, mode='fuglevand'
    ... )
    >>> rt, rtz = thresholds  # Tuple unpacking works
    >>> # Or access directly
    >>> rt = thresholds.rt
    >>> rtz = thresholds.rtz
    >>>
    >>> # Generate thresholds using Konstantin model with explicit max threshold
    >>> thresholds = RecruitmentThresholds(
    ...     N=100, recruitment_range__ratio=50.0, konstantin__max_threshold__ratio=1.0, mode='konstantin'
    ... )
    >>> rt, rtz = thresholds
    """
    # Store immutable public parameters (for joblib serialization)
    self.N = N
    self.recruitment_range__ratio = recruitment_range__ratio
    self.deluca__slope = deluca__slope
    self.konstantin__max_threshold__ratio = konstantin__max_threshold__ratio
    self.mode = mode

    # Generate and store thresholds
    self.rt, self.rtz = self._generate_thresholds()

Neuron populations

AlphaMN__Pool

AlphaMN__Pool(n: int | None = None, recruitment_thresholds__array: RECRUITMENT_THRESHOLDS__ARRAY | None = None, config_file: Union[str, Path, None] = None, model: str | None = None, mode: str | None = None, axon_velocities: tuple[float, float] | None = None, axon_length: float | None = None, gamma: float | None = None, cell_index: Optional[int] = None, lambda_factor: float | None = None, initial_voltage__mV: Union[float, list[float], None] = None, spike_threshold__mV: float | None = None, soma_length_range: tuple[float, float, float] | None = None, soma_diameter_range: tuple[float, float, float] | None = None, soma_capacitance_range: tuple[float, float, float] | None = None, soma_passive_conductance_range: tuple[float, float, float] | None = None, soma_passive_reversal_range: tuple[float, float, float] | None = None, soma_na3rp_conductance_range: tuple[float, float, float] | None = None, soma_naps_conductance_range: tuple[float, float, float] | None = None, soma_kdrrl_conductance_range: tuple[float, float, float] | None = None, soma_mahp_ca_conductance_range: tuple[float, float, float] | None = None, soma_mahp_k_conductance_range: tuple[float, float, float] | None = None, soma_mahp_tau_range: tuple[float, float, float] | None = None, soma_gh_conductance_range: tuple[float, float, float] | None = None, dendrite_length_range: tuple[float, float, float] | None = None, dendrite_diameter_range: tuple[float, float, float] | None = None, dendrite_passive_conductance_range: tuple[float, float, float] | None = None, dendrite_passive_reversal_range: tuple[float, float, float] | None = None, dendrite_resistance_range: tuple[float, float, float] | None = None, dendrite_capacitance_range: tuple[float, float, float] | None = None, dendrite_gh_conductance_range: tuple[float, float, float] | None = None, dendrite_ca_conductance_ranges: tuple[tuple[float, float, float], ...] | None = None, dendrite_ca_theta_m_range: tuple[float, float, float] | None = None, dendrite_ca_theta_h_range: tuple[float, float, float] | None = None)

Bases: _Pool

Container for a population of alpha motor neurons.

Manages a collection of AlphaMN (alpha motor neuron) cells with different biophysical models: ModALS or Powers2017. These cells form the final common pathway for motor control.

Parameters:

Name Type Description Default
n int

Number of alpha motor neurons to create.

None
recruitment_thresholds__array RECRUITMENT_THRESHOLDS__ARRAY

Array of recruitment thresholds for each motor neuron, by default None.

None
config_file str or Path

Path to YAML configuration file containing model parameters. If provided, parameters from this file will be used as defaults, which can be overridden by explicitly passed parameters. Can be a filename (searches in myogen/config/), relative path, or absolute path. By default uses "alpha_mn_default.yaml".

None
model str

Motor neuron model type ("NERLab" or "Powers2017"), by default "NERLab". When model="Powers2017", the soma/dendrite parameters listed under Other Parameters below become available.

None
mode str

Simulation mode ("active" or "passive"), by default "active".

None
axon_velocities tuple[float, float]

Min and max axon conduction velocities (m/s), by default (50, 65).

None
axon_length float

Length of the axon (mm), by default 0.6.

None
gamma float

Neuromodulation level (a.u.), by default 0.2.

None
cell_index Optional[int]

Specific cell index to create (creates only one cell), by default None.

None
lambda_factor float

Lambda factor for Powers2017 model persistent sodium scaling, by default 1.0.

None
initial_voltage__mV float or list[float]

Initial membrane voltage (mV), by default -67.

None
spike_threshold__mV float

Spike detection threshold for recording motor neuron spikes, by default 50.0. Motor neurons have large action potentials (80-100 mV) requiring higher thresholds.

None

Other Parameters:

Name Type Description
soma_length_range tuple[float, float, float]

Soma length [min, max, curve] (um). Only used when model="Powers2017".

soma_diameter_range tuple[float, float, float]

Soma diameter [min, max, curve] (um).

soma_capacitance_range tuple[float, float, float]

Soma capacitance [min, max, curve] (uF/cm²).

soma_passive_conductance_range tuple[float, float, float]

Soma passive conductance [min, max, curve] (S/cm²).

soma_passive_reversal_range tuple[float, float, float]

Soma passive reversal potential [min, max, curve] (mV).

soma_na3rp_conductance_range tuple[float, float, float]

Soma Na3RP conductance [min, max, curve] (S/cm²).

soma_naps_conductance_range tuple[float, float, float]

Soma NaPS conductance [min, max, curve] (S/cm²).

soma_kdrrl_conductance_range tuple[float, float, float]

Soma KDRRL conductance [min, max, curve] (S/cm²).

soma_mahp_ca_conductance_range tuple[float, float, float]

Soma mAHP calcium conductance [min, max, curve] (S/cm²).

soma_mahp_k_conductance_range tuple[float, float, float]

Soma mAHP potassium conductance [min, max, curve] (S/cm²).

soma_mahp_tau_range tuple[float, float, float]

Soma mAHP time constant [min, max, curve] (ms).

soma_gh_conductance_range tuple[float, float, float]

Soma h-current conductance [min, max, curve] (S/cm²).

dendrite_length_range tuple[float, float, float]

Dendrite length [min, max, curve] (um).

dendrite_diameter_range tuple[float, float, float]

Dendrite diameter [min, max, curve] (um).

dendrite_passive_conductance_range tuple[float, float, float]

Dendrite passive conductance [min, max, curve] (S/cm²).

dendrite_passive_reversal_range tuple[float, float, float]

Dendrite passive reversal potential [min, max, curve] (mV).

dendrite_resistance_range tuple[float, float, float]

Dendrite axial resistance [min, max, curve] (Ω·cm).

dendrite_capacitance_range tuple[float, float, float]

Dendrite capacitance [min, max, curve] (uF/cm²).

dendrite_gh_conductance_range tuple[float, float, float]

Dendrite h-current conductance [min, max, curve] (S/cm²).

dendrite_ca_conductance_ranges tuple[tuple[float, float, float], ...]

L-type Ca conductance ranges for each dendrite (4 tuples).

dendrite_ca_theta_m_range tuple[float, float, float]

Ca channel activation threshold [min, max, curve] (mV).

dendrite_ca_theta_h_range tuple[float, float, float]

Ca channel inactivation threshold [min, max, curve] (mV).

Source code in myogen/simulator/neuron/populations/motor_neurons.py
def __init__(
    self,
    n: int | None = None,
    recruitment_thresholds__array: RECRUITMENT_THRESHOLDS__ARRAY | None = None,
    config_file: Union[str, Path, None] = None,
    model: str | None = None,
    mode: str | None = None,
    axon_velocities: tuple[float, float] | None = None,
    axon_length: float | None = None,
    gamma: float | None = None,
    cell_index: Optional[int] = None,
    lambda_factor: float | None = None,
    initial_voltage__mV: Union[float, list[float], None] = None,
    spike_threshold__mV: float | None = None,
    # Powers2017 parameters
    # Soma parameters
    soma_length_range: tuple[float, float, float] | None = None,
    soma_diameter_range: tuple[float, float, float] | None = None,
    soma_capacitance_range: tuple[float, float, float] | None = None,
    soma_passive_conductance_range: tuple[float, float, float] | None = None,
    soma_passive_reversal_range: tuple[float, float, float] | None = None,
    soma_na3rp_conductance_range: tuple[float, float, float] | None = None,
    soma_naps_conductance_range: tuple[float, float, float] | None = None,
    soma_kdrrl_conductance_range: tuple[float, float, float] | None = None,
    soma_mahp_ca_conductance_range: tuple[float, float, float] | None = None,
    soma_mahp_k_conductance_range: tuple[float, float, float] | None = None,
    soma_mahp_tau_range: tuple[float, float, float] | None = None,
    soma_gh_conductance_range: tuple[float, float, float] | None = None,
    # Dendrite parameters
    dendrite_length_range: tuple[float, float, float] | None = None,
    dendrite_diameter_range: tuple[float, float, float] | None = None,
    dendrite_passive_conductance_range: tuple[float, float, float] | None = None,
    dendrite_passive_reversal_range: tuple[float, float, float] | None = None,
    dendrite_resistance_range: tuple[float, float, float] | None = None,
    dendrite_capacitance_range: tuple[float, float, float] | None = None,
    dendrite_gh_conductance_range: tuple[float, float, float] | None = None,
    dendrite_ca_conductance_ranges: tuple[tuple[float, float, float], ...] | None = None,
    dendrite_ca_theta_m_range: tuple[float, float, float] | None = None,
    dendrite_ca_theta_h_range: tuple[float, float, float] | None = None,
):
    # Load configuration from YAML files with fallback mechanism
    if config_file is None:
        config_file = "alpha_mn_default.yaml"

    try:
        # Always load the default configuration first
        default_config = load_yaml_config("alpha_mn_default.yaml")

        if config_file == "alpha_mn_default.yaml":
            config = default_config
        else:
            # Load the specific config file and merge with defaults
            specific_config = load_yaml_config(config_file)
            config = merge_configs(default_config, specific_config)

    except ImportError:
        # PyYAML not installed
        raise ImportError(
            f"PyYAML is required to load config file '{config_file}'. "
            "Install it with: pip install pyyaml"
        )

    # Helper function to get parameter value (explicit param > config > required)
    def get_param(param_value, config_key):
        if param_value is not None:
            return param_value
        if config_key not in config:
            # Debug information
            available_keys = (
                list(config.keys()) if isinstance(config, dict) else "config is not a dict"
            )
            raise ValueError(
                f"Parameter '{config_key}' not found in merged config and not provided explicitly. "
                f"Available top-level keys: {available_keys}. "
                f"Config file used: '{config_file}'"
            )
        return config[config_key]

    def get_nested_param(param_value, *config_keys):
        """Get parameter from nested config structure."""
        if param_value is not None:
            return param_value
        value = config
        for key in config_keys:
            if isinstance(value, dict):
                value = value.get(key, {})
            else:
                raise ValueError(
                    f"Parameter path {' -> '.join(config_keys)} not found in merged config and not provided explicitly. "
                    f"Failed at key '{key}'. Config file used: '{config_file}'"
                )
        if value == {} or value is None:
            raise ValueError(
                f"Parameter path {' -> '.join(config_keys)} not found in merged config and not provided explicitly. "
                f"Config file used: '{config_file}'"
            )
        return value

    # Set basic parameters
    self.n = n
    self.recruitment_thresholds__array = recruitment_thresholds__array

    if self.recruitment_thresholds__array is not None:
        self.n = len(self.recruitment_thresholds__array)

    if self.n is None and self.recruitment_thresholds__array is None:
        raise ValueError("Either n or recruitment_thresholds__array must be provided.")

    self.model = get_param(model, "model")
    self.mode = get_param(mode, "mode")
    self.axon_velocities = get_param(axon_velocities, "axon_velocities")
    self.axon_length = get_param(axon_length, "axon_length")
    self.gamma = get_param(gamma, "gamma")
    self.cell_index = cell_index
    self.lambda_factor = get_param(lambda_factor, "lambda_factor")

    # Store Powers2017 parameters
    self.soma_length_range = get_nested_param(
        soma_length_range, "powers2017", "soma", "length_range"
    )
    self.soma_diameter_range = get_nested_param(
        soma_diameter_range, "powers2017", "soma", "diameter_range"
    )
    self.soma_capacitance_range = get_nested_param(
        soma_capacitance_range, "powers2017", "soma", "capacitance_range"
    )
    self.soma_passive_conductance_range = get_nested_param(
        soma_passive_conductance_range, "powers2017", "soma", "passive_conductance_range"
    )
    self.soma_passive_reversal_range = get_nested_param(
        soma_passive_reversal_range, "powers2017", "soma", "passive_reversal_range"
    )
    self.soma_na3rp_conductance_range = get_nested_param(
        soma_na3rp_conductance_range, "powers2017", "soma", "na3rp_conductance_range"
    )
    self.soma_naps_conductance_range = get_nested_param(
        soma_naps_conductance_range, "powers2017", "soma", "naps_conductance_range"
    )
    self.soma_kdrrl_conductance_range = get_nested_param(
        soma_kdrrl_conductance_range, "powers2017", "soma", "kdrrl_conductance_range"
    )
    self.soma_mahp_ca_conductance_range = get_nested_param(
        soma_mahp_ca_conductance_range, "powers2017", "soma", "mahp_ca_conductance_range"
    )
    self.soma_mahp_k_conductance_range = get_nested_param(
        soma_mahp_k_conductance_range, "powers2017", "soma", "mahp_k_conductance_range"
    )
    self.soma_mahp_tau_range = get_nested_param(
        soma_mahp_tau_range, "powers2017", "soma", "mahp_tau_range"
    )
    self.soma_gh_conductance_range = get_nested_param(
        soma_gh_conductance_range, "powers2017", "soma", "gh_conductance_range"
    )
    self.dendrite_length_range = get_nested_param(
        dendrite_length_range, "powers2017", "dendrite", "length_range"
    )
    self.dendrite_diameter_range = get_nested_param(
        dendrite_diameter_range, "powers2017", "dendrite", "diameter_range"
    )
    self.dendrite_passive_conductance_range = get_nested_param(
        dendrite_passive_conductance_range,
        "powers2017",
        "dendrite",
        "passive_conductance_range",
    )
    self.dendrite_passive_reversal_range = get_nested_param(
        dendrite_passive_reversal_range, "powers2017", "dendrite", "passive_reversal_range"
    )
    self.dendrite_resistance_range = get_nested_param(
        dendrite_resistance_range, "powers2017", "dendrite", "resistance_range"
    )
    self.dendrite_capacitance_range = get_nested_param(
        dendrite_capacitance_range, "powers2017", "dendrite", "capacitance_range"
    )
    self.dendrite_gh_conductance_range = get_nested_param(
        dendrite_gh_conductance_range, "powers2017", "dendrite", "gh_conductance_range"
    )
    self.dendrite_ca_conductance_ranges = get_nested_param(
        dendrite_ca_conductance_ranges, "powers2017", "dendrite", "ca_conductance_ranges"
    )
    self.dendrite_ca_theta_m_range = get_nested_param(
        dendrite_ca_theta_m_range, "powers2017", "dendrite", "ca_theta_m_range"
    )
    self.dendrite_ca_theta_h_range = get_nested_param(
        dendrite_ca_theta_h_range, "powers2017", "dendrite", "ca_theta_h_range"
    )

    # Store NERLab napp parameters
    self.napp_m_alpha_A = get_nested_param(None, "nerlab", "napp", "m_alpha_A")
    self.napp_m_alpha_v_offset = get_nested_param(None, "nerlab", "napp", "m_alpha_v_offset")
    self.napp_m_alpha_k = get_nested_param(None, "nerlab", "napp", "m_alpha_k")
    self.napp_m_beta_A = get_nested_param(None, "nerlab", "napp", "m_beta_A")
    self.napp_m_beta_v_offset = get_nested_param(None, "nerlab", "napp", "m_beta_v_offset")
    self.napp_m_beta_k = get_nested_param(None, "nerlab", "napp", "m_beta_k")

    self.napp_h_alpha_A = get_nested_param(None, "nerlab", "napp", "h_alpha_A")
    self.napp_h_alpha_v_offset = get_nested_param(None, "nerlab", "napp", "h_alpha_v_offset")
    self.napp_h_alpha_tau = get_nested_param(None, "nerlab", "napp", "h_alpha_tau")
    self.napp_h_beta_A = get_nested_param(None, "nerlab", "napp", "h_beta_A")
    self.napp_h_beta_v_offset = get_nested_param(None, "nerlab", "napp", "h_beta_v_offset")
    self.napp_h_beta_k = get_nested_param(None, "nerlab", "napp", "h_beta_k")

    self.napp_p_alpha_A = get_nested_param(None, "nerlab", "napp", "p_alpha_A")
    self.napp_p_alpha_v_offset = get_nested_param(None, "nerlab", "napp", "p_alpha_v_offset")
    self.napp_p_alpha_k = get_nested_param(None, "nerlab", "napp", "p_alpha_k")
    self.napp_p_beta_A = get_nested_param(None, "nerlab", "napp", "p_beta_A")
    self.napp_p_beta_v_offset = get_nested_param(None, "nerlab", "napp", "p_beta_v_offset")
    self.napp_p_beta_k = get_nested_param(None, "nerlab", "napp", "p_beta_k")

    self.napp_n_alpha_A = get_nested_param(None, "nerlab", "napp", "n_alpha_A")
    self.napp_n_alpha_v_offset = get_nested_param(None, "nerlab", "napp", "n_alpha_v_offset")
    self.napp_n_alpha_k = get_nested_param(None, "nerlab", "napp", "n_alpha_k")
    self.napp_n_beta_A = get_nested_param(None, "nerlab", "napp", "n_beta_A")
    self.napp_n_beta_v_offset = get_nested_param(None, "nerlab", "napp", "n_beta_v_offset")
    self.napp_n_beta_tau = get_nested_param(None, "nerlab", "napp", "n_beta_tau")

    self.napp_r_alpha_A = get_nested_param(None, "nerlab", "napp", "r_alpha_A")
    self.napp_r_alpha_v_offset = get_nested_param(None, "nerlab", "napp", "r_alpha_v_offset")
    self.napp_r_alpha_k = get_nested_param(None, "nerlab", "napp", "r_alpha_k")

    # Store NERLab soma parameters
    self.nerlab_soma_diameter_range = get_nested_param(None, "nerlab", "soma", "diameter_range")
    self.nerlab_soma_gnabar_range = get_nested_param(None, "nerlab", "soma", "gnabar_range")
    self.nerlab_soma_gnapbar_range = get_nested_param(None, "nerlab", "soma", "gnapbar_range")
    self.nerlab_soma_gkfbar_range = get_nested_param(None, "nerlab", "soma", "gkfbar_range")
    self.nerlab_soma_gksbar_range = get_nested_param(None, "nerlab", "soma", "gksbar_range")
    self.nerlab_soma_mact_range = get_nested_param(None, "nerlab", "soma", "mact_range")
    self.nerlab_soma_rinact_range = get_nested_param(None, "nerlab", "soma", "rinact_range")
    self.nerlab_soma_gls_range = get_nested_param(None, "nerlab", "soma", "gls_range")
    self.nerlab_soma_ena = get_nested_param(None, "nerlab", "soma", "ena")
    self.nerlab_soma_ek = get_nested_param(None, "nerlab", "soma", "ek")
    self.nerlab_soma_el_napp = get_nested_param(None, "nerlab", "soma", "el_napp")
    self.nerlab_soma_vtraub_napp = get_nested_param(None, "nerlab", "soma", "vtraub_napp")
    self.nerlab_soma_Ra = get_nested_param(None, "nerlab", "soma", "Ra")
    self.nerlab_soma_cm = get_nested_param(None, "nerlab", "soma", "cm")

    # Store NERLab dendrite parameters
    self.nerlab_dendrite_diameter_range = get_nested_param(
        None, "nerlab", "dendrite", "diameter_range"
    )
    self.nerlab_dendrite_length_range = get_nested_param(
        None, "nerlab", "dendrite", "length_range"
    )
    self.nerlab_dendrite_gcaLbar_range = get_nested_param(
        None, "nerlab", "dendrite", "gcaLbar_range"
    )
    self.nerlab_dendrite_vtraub_caL_range = get_nested_param(
        None, "nerlab", "dendrite", "vtraub_caL_range"
    )
    self.nerlab_dendrite_ltau_caL_range = get_nested_param(
        None, "nerlab", "dendrite", "ltau_caL_range"
    )
    self.nerlab_dendrite_gl_caL_range = get_nested_param(
        None, "nerlab", "dendrite", "gl_caL_range"
    )
    self.nerlab_dendrite_Ra = get_nested_param(None, "nerlab", "dendrite", "Ra")
    self.nerlab_dendrite_cm = get_nested_param(None, "nerlab", "dendrite", "cm")
    self.nerlab_dendrite_ecaL = get_nested_param(None, "nerlab", "dendrite", "ecaL")
    self.nerlab_dendrite_el_caL = get_nested_param(None, "nerlab", "dendrite", "el_caL")

    if self.model == "NERLab":
        _cells = self._create_nerlab_cells()
    elif self.model == "Powers2017":
        _cells = self._create_powers2017_cells()
    else:
        raise ValueError("Could not find the specific model for alpha MNs.")

    # Get initial voltage and spike threshold from config if not provided
    _initial_voltage = get_param(initial_voltage__mV, "initial_voltage__mV")
    _spike_threshold = get_param(spike_threshold__mV, "spike_threshold__mV")

    super().__init__(
        cells=_cells,
        initial_voltage__mV=_initial_voltage,
        spike_threshold__mV=_spike_threshold,
    )

DescendingDrive__Pool

DescendingDrive__Pool(n: int, timestep__ms: Quantity__ms | None = None, process_type: str = 'poisson', shape: float = 3.0)

Bases: _Pool

Container for a population of descending drive neurons.

Manages a collection of DD cells that generate spike trains using either Poisson or Gamma point processes for cortical input to spinal circuits.

Parameters:

Name Type Description Default
n int

Number of descending drive neurons to create.

required
timestep__ms Quantity__ms

Time step for simulation as a Quantity with units of milliseconds (required).

None
process_type str

Type of point process: "poisson" or "gamma", by default "poisson". - "poisson": discrete-time Poisson process, irregular firing (CV=1.0, exact as dt->0) - "gamma": More regular firing with CV controlled by shape parameter

'poisson'
shape float

Shape parameter for Gamma process (only used when process_type="gamma"), by default 3.0. Controls spike regularity: - shape=1: Poisson-like (CV=1.0) - shape=2-5: Typical cortical neuron regularity (CV=0.45-0.71) - Higher values: More regular firing (CV=1/sqrt(shape))

3.0
Source code in myogen/simulator/neuron/populations/descending_drive.py
def __init__(
    self,
    n: int,
    timestep__ms: Quantity__ms | None = None,
    process_type: str = "poisson",
    shape: float = 3.0,
):
    if timestep__ms is None:
        raise ValueError("timestep__ms is required")

    self.n = n
    self.timestep__ms = timestep__ms
    self.process_type = process_type
    self.shape = shape

    if process_type.lower() == "gamma":
        _cells = [
            cells.DD_Gamma(
                timestep__ms=timestep__ms,
                shape=shape,
                pool__ID=i,
            )
            for i in range(n)
        ]
    elif process_type.lower() == "poisson":
        _cells = [cells.DD(dt=timestep__ms, pool__ID=i) for i in range(n)]
    else:
        raise ValueError(
            f"Invalid process_type '{process_type}'. Must be 'poisson' or 'gamma'."
        )

    super().__init__(cells=_cells)

DescendingDrive_Gamma__Pool is a backward-compatibility alias. New code should prefer DescendingDrive__Pool(process_type="gamma", shape=...).

DescendingDrive_Gamma__Pool

DescendingDrive_Gamma__Pool(n: int, timestep__ms: Quantity__ms, shape: float = 3.0)

Bases: _Pool

Container for a population of descending drive neurons using Gamma process.

Manages a collection of DD_Gamma cells that generate Gamma-distributed spike trains for more regular cortical input to spinal circuits, typical of cortical neuron firing patterns.

Note: This class is kept for backward compatibility. Consider using DescendingDrive__Pool with process_type='gamma' instead.

Parameters:

Name Type Description Default
n int

Number of descending drive neurons to create.

required
timestep__ms Quantity__ms

Time step for simulation as a Quantity with units of milliseconds.

required
shape float

Shape parameter controlling spike regularity, by default 3.0. - shape=1: Poisson-like (irregular) firing - shape=2-5: Typical cortical neuron regularity - Higher values: More regular, clock-like firing

3.0
Source code in myogen/simulator/neuron/populations/descending_drive.py
def __init__(
    self,
    n: int,
    timestep__ms: Quantity__ms,
    shape: float = 3.0,
):
    self.n = n
    self.timestep__ms = timestep__ms
    self.shape = shape

    super().__init__(
        cells=[
            cells.DD_Gamma(
                timestep__ms=timestep__ms,
                shape=shape,
                pool__ID=i,
            )
            for i in range(n)
        ]
    )

AffIa__Pool

AffIa__Pool(n: int, timestep__ms: Quantity__ms, recruitment_thresholds: tuple[float, float] = (0, 40), axon_velocities__m_per_s: tuple[Quantity__m_per_s, Quantity__m_per_s] = (61 * m / s, 75 * m / s), axon_length__m: Quantity__m = 0.6 * m, shape: int = 145, init_order: int = 0)

Bases: _Pool

Container for a population of afferent Ia neurons.

Manages a collection of AffIa (type Ia afferent) cells that provide proprioceptive feedback from muscle spindles to spinal circuits.

Parameters:

Name Type Description Default
n int

Number of type Ia afferent neurons to create.

required
recruitment_thresholds tuple[float, float]

Min and max recruitment thresholds (Hz).

(0, 40)
axon_velocities__m_per_s tuple[Quantity__m_per_s, Quantity__m_per_s]

Min and max axon conduction velocities (m/s).

(61 * m / s, 75 * m / s)
axon_length__m Quantity__m

Length of the axon (m).

0.6 * m
shape int

Gamma shape parameter controlling ISI regularity: CV = 1/sqrt(shape). Larger values give more regular firing. Does not set the firing rate.

145
timestep__ms Quantity__ms

Time step for simulation (ms).

required
init_order int

Initial order parameter for afferent initialization.

0
Source code in myogen/simulator/neuron/populations/afferents.py
def __init__(
    self,
    n: int,
    timestep__ms: Quantity__ms,
    recruitment_thresholds: tuple[float, float] = (0, 40),
    axon_velocities__m_per_s: tuple[Quantity__m_per_s, Quantity__m_per_s] = (
        61 * pq.m / pq.s,
        75 * pq.m / pq.s,
    ),
    axon_length__m: Quantity__m = 0.6 * pq.m,
    shape: int = 145,  # Gamma shape: CV = 1/sqrt(145) = 8.3%
    init_order: int = 0,
):
    self.n = n
    self.recruitment_thresholds = recruitment_thresholds
    self.axon_velocities = axon_velocities__m_per_s
    self.axon_length = axon_length__m
    self.shape = shape
    self.timestep__ms = timestep__ms
    self.init_order = init_order

    rt = np.linspace(*recruitment_thresholds, n)
    vcon = np.linspace(*axon_velocities__m_per_s, n)

    _cells = []
    for i, (rt_i, vcon_i) in enumerate(zip(rt, vcon)):
        ia = cells.AffIa(
            RT=rt_i,
            shape=shape,
            timestep__ms=timestep__ms,
            initN=init_order,
            pool__ID=i,
        )
        ia.create_axon(length__m=axon_length__m, conduction_velocity__m_per_s=vcon_i)
        _cells.append(ia)

    super().__init__(cells=_cells)

AffII__Pool

AffII__Pool(n: int, timestep__ms: Quantity__ms, recruitment_thresholds: tuple[float, float] = (0, 40), axon_velocities__m_per_s: tuple[Quantity__m_per_s, Quantity__m_per_s] = (30 * m / s, 50 * m / s), axon_length__m: Quantity__m = 0.6 * m, shape: int = 772, init_order: int = 0)

Bases: _Pool

Container for a population of afferent II neurons.

Manages a collection of AffII (type II afferent) cells that provide secondary proprioceptive feedback from muscle spindles to spinal circuits.

Parameters:

Name Type Description Default
n int

Number of type II afferent neurons to create.

required
recruitment_thresholds tuple[float, float]

Min and max recruitment thresholds (Hz).

(0, 40)
axon_velocities__m_per_s tuple[Quantity__m_per_s, Quantity__m_per_s]

Min and max axon conduction velocities (m/s).

(30 * m / s, 50 * m / s)
axon_length__m Quantity__m

Length of the axon (m).

0.6 * m
shape int

Gamma shape parameter controlling ISI regularity: CV = 1/sqrt(shape). Larger values give more regular firing. Does not set the firing rate.

772
timestep__ms Quantity__ms

Time step for simulation (ms).

required
init_order int

Initial order parameter for afferent initialization.

0
Source code in myogen/simulator/neuron/populations/afferents.py
def __init__(
    self,
    n: int,
    timestep__ms: Quantity__ms,
    recruitment_thresholds: tuple[float, float] = (0, 40),
    axon_velocities__m_per_s: tuple[Quantity__m_per_s, Quantity__m_per_s] = (
        30 * pq.m / pq.s,
        50 * pq.m / pq.s,
    ),
    axon_length__m: Quantity__m = 0.6 * pq.m,
    shape: int = 772,  # Gamma shape: CV = 1/sqrt(772) = 3.6%
    init_order: int = 0,
):
    self.n = n
    self.recruitment_thresholds = recruitment_thresholds
    self.axon_velocities = axon_velocities__m_per_s
    self.axon_length = axon_length__m
    self.shape = shape
    self.timestep__ms = timestep__ms
    self.init_order = init_order

    rt = np.linspace(*recruitment_thresholds, n)
    vcon = np.linspace(*axon_velocities__m_per_s, n)

    _cells = []
    for i, (rt_i, vcon_i) in enumerate(zip(rt, vcon)):
        ii = cells.AffII(
            RT=rt_i,
            shape=shape,
            timestep__ms=timestep__ms,
            initN=init_order,
            pool__ID=i,
        )
        ii.create_axon(length__m=axon_length__m, conduction_velocity__m_per_s=vcon_i)
        _cells.append(ii)

    super().__init__(cells=_cells)

AffIb__Pool

AffIb__Pool(n: int, timestep__ms: Quantity__ms, recruitment_thresholds: tuple[float, float] = (0, 40), axon_velocities__m_per_s: tuple[Quantity__m_per_s, Quantity__m_per_s] = (64 * m / s, 72 * m / s), axon_length__mm: Quantity__mm = 0.6 * mm, shape: int = 145, init_order: int = 0)

Bases: _Pool

Container for a population of afferent Ib neurons.

Manages a collection of AffIb (type Ib afferent) cells that provide primary proprioceptive feedback from Golgi tendon organs to spinal circuits.

Parameters:

Name Type Description Default
n int

Number of type Ib afferent neurons to create.

required
recruitment_thresholds tuple[float, float]

Min and max recruitment thresholds (Hz).

(0, 40)
axon_velocities__m_per_s tuple[Quantity__m_per_s, Quantity__m_per_s]

Min and max axon conduction velocities (m/s).

(64 * m / s, 72 * m / s)
axon_length__mm Quantity__mm

Length of the axon (mm).

0.6 * mm
shape int

Gamma shape parameter controlling ISI regularity: CV = 1/sqrt(shape). Larger values give more regular firing. Does not set the firing rate.

145
timestep__ms Quantity__ms

Time step for simulation (ms).

required
init_order int

Initial order parameter for afferent initialization.

0
Source code in myogen/simulator/neuron/populations/afferents.py
def __init__(
    self,
    n: int,
    timestep__ms: Quantity__ms,
    recruitment_thresholds: tuple[float, float] = (0, 40),
    axon_velocities__m_per_s: tuple[Quantity__m_per_s, Quantity__m_per_s] = (
        64 * pq.m / pq.s,
        72 * pq.m / pq.s,
    ),
    axon_length__mm: Quantity__mm = 0.6 * pq.mm,
    shape: int = 145,  # Gamma shape: CV = 1/sqrt(145) = 8.3%
    init_order: int = 0,
):
    self.n = n
    self.recruitment_thresholds = recruitment_thresholds
    self.axon_velocities = axon_velocities__m_per_s
    self.axon_length = axon_length__mm
    self.shape = shape
    self.timestep__ms = timestep__ms
    self.init_order = init_order

    rt = np.linspace(*recruitment_thresholds, n)
    vcon = np.linspace(*axon_velocities__m_per_s, n)

    _cells = []
    for i, (rt_i, vcon_i) in enumerate(zip(rt, vcon)):
        ib = cells.AffIb(
            RT=rt_i,
            shape=shape,
            timestep__ms=timestep__ms,
            initN=init_order,
            pool__ID=i,
        )
        ib.create_axon(length__m=axon_length__mm.rescale(pq.m), conduction_velocity__m_per_s=vcon_i)
        _cells.append(ib)

    super().__init__(cells=_cells)

GII__Pool

GII__Pool(n: int, soma_length_range__um: tuple[float, float] = _get_interneuron_diameter_range__um(), soma_diameter_range: tuple[float, float] = _get_interneuron_diameter_range__um(), passive_conductance_range: tuple[float, float] = (3e-05, 7e-05), na3rp_conductance_range: tuple[float, float] = (0.003, 0.01), kdrrl_conductance_range: tuple[float, float] = (0.015, 0.015), mahp_ca_conductance_range: tuple[float, float] = (3e-06, 3e-06), mahp_k_conductance_range: tuple[float, float] = (0.0005, 0.0005), mahp_tau_range: tuple[float, float] = (60, 70), gh_conductance_range: tuple[float, float] = (2.5e-05, 2.5e-05), axon_velocities: tuple[float, float] = (10, 10), axon_length: float = 0.05, cell_index: Optional[int] = None, initial_voltage__mV: Union[float, list[float]] = -70.0)

Bases: _Pool

Container for a population of group II interneurons.

Manages a collection of INgII (group II interneuron) cells that provide inhibitory feedback in spinal circuits, processing type II afferent input.

Parameters:

Name Type Description Default
n int

Number of group II interneurons to create.

required
soma_length_range__um tuple[float, float]

Min and max soma length (um). By default, it is set to the estimated range for interneurons from Bui et al. 2003 [1].

_get_interneuron_diameter_range__um()
soma_diameter_range tuple[float, float]

Min and max soma diameter (um). By default, it is set to the estimated range for interneurons from Bui et al. 2003 [1].

_get_interneuron_diameter_range__um()
passive_conductance_range tuple[float, float]

Min and max passive membrane conductance (S/cm²).

(3e-05, 7e-05)
na3rp_conductance_range tuple[float, float]

Min and max Na3RP sodium channel conductance (S/cm²).

(0.003, 0.01)
kdrrl_conductance_range tuple[float, float]

Min and max KDRRL potassium channel conductance (S/cm²).

(0.015, 0.015)
mahp_ca_conductance_range tuple[float, float]

Min and max mAHP calcium conductance (S/cm²).

(3e-06, 3e-06)
mahp_k_conductance_range tuple[float, float]

Min and max mAHP potassium conductance (S/cm²).

(0.0005, 0.0005)
mahp_tau_range tuple[float, float]

Min and max mAHP time constant (ms).

(60, 70)
gh_conductance_range tuple[float, float]

Min and max h-current conductance (S/cm²).

(2.5e-05, 2.5e-05)
axon_velocities tuple[float, float]

Min and max axon conduction velocities (m/s).

(10, 10)
axon_length float

Length of the axon (mm).

0.05
cell_index int

Specific cell index to create (creates only one cell), by default None.

None
References

[1] Bui, T.V., Cushing, S., Dewey, D., Fyffe, R.E., Rose, P.K., 2003. Comparison of the Morphological and Electrotonic Properties of Renshaw Cells, Ia Inhibitory Interneurons, and Motoneurons in the Cat. Journal of Neurophysiology 90, 2900–2918. https://doi.org/10.1152/jn.00533.2003

Source code in myogen/simulator/neuron/populations/interneurons.py
def __init__(
    self,
    n: int,
    soma_length_range__um: tuple[float, float] = _get_interneuron_diameter_range__um(),
    soma_diameter_range: tuple[float, float] = _get_interneuron_diameter_range__um(),
    passive_conductance_range: tuple[float, float] = (3e-5, 7e-5),
    na3rp_conductance_range: tuple[float, float] = (0.003, 0.01),
    kdrrl_conductance_range: tuple[float, float] = (0.015, 0.015),
    mahp_ca_conductance_range: tuple[float, float] = (3e-6, 3e-6),
    mahp_k_conductance_range: tuple[float, float] = (5e-4, 5e-4),
    mahp_tau_range: tuple[float, float] = (60, 70),
    gh_conductance_range: tuple[float, float] = (2.5e-5, 2.5e-5),
    axon_velocities: tuple[float, float] = (10, 10),
    axon_length: float = 0.05,
    cell_index: Optional[int] = None,
    initial_voltage__mV: Union[float, list[float]] = -70.0,
):
    self.n = n
    self.soma_length_range__um = soma_length_range__um
    self.soma_diameter_range = soma_diameter_range
    self.passive_conductance_range = passive_conductance_range
    self.na3rp_conductance_range = na3rp_conductance_range
    self.kdrrl_conductance_range = kdrrl_conductance_range
    self.mahp_ca_conductance_range = mahp_ca_conductance_range
    self.mahp_k_conductance_range = mahp_k_conductance_range
    self.mahp_tau_range = mahp_tau_range
    self.gh_conductance_range = gh_conductance_range
    self.axon_velocities = axon_velocities
    self.axon_length = axon_length
    self.cell_index = cell_index

    sL = np.linspace(*soma_length_range__um, n)
    sdiam = np.linspace(*soma_diameter_range, n)
    sg_pas = np.linspace(*passive_conductance_range, n)
    sgbar_na3rp = np.linspace(*na3rp_conductance_range, n)
    sgMax_kdrRL = np.linspace(*kdrrl_conductance_range, n)
    sgcamax_mAHP = np.linspace(*mahp_ca_conductance_range, n)
    sgkcamax_mAHP = np.linspace(*mahp_k_conductance_range, n)
    stau_mAHP = np.linspace(*mahp_tau_range, n)
    sghbar_gh = np.linspace(*gh_conductance_range, n)
    vcon = np.linspace(*axon_velocities, n)

    if cell_index is not None:
        init, end = cell_index, cell_index + 1
    else:
        init, end = 0, n

    _cells = []
    for i, (
        sL_i,
        sdiam_i,
        sg_pas_i,
        sgbar_na3rp_i,
        sgMax_kdrRL_i,
        sgcamax_mAHP_i,
        sgkcamax_mAHP_i,
        stau_mAHP_i,
        sghbar_gh_i,
        vcon_i,
    ) in enumerate(
        zip(
            sL[init:end],
            sdiam[init:end],
            sg_pas[init:end],
            sgbar_na3rp[init:end],
            sgMax_kdrRL[init:end],
            sgcamax_mAHP[init:end],
            sgkcamax_mAHP[init:end],
            stau_mAHP[init:end],
            sghbar_gh[init:end],
            vcon[init:end],
        )
    ):
        gII = cells.INgII(pool__ID=i)

        gII.soma.L = sL_i
        gII.soma.diam = sdiam_i
        gII.soma.g_pas = sg_pas_i
        gII.soma.gbar_na3rp = sgbar_na3rp_i
        gII.soma.gMax_kdrRL = sgMax_kdrRL_i
        gII.soma.gcamax_mAHP = sgcamax_mAHP_i
        gII.soma.gkcamax_mAHP = sgkcamax_mAHP_i
        gII.soma.tau_mAHP = stau_mAHP_i
        gII.soma.ghbar_gh = sghbar_gh_i

        import quantities as pq

        gII.create_axon(
            length__m=axon_length * pq.m, conduction_velocity__m_per_s=vcon_i * pq.m / pq.s
        )
        _cells.append(gII)

    super().__init__(cells=_cells, initial_voltage__mV=initial_voltage__mV)

GIb__Pool

GIb__Pool(n: int, soma_length_range: tuple[float, float] = _get_interneuron_diameter_range__um(), soma_diameter_range: tuple[float, float] = _get_interneuron_diameter_range__um(), passive_conductance_range: tuple[float, float] = (3e-05, 8e-05), na3rp_conductance_range: tuple[float, float] = (0.01, 0.03), kdrrl_conductance_range: tuple[float, float] = (0.035, 0.028), mahp_ca_conductance_range: tuple[float, float] = (1e-06, 6e-06), mahp_k_conductance_range: tuple[float, float] = (0.0003, 0.00045), mahp_tau_range: tuple[float, float] = (120, 90), gh_conductance_range: tuple[float, float] = (2.5e-05, 2.5e-05), axon_velocities: tuple[float, float] = (10, 10), axon_length: float = 0.05, cell_index: int | None = None, initial_voltage__mV: float | list[float] = -70.0)

Bases: _Pool

Container for a population of group Ib interneurons.

Manages a collection of INgIb (group Ib interneuron) cells that provide inhibitory feedback in spinal circuits, processing type Ib afferent input from Golgi tendon organs.

Parameters:

Name Type Description Default
n int

Number of group Ib interneurons to create.

required
soma_length_range tuple[float, float]

Min and max soma length (um).

_get_interneuron_diameter_range__um()
soma_diameter_range tuple[float, float]

Min and max soma diameter (um).

_get_interneuron_diameter_range__um()
passive_conductance_range tuple[float, float]

Min and max passive membrane conductance (S/cm²).

(3e-05, 8e-05)
na3rp_conductance_range tuple[float, float]

Min and max Na3RP sodium channel conductance (S/cm²).

(0.01, 0.03)
kdrrl_conductance_range tuple[float, float]

Min and max KDRRL potassium channel conductance (S/cm²).

(0.035, 0.028)
mahp_ca_conductance_range tuple[float, float]

Min and max mAHP calcium conductance (S/cm²).

(1e-06, 6e-06)
mahp_k_conductance_range tuple[float, float]

Min and max mAHP potassium conductance (S/cm²).

(0.0003, 0.00045)
mahp_tau_range tuple[float, float]

Min and max mAHP time constant (ms).

(120, 90)
gh_conductance_range tuple[float, float]

Min and max h-current conductance (S/cm²).

(2.5e-05, 2.5e-05)
axon_velocities tuple[float, float]

Min and max axon conduction velocities (m/s).

(10, 10)
axon_length float

Length of the axon (mm).

0.05
cell_index Optional[int]

Specific cell index to create (creates only one cell), by default None.

None
Source code in myogen/simulator/neuron/populations/interneurons.py
def __init__(
    self,
    n: int,
    soma_length_range: tuple[float, float] = _get_interneuron_diameter_range__um(),
    soma_diameter_range: tuple[float, float] = _get_interneuron_diameter_range__um(),
    passive_conductance_range: tuple[float, float] = (3e-5, 8e-5),
    na3rp_conductance_range: tuple[float, float] = (0.01, 0.03),
    kdrrl_conductance_range: tuple[float, float] = (0.035, 0.028),
    mahp_ca_conductance_range: tuple[float, float] = (1e-6, 6e-6),
    mahp_k_conductance_range: tuple[float, float] = (3e-4, 4.5e-4),
    mahp_tau_range: tuple[float, float] = (120, 90),
    gh_conductance_range: tuple[float, float] = (2.5e-5, 2.5e-5),
    axon_velocities: tuple[float, float] = (10, 10),
    axon_length: float = 0.05,
    cell_index: int | None = None,
    initial_voltage__mV: float | list[float] = -70.0,
):
    self.n = n
    self.soma_length_range = soma_length_range
    self.soma_diameter_range = soma_diameter_range
    self.passive_conductance_range = passive_conductance_range
    self.na3rp_conductance_range = na3rp_conductance_range
    self.kdrrl_conductance_range = kdrrl_conductance_range
    self.mahp_ca_conductance_range = mahp_ca_conductance_range
    self.mahp_k_conductance_range = mahp_k_conductance_range
    self.mahp_tau_range = mahp_tau_range
    self.gh_conductance_range = gh_conductance_range
    self.axon_velocities = axon_velocities
    self.axon_length = axon_length
    self.cell_index = cell_index

    sL = np.linspace(*soma_length_range, n)
    sdiam = np.linspace(*soma_diameter_range, n)
    sg_pas = np.linspace(*passive_conductance_range, n)
    sgbar_na3rp = np.linspace(*na3rp_conductance_range, n)
    sgMax_kdrRL = np.linspace(*kdrrl_conductance_range, n)
    sgcamax_mAHP = np.linspace(*mahp_ca_conductance_range, n)
    sgkcamax_mAHP = np.linspace(*mahp_k_conductance_range, n)
    stau_mAHP = np.linspace(*mahp_tau_range, n)
    sghbar_gh = np.linspace(*gh_conductance_range, n)
    vcon = np.linspace(*axon_velocities, n)

    if cell_index is not None:
        init, end = cell_index, cell_index + 1
    else:
        init, end = 0, n

    _cells = []
    for i, (
        sL_i,
        sdiam_i,
        sg_pas_i,
        sgbar_na3rp_i,
        sgMax_kdrRL_i,
        sgcamax_mAHP_i,
        sgkcamax_mAHP_i,
        stau_mAHP_i,
        sghbar_gh_i,
        vcon_i,
    ) in enumerate(
        zip(
            sL[init:end],
            sdiam[init:end],
            sg_pas[init:end],
            sgbar_na3rp[init:end],
            sgMax_kdrRL[init:end],
            sgcamax_mAHP[init:end],
            sgkcamax_mAHP[init:end],
            stau_mAHP[init:end],
            sghbar_gh[init:end],
            vcon[init:end],
        )
    ):
        gIb = cells.INgIb(pool__ID=i)

        gIb.soma.L = sL_i
        gIb.soma.diam = sdiam_i
        gIb.soma.g_pas = sg_pas_i
        gIb.soma.gbar_na3rp = sgbar_na3rp_i
        gIb.soma.gMax_kdrRL = sgMax_kdrRL_i
        gIb.soma.gcamax_mAHP = sgcamax_mAHP_i
        gIb.soma.gkcamax_mAHP = sgkcamax_mAHP_i
        gIb.soma.tau_mAHP = stau_mAHP_i
        gIb.soma.ghbar_gh = sghbar_gh_i

        import quantities as pq

        gIb.create_axon(
            length__m=axon_length * pq.m, conduction_velocity__m_per_s=vcon_i * pq.m / pq.s
        )

        _cells.append(gIb)

    super().__init__(cells=_cells, initial_voltage__mV=initial_voltage__mV)

Network & runner

Network

Network(populations: dict[str, _Pool], spike_recording: Optional[dict] = None)

Modern neural network builder with intuitive connection API.

Provides a clean, discoverable interface for creating neural network connections while maintaining compatibility with existing NEURON-based infrastructure.

Initialize network with neural populations.

Parameters:

Name Type Description Default
populations dict[str, Union[list, Any]]

Dictionary mapping population names to Pool objects or lists of neuron objects. Pool objects will have .neurons extracted, lists used directly. Example: {"alpha_mn": alphaMN_pool, "ia": ia_pool}

required
spike_recording dict

Dictionary containing 'idvec' and 'spkvec' for spike recording. Example: {"idvec": {"aMN": h.Vector()}, "spkvec": {"aMN": h.Vector()}}

None
Source code in myogen/simulator/neuron/network.py
def __init__(self, populations: dict[str, _Pool], spike_recording: Optional[dict] = None):
    """
    Initialize network with neural populations.

    Parameters
    ----------
    populations : dict[str, Union[list, Any]]
        Dictionary mapping population names to Pool objects or lists of neuron objects.
        Pool objects will have .neurons extracted, lists used directly.
        Example: {"alpha_mn": alphaMN_pool, "ia": ia_pool}
    spike_recording : dict, optional
        Dictionary containing 'idvec' and 'spkvec' for spike recording.
        Example: {"idvec": {"aMN": h.Vector()}, "spkvec": {"aMN": h.Vector()}}
    """
    self.populations: dict[str, _Pool] = populations
    self.connections = []
    self._netcons_by_connection: dict[tuple[str, str], list[h.NetCon]] = {}
    self.spike_recording = spike_recording

setup_spike_recording

setup_spike_recording()

Set up spike recording NetCons for all neurons in all populations.

This creates additional NetCons specifically for recording spikes from neurons that might not have outgoing connections but still need spike recording for analysis.

Source code in myogen/simulator/neuron/network.py
def setup_spike_recording(self):
    """
    Set up spike recording NetCons for all neurons in all populations.

    This creates additional NetCons specifically for recording spikes from neurons
    that might not have outgoing connections but still need spike recording for analysis.
    """
    if not self.spike_recording:
        return

    for pop_name, population in self.populations.items():
        # Skip non-neuron populations (like gMN which is a config dict)
        if not hasattr(population, "__iter__") or isinstance(population, dict):
            continue

        # Get spike recording vectors for this population
        id_vector = self.spike_recording.get("idvec", {}).get(pop_name)
        spike_vector = self.spike_recording.get("spkvec", {}).get(pop_name)

        if id_vector is not None and spike_vector is not None:
            # Create NetCons for spike recording (no target, just recording)
            recording_netcons = []
            for neuron in population:
                # Skip if this is not actually a neuron object
                if not hasattr(neuron, "soma") and not hasattr(neuron, "ns"):
                    continue

                # Create a NetCon from the neuron to None (just for recording)
                if hasattr(neuron, "soma"):
                    # Compartmental neuron
                    nc = h.NetCon(neuron.soma(0.5)._ref_v, None, sec=neuron.soma)
                else:
                    # Point process neuron
                    nc = h.NetCon(neuron.ns, None)

                # Set up spike recording with population-specific threshold
                if hasattr(population, "spike_threshold__mV"):
                    nc.threshold = _to_float(population.spike_threshold__mV, pq.mV)
                else:
                    nc.threshold = _to_float(DEFAULT_SPIKE_THRESHOLD, pq.mV)  # Fallback for populations without explicit threshold
                nc.record(spike_vector, id_vector, neuron.global__ID)
                recording_netcons.append(nc)

            # Store these recording NetCons separately
            connection_key = (pop_name, "spike_recording")
            self._netcons_by_connection[connection_key] = recording_netcons

connect

connect(source: str, target: str, probability: float = 1.0, weight__uS: Quantity__uS = DEFAULT_SYNAPTIC_WEIGHT, delay__ms: Quantity__ms = DEFAULT_SYNAPTIC_DELAY, threshold__mV: Quantity__mV = DEFAULT_SPIKE_THRESHOLD, deterministic: bool = False, inhibitory: bool = False) -> list

Connect two neural populations with specified parameters.

Parameters:

Name Type Description Default
source str

Name of source population (must exist in populations dict).

required
target str

Name of target population (must exist in populations dict).

required
probability float

Connection probability between 0.0 and 1.0, by default 1.0. Each source-target neuron pair connects with this probability (if deterministic=False). If deterministic=True, each source connects to exactly int(probability × n_targets) targets.

1.0
weight__uS float

Synaptic weight in microsiemens, by default 0.6.

DEFAULT_SYNAPTIC_WEIGHT
delay__ms float

Synaptic delay in milliseconds, by default 1.0.

DEFAULT_SYNAPTIC_DELAY
threshold__mV float

Spike threshold in millivolts, by default -10.0.

DEFAULT_SPIKE_THRESHOLD
deterministic bool

If True, each source neuron connects to exactly int(probability × n_targets) randomly selected target neurons. If False, uses probabilistic sampling. Default False.

False
inhibitory bool

If True, connect to inhibitory synapses on target neurons (reversal < -40 mV). If False, connect to excitatory synapses (reversal >= -40 mV). Default False. Use inhibitory=True for connections from inhibitory interneurons (e.g., gII→aMN, gIb→aMN).

False

Returns:

Type Description
list[NetCon]

List of created NEURON NetCon objects for this connection group.

Raises:

Type Description
ValueError

If source or target populations don't exist, or probability out of range.

Source code in myogen/simulator/neuron/network.py
def connect(
    self,
    source: str,
    target: str,
    probability: float = 1.0,
    weight__uS: Quantity__uS = DEFAULT_SYNAPTIC_WEIGHT,
    delay__ms: Quantity__ms = DEFAULT_SYNAPTIC_DELAY,
    threshold__mV: Quantity__mV = DEFAULT_SPIKE_THRESHOLD,
    deterministic: bool = False,
    inhibitory: bool = False,
) -> list:
    """
    Connect two neural populations with specified parameters.

    Parameters
    ----------
    source : str
        Name of source population (must exist in populations dict).
    target : str
        Name of target population (must exist in populations dict).
    probability : float, optional
        Connection probability between 0.0 and 1.0, by default 1.0.
        Each source-target neuron pair connects with this probability (if deterministic=False).
        If deterministic=True, each source connects to exactly int(probability × n_targets) targets.
    weight__uS : float, optional
        Synaptic weight in microsiemens, by default 0.6.
    delay__ms : float, optional
        Synaptic delay in milliseconds, by default 1.0.
    threshold__mV : float, optional
        Spike threshold in millivolts, by default -10.0.
    deterministic : bool, optional
        If True, each source neuron connects to exactly int(probability × n_targets) randomly
        selected target neurons. If False, uses probabilistic sampling. Default False.
    inhibitory : bool, optional
        If True, connect to inhibitory synapses on target neurons (reversal < -40 mV).
        If False, connect to excitatory synapses (reversal >= -40 mV). Default False.
        Use inhibitory=True for connections from inhibitory interneurons (e.g., gII→aMN, gIb→aMN).

    Returns
    -------
    list[h.NetCon]
        List of created NEURON NetCon objects for this connection group.

    Raises
    ------
    ValueError
        If source or target populations don't exist, or probability out of range.
    """
    # Validation
    if source not in self.populations:
        raise ValueError(f"Source population '{source}' not found")
    if target not in self.populations:
        raise ValueError(f"Target population '{target}' not found")
    if not 0.0 <= probability <= 1.0:
        raise ValueError(f"Probability must be 0.0-1.0, got {probability}")

    # Rescale Quantities to NEURON's native units (or pass plain floats through)
    weight_value = _to_float(weight__uS, pq.uS)
    threshold_value = _to_float(threshold__mV, pq.mV)
    delay_value = _to_float(delay__ms, pq.ms)

    # Extract spike recording vectors for source population
    id_vector = None
    spike_vector = None
    if self.spike_recording:
        id_vector = self.spike_recording.get("idvec", {}).get(source)
        spike_vector = self.spike_recording.get("spkvec", {}).get(source)

    # Create connections using existing infrastructure
    netcons = _connect_populations(
        populations=self.populations,
        source_pop=source,
        target_pop=target,
        connection_probability=probability,
        synaptic_delay=delay_value,
        synaptic_weight=weight_value,
        spike_threshold=threshold_value,
        id_vector=id_vector,
        spike_vector=spike_vector,
        deterministic=deterministic,
        inhibitory=inhibitory,
    )

    self.connections.append(
        {
            "type": "neural",
            "source": source,
            "target": target,
            "probability": probability,
            "weight__uS": weight_value,
            "delay__ms": delay_value,
            "threshold__mV": threshold_value,
            "inhibitory": inhibitory,
        }
    )
    self._netcons_by_connection[(source, target)] = netcons

    return netcons

connect_to_muscle

connect_to_muscle(source: str, muscle, activation_callback: Callable, weight__uS: Quantity__uS = DEFAULT_SYNAPTIC_WEIGHT, delay__ms: Quantity__ms = DEFAULT_SYNAPTIC_DELAY, threshold__mV: Quantity__mV = DEFAULT_SPIKE_THRESHOLD) -> list

Connect a neural population to a muscle with activation callback.

Parameters:

Name Type Description Default
source str

Name of motor neuron population.

required
muscle object

Muscle object for force generation.

required
activation_callback Callable

Function called when motor neurons fire. Expected signature: callback(neuron_id, muscle, delay_time)

required
weight__uS float

Synaptic weight in microsiemens, by default 1.0.

DEFAULT_SYNAPTIC_WEIGHT
threshold__mV float

Spike threshold in millivolts, by default -10.0.

DEFAULT_SPIKE_THRESHOLD

Returns:

Type Description
list[NetCon]

List of motor neuron to muscle NetCon objects.

Source code in myogen/simulator/neuron/network.py
def connect_to_muscle(
    self,
    source: str,
    muscle,
    activation_callback: Callable,
    weight__uS: Quantity__uS = DEFAULT_SYNAPTIC_WEIGHT,
    delay__ms: Quantity__ms = DEFAULT_SYNAPTIC_DELAY,
    threshold__mV: Quantity__mV = DEFAULT_SPIKE_THRESHOLD,
) -> list:
    """
    Connect a neural population to a muscle with activation callback.

    Parameters
    ----------
    source : str
        Name of motor neuron population.
    muscle : object
        Muscle object for force generation.
    activation_callback : Callable
        Function called when motor neurons fire.
        Expected signature: callback(neuron_id, muscle, delay_time)
    weight__uS : float, optional
        Synaptic weight in microsiemens, by default 1.0.
    threshold__mV : float, optional
        Spike threshold in millivolts, by default -10.0.

    Returns
    -------
    list[h.NetCon]
        List of motor neuron to muscle NetCon objects.
    """
    if source not in self.populations:
        raise ValueError(f"Source population '{source}' not found")

    # Rescale Quantities to NEURON's native units (or pass plain floats through)
    weight_value = _to_float(weight__uS, pq.uS)
    delay_value = _to_float(delay__ms, pq.ms)
    threshold_value = _to_float(threshold__mV, pq.mV)

    # Extract spike recording vectors for source population
    id_vector = None
    spike_vector = None
    if self.spike_recording:
        id_vector = self.spike_recording.get("idvec", {}).get(source)
        spike_vector = self.spike_recording.get("spkvec", {}).get(source)

    # Create muscle connections using existing infrastructure
    netcons = _connect_populations(
        populations=self.populations,
        source_pop=source,
        target_pop=None,  # External target
        connection_probability=1.0,  # All motor neurons connect
        muscle_callback=activation_callback,
        muscle=muscle,
        synaptic_delay=delay_value,
        synaptic_weight=weight_value,
        spike_threshold=threshold_value,
        id_vector=id_vector,
        spike_vector=spike_vector,
    )

    self.connections.append(
        {
            "type": "muscle",
            "source": source,
            "target": "muscle",
            "muscle": muscle,
            "callback": activation_callback,
            "weight__uS": weight_value,
            "delay__ms": delay_value,
            "threshold__mV": threshold_value,
        }
    )
    self._netcons_by_connection[(source, "muscle")] = netcons

    return netcons

connect_from_external

connect_from_external(source: str, target: str, weight__uS: Quantity__uS = DEFAULT_SYNAPTIC_WEIGHT, delay__ms: Quantity__ms = DEFAULT_SYNAPTIC_DELAY, threshold__mV: Quantity__mV = DEFAULT_SPIKE_THRESHOLD) -> list

Connect external input source to a neural population.

Parameters:

Name Type Description Default
source str

Name/label for external input source (e.g., "spindle", "cortical_drive").

required
target str

Name of target neural population.

required
weight__uS Quantity__uS

Synaptic weight in microsiemens, by default 0.8.

DEFAULT_SYNAPTIC_WEIGHT
delay__ms float

Synaptic delay in milliseconds, by default 1.0.

DEFAULT_SYNAPTIC_DELAY
threshold__mV float

Spike threshold in millivolts, by default -10.0.

DEFAULT_SPIKE_THRESHOLD

Returns:

Type Description
list[NetCon]

List of external to neural NetCon objects.

Source code in myogen/simulator/neuron/network.py
def connect_from_external(
    self,
    source: str,
    target: str,
    weight__uS: Quantity__uS = DEFAULT_SYNAPTIC_WEIGHT,
    delay__ms: Quantity__ms = DEFAULT_SYNAPTIC_DELAY,
    threshold__mV: Quantity__mV = DEFAULT_SPIKE_THRESHOLD,
) -> list:
    """
    Connect external input source to a neural population.

    Parameters
    ----------
    source : str
        Name/label for external input source (e.g., "spindle", "cortical_drive").
    target : str
        Name of target neural population.
    weight__uS : Quantity__uS, optional
        Synaptic weight in microsiemens, by default 0.8.
    delay__ms : float, optional
        Synaptic delay in milliseconds, by default 1.0.
    threshold__mV : float, optional
        Spike threshold in millivolts, by default -10.0.

    Returns
    -------
    list[h.NetCon]
        List of external to neural NetCon objects.
    """
    if target not in self.populations:
        raise ValueError(f"Target population '{target}' not found")

    # Rescale Quantities to NEURON's native units (or pass plain floats through)
    weight_value = _to_float(weight__uS, pq.uS)
    delay_value = _to_float(delay__ms, pq.ms)
    threshold_value = _to_float(threshold__mV, pq.mV)

    # Create external NetCons manually to maintain individual access
    from neuron import h

    netcons = []
    target_neurons = self.populations[target]

    for target_neuron in target_neurons:
        # Create NetCon from None (external source) to target neuron
        nc = h.NetCon(None, target_neuron.ns)
        nc.weight[0] = weight_value
        nc.delay = delay_value
        nc.threshold = threshold_value
        netcons.append(nc)

    self.connections.append(
        {
            "type": "external",
            "source": source,
            "target": target,
            "weight__uS": weight_value,
            "delay__ms": delay_value,
            "threshold__mV": threshold_value,
        }
    )
    self._netcons_by_connection[(source, target)] = netcons

    return netcons

connect_one_to_one

connect_one_to_one(source: str, target: str, probability: float = 1.0, weight__uS: Quantity__uS = DEFAULT_SYNAPTIC_WEIGHT, delay__ms: Quantity__ms = DEFAULT_SYNAPTIC_DELAY, threshold__mV: Quantity__mV = DEFAULT_SPIKE_THRESHOLD, inhibitory: bool = False) -> list

Connect two neural populations with one-to-one mapping.

Creates individual connections between source[i] and target[i] for each neuron pair at matching indices with specified probability. This is particularly useful for modeling independent noise sources (e.g., independent Poisson drives) where each target neuron should receive input from exactly one source neuron.

Parameters:

Name Type Description Default
source str

Name of source population (must exist in populations dict).

required
target str

Name of target population (must exist in populations dict).

required
probability float

Probability that each source[i] -> target[i] connection is made, by default 1.0. Must be between 0.0 and 1.0.

1.0
weight__uS Quantity__uS

Synaptic weight in microsiemens, by default 0.6.

DEFAULT_SYNAPTIC_WEIGHT
delay__ms Quantity__ms

Synaptic delay in milliseconds, by default 1.0.

DEFAULT_SYNAPTIC_DELAY
threshold__mV Quantity__mV

Spike threshold in millivolts, by default -10.0.

DEFAULT_SPIKE_THRESHOLD
inhibitory bool

If True, connect to inhibitory synapses on target neurons (reversal < -40 mV). If False, connect to excitatory synapses (reversal >= -40 mV). Default False.

False

Returns:

Type Description
list[NetCon]

List of created NEURON NetCon objects for connections that were made.

Raises:

Type Description
ValueError

If source or target populations don't exist, have different sizes, or probability is not in [0.0, 1.0].

Examples:

>>> # Create independent noise for each motor neuron
>>> noise_pool = DescendingDrive__Pool(n=10, timestep__ms=0.05)
>>> mn_pool = AlphaMN__Pool(n=10)
>>> network = Network({"noise": noise_pool, "mn": mn_pool})
>>> network.connect_one_to_one("noise", "mn", weight__uS=0.5)
Source code in myogen/simulator/neuron/network.py
def connect_one_to_one(
    self,
    source: str,
    target: str,
    probability: float = 1.0,
    weight__uS: Quantity__uS = DEFAULT_SYNAPTIC_WEIGHT,
    delay__ms: Quantity__ms = DEFAULT_SYNAPTIC_DELAY,
    threshold__mV: Quantity__mV = DEFAULT_SPIKE_THRESHOLD,
    inhibitory: bool = False,
) -> list:
    """
    Connect two neural populations with one-to-one mapping.

    Creates individual connections between source[i] and target[i] for each neuron
    pair at matching indices with specified probability. This is particularly useful
    for modeling independent noise sources (e.g., independent Poisson drives) where
    each target neuron should receive input from exactly one source neuron.

    Parameters
    ----------
    source : str
        Name of source population (must exist in populations dict).
    target : str
        Name of target population (must exist in populations dict).
    probability : float, optional
        Probability that each source[i] -> target[i] connection is made, by default 1.0.
        Must be between 0.0 and 1.0.
    weight__uS : Quantity__uS, optional
        Synaptic weight in microsiemens, by default 0.6.
    delay__ms : Quantity__ms, optional
        Synaptic delay in milliseconds, by default 1.0.
    threshold__mV : Quantity__mV, optional
        Spike threshold in millivolts, by default -10.0.
    inhibitory : bool, optional
        If True, connect to inhibitory synapses on target neurons (reversal < -40 mV).
        If False, connect to excitatory synapses (reversal >= -40 mV). Default False.

    Returns
    -------
    list[h.NetCon]
        List of created NEURON NetCon objects for connections that were made.

    Raises
    ------
    ValueError
        If source or target populations don't exist, have different sizes,
        or probability is not in [0.0, 1.0].

    Examples
    --------
    >>> # Create independent noise for each motor neuron
    >>> noise_pool = DescendingDrive__Pool(n=10, timestep__ms=0.05)
    >>> mn_pool = AlphaMN__Pool(n=10)
    >>> network = Network({"noise": noise_pool, "mn": mn_pool})
    >>> network.connect_one_to_one("noise", "mn", weight__uS=0.5)
    """
    # Validation
    if source not in self.populations:
        raise ValueError(f"Source population '{source}' not found")
    if target not in self.populations:
        raise ValueError(f"Target population '{target}' not found")
    if not 0.0 <= probability <= 1.0:
        raise ValueError(f"Probability must be 0.0-1.0, got {probability}")

    # Rescale Quantities to NEURON's native units (or pass plain floats through)
    weight_value = _to_float(weight__uS, pq.uS)
    delay_value = _to_float(delay__ms, pq.ms)
    threshold_value = _to_float(threshold__mV, pq.mV)

    # Extract spike recording vectors for source population
    id_vector = None
    spike_vector = None
    if self.spike_recording:
        id_vector = self.spike_recording.get("idvec", {}).get(source)
        spike_vector = self.spike_recording.get("spkvec", {}).get(source)

    # Create one-to-one connections using new helper function
    netcons = _connect_one_to_one(
        source_pop=source,
        target_pop=target,
        populations=self.populations,
        connection_probability=probability,
        synaptic_delay=delay_value,
        synaptic_weight=weight_value,
        spike_threshold=threshold_value,
        id_vector=id_vector,
        spike_vector=spike_vector,
        inhibitory=inhibitory,
    )

    self.connections.append(
        {
            "type": "one_to_one",
            "source": source,
            "target": target,
            "probability": probability,
            "weight__uS": weight_value,
            "delay__ms": delay_value,
            "threshold__mV": threshold_value,
            "inhibitory": inhibitory,
        }
    )
    self._netcons_by_connection[(source, target)] = netcons

    return netcons

get_connections

get_connections() -> list[dict]

Get list of all connection specifications.

Source code in myogen/simulator/neuron/network.py
def get_connections(self) -> list[dict]:
    """Get list of all connection specifications."""
    return self.connections.copy()

get_netcons

get_netcons(source: Optional[str] = None, target: Optional[str] = None) -> list

Get NEURON NetCon objects with optional filtering by source and target.

Parameters:

Name Type Description Default
source str

Filter by source population/input name. If None, returns NetCons from all sources.

None
target str

Filter by target population name. If None, returns NetCons to all targets.

None

Returns:

Type Description
list[NetCon]

List of matching NetCon objects.

Source code in myogen/simulator/neuron/network.py
def get_netcons(self, source: Optional[str] = None, target: Optional[str] = None) -> list:
    """
    Get NEURON NetCon objects with optional filtering by source and target.

    Parameters
    ----------
    source : str, optional
        Filter by source population/input name. If None, returns NetCons from all sources.
    target : str, optional
        Filter by target population name. If None, returns NetCons to all targets.

    Returns
    -------
    list[h.NetCon]
        List of matching NetCon objects.
    """
    if source is None and target is None:
        all_netcons = []
        for netcon_list in self._netcons_by_connection.values():
            all_netcons.extend(netcon_list)
        return all_netcons

    matching_netcons = []

    for (
        conn_source,
        conn_target,
    ), netcon_list in self._netcons_by_connection.items():
        source_matches = source is None or conn_source == source
        target_matches = target is None or conn_target == target

        if source_matches and target_matches:
            matching_netcons.extend(netcon_list)

    return matching_netcons

print_network

print_network()

Print a summary of network structure.

Source code in myogen/simulator/neuron/network.py
def print_network(self):
    """Print a summary of network structure."""
    print(f"Network with {len(self.populations)} populations:")
    for name, neurons in self.populations.items():
        print(f"  {name}: {len(neurons)} neurons")

    print(f"\nConnections ({len(self.connections)}):")
    for i, conn in enumerate(self.connections):
        if conn["type"] == "neural":
            print(
                f"  {i + 1}. {conn['source']}{conn['target']} "
                f"(p={conn['probability']}, w={conn['weight__uS']}uS)"
            )
        elif conn["type"] == "muscle":
            print(f"  {i + 1}. {conn['source']} → muscle (w={conn['weight__uS']}uS)")
        elif conn["type"] == "external":
            print(f"  {i + 1}. {conn['source']}{conn['target']} (w={conn['weight__uS']}uS)")
        elif conn["type"] == "one_to_one":
            print(
                f"  {i + 1}. {conn['source']}{conn['target']} [1-to-1] "
                f"(p={conn['probability']}, w={conn['weight__uS']}uS)"
            )

SimulationRunner

SimulationRunner(network: Network, models: dict[str, Any], step_callback: Callable[[Any], Any], model_outputs: Optional[dict[str, Union[list[str], None]]] = None, temperature__celsius: float = 36.0)

Manages NEURON simulation execution with automated setup, initialization, and result collection for neuromuscular simulations.

Provides a clean interface for running complex neuromuscular simulations while maintaining full user control over populations, connections, and step-by-step simulation logic. Automatically handles NEURON environment setup, voltage initialization, and structured result collection.

Separates simulation control from plotting and analysis concerns.

Initialize SimulationRunner with network, models, and step callback.

Parameters:

Name Type Description Default
network Network

Configured Network instance with populations and connections.

required
models Dict[str, Any]

Physiological models (e.g., {"hill": hill_model, "spin": spindle_model}).

required
step_callback Callable

User-defined function called at each simulation timestep.

required
model_outputs Optional[Dict[str, Union[List[str], None]]]

Explicit model output attributes to collect. None uses smart defaults. Format: {"model_name": ["attr1", "attr2"]} or {"model_name": None} for defaults, by default None.

None
temperature__celsius float

NEURON simulation temperature, by default 36.0.

36.0
Source code in myogen/simulator/neuron/simulation_runner.py
def __init__(
    self,
    network: Network,
    models: dict[str, Any],
    step_callback: Callable[[Any], Any],
    model_outputs: Optional[dict[str, Union[list[str], None]]] = None,
    temperature__celsius: float = 36.0,
):
    """
    Initialize SimulationRunner with network, models, and step callback.

    Parameters
    ----------
    network : Network
        Configured Network instance with populations and connections.
    models : Dict[str, Any]
        Physiological models (e.g., {"hill": hill_model, "spin": spindle_model}).
    step_callback : Callable
        User-defined function called at each simulation timestep.
    model_outputs : Optional[Dict[str, Union[List[str], None]]], optional
        Explicit model output attributes to collect. None uses smart defaults.
        Format: {"model_name": ["attr1", "attr2"]} or {"model_name": None}
        for defaults, by default None.
    temperature__celsius : float, optional
        NEURON simulation temperature, by default 36.0.
    """
    # Store immutable parameters following project pattern
    self.network = network
    self.populations = network.populations  # Expose populations from network
    self.models = models
    self.step_callback = step_callback
    self.model_outputs = model_outputs
    self.temperature__celsius = temperature__celsius

    # Private working copies
    self._network = network
    self._populations = network.populations  # Get populations from network
    self._models = models
    self._step_callback = step_callback
    self._model_outputs = self._resolve_model_outputs()
    self._temperature__celsius = temperature__celsius

    # Runtime state
    self._trace_vectors: dict[str, dict[int, Any]] = {}
    self._step_counter = None
    self._progress_bar = None
    self._total_steps = None

    # Setup internal spike recording vectors
    self._spike_recording = self._setup_spike_recording()

run

run(duration__ms: Quantity__ms, timestep__ms: Quantity__ms, membrane_recording: Optional[dict[str, list[int]]] = None, verbose: bool = True) -> Block

Execute NEURON simulation with automated setup and result collection.

Parameters:

Name Type Description Default
duration__ms Quantity__ms

Total simulation duration in milliseconds.

required
timestep__ms Quantity__ms

Integration timestep in milliseconds.

required
membrane_recording Optional[Dict[str, List[int]]]

Populations and cell indices for membrane potential recording. Format: {"population_name": [cell_id1, cell_id2, ...]}, by default None.

None
verbose bool

If True, display progress bar and status messages. Set to False to disable.

True

Returns:

Type Description
Block

Structured simulation results containing: - spikes: Spike timing and ID data for all populations - membrane: Membrane potential traces (if requested) - models: Output data from all physiological models - simulation: Time vector and simulation metadata

Raises:

Type Description
ValueError

If model output attributes don't exist on model instances.

RuntimeError

If NEURON simulation fails to complete.

Source code in myogen/simulator/neuron/simulation_runner.py
def run(
    self,
    duration__ms: Quantity__ms,
    timestep__ms: Quantity__ms,
    membrane_recording: Optional[dict[str, list[int]]] = None,
    verbose: bool = True,
) -> Block:
    """
    Execute NEURON simulation with automated setup and result collection.

    Parameters
    ----------
    duration__ms : Quantity__ms
        Total simulation duration in milliseconds.
    timestep__ms : Quantity__ms
        Integration timestep in milliseconds.
    membrane_recording : Optional[Dict[str, List[int]]], optional
        Populations and cell indices for membrane potential recording.
        Format: {"population_name": [cell_id1, cell_id2, ...]}, by default None.
    verbose : bool, default=True
        If True, display progress bar and status messages. Set to False to disable.

    Returns
    -------
    Block
        Structured simulation results containing:
        - spikes: Spike timing and ID data for all populations
        - membrane: Membrane potential traces (if requested)
        - models: Output data from all physiological models
        - simulation: Time vector and simulation metadata

    Raises
    ------
    ValueError
        If model output attributes don't exist on model instances.
    RuntimeError
        If NEURON simulation fails to complete.
    """
    try:
        # Setup NEURON environment
        self._setup_neuron_environment(duration__ms, timestep__ms, verbose=verbose)

        # Setup optional membrane recording
        if membrane_recording:
            self._setup_membrane_recording(membrane_recording)

        # Initialize population voltages
        self._initialize_voltages()

        # Register step callback for closed-loop dynamics
        self._register_step_callback()

        # Validate model outputs before simulation
        self._validate_model_outputs()

        # Setup spike recording on network
        self._setup_network_spike_recording()

        h.run()

        # Close progress bar (with error handling)
        if self._progress_bar is not None:
            try:
                self._progress_bar.close()
            except (TypeError, AttributeError):
                # Ignore progress bar closing errors
                pass

        if verbose:
            print("Simulation completed")

        # Collect and structure results
        results = self._collect_results(duration__ms, timestep__ms)

        return results

    except Exception as e:
        # Close progress bar in case of error (with error handling)
        if self._progress_bar is not None:
            try:
                self._progress_bar.close()
            except (TypeError, AttributeError):
                # Ignore progress bar closing errors
                pass
        raise RuntimeError(f"Simulation failed: {str(e)}") from e

get_model_outputs

get_model_outputs(model_name: str) -> list[str]

Get the list of output attributes that will be collected for a model.

Parameters:

Name Type Description Default
model_name str

Name of the model as specified in the models dictionary.

required

Returns:

Type Description
List[str]

List of attribute names that will be collected from this model.

Source code in myogen/simulator/neuron/simulation_runner.py
def get_model_outputs(self, model_name: str) -> list[str]:
    """
    Get the list of output attributes that will be collected for a model.

    Parameters
    ----------
    model_name : str
        Name of the model as specified in the models dictionary.

    Returns
    -------
    List[str]
        List of attribute names that will be collected from this model.
    """
    return self._model_outputs.get(model_name, [])

set_model_outputs

set_model_outputs(model_name: str, output_attrs: list[str]) -> None

Override the output attributes for a specific model.

Parameters:

Name Type Description Default
model_name str

Name of the model as specified in the models dictionary.

required
output_attrs List[str]

List of attribute names to collect from this model.

required

Raises:

Type Description
ValueError

If model_name is not found in the models dictionary.

Source code in myogen/simulator/neuron/simulation_runner.py
def set_model_outputs(self, model_name: str, output_attrs: list[str]) -> None:
    """
    Override the output attributes for a specific model.

    Parameters
    ----------
    model_name : str
        Name of the model as specified in the models dictionary.
    output_attrs : List[str]
        List of attribute names to collect from this model.

    Raises
    ------
    ValueError
        If model_name is not found in the models dictionary.
    """
    if model_name not in self._models:
        raise ValueError(f"Model '{model_name}' not found in models")

    self._model_outputs[model_name] = output_attrs

Muscle & force

Muscle

Muscle(recruitment_thresholds: RECRUITMENT_THRESHOLDS__ARRAY, radius__mm: Quantity__mm = 6.91 * mm, length__mm: Quantity__mm = 30.0 * mm, fiber_density__fibers_per_mm2: Quantity__per_mm2 = 350 * mm ** -2, max_innervation_area_to_total_muscle_area__ratio: float = 1 / 4, mean_conduction_velocity__m_per_s: Quantity__m_per_s = 4.2 * m / s, mean_fiber_length__mm: Quantity__mm = 31.7 * mm, var_fiber_length__mm: Quantity__mm = 2.8 * mm, radius_bone__mm: Quantity__mm = 0 * mm, fat_thickness__mm: Quantity__mm = 0.3 * mm, skin_thickness__mm: Quantity__mm = 1.29 * mm, muscle_conductivity_radial__S_per_m: Quantity__S_per_m = 0.09 * S / m, muscle_conductivity_longitudinal__S_per_m: Quantity__S_per_m = 0.4 * S / m, fat_conductivity__S_per_m: Quantity__S_per_m = 0.0407 * S / m, skin_conductivity__S_per_m: Quantity__S_per_m = 0.000488 * S / m, grid_resolution: int = 256, autorun: bool = False)

A muscle model based on the cylindrical description of the volume conductor by Farina et al. 2004 [1] and the motor unit distribution by Konstantin et al. 2020 [2].

All default values are set to simulate the First Dorsal Interosseous (FDI) muscle. Values are pulled from the literature.

Parameters:

Name Type Description Default
recruitment_thresholds RECRUITMENT_THRESHOLDS__ARRAY

Array of recruitment thresholds for each motor unit (see myogen.simulator.RecruitmentThresholds). Values range from 0 to 1 with the largest motor units having thresholds near 1.

required
radius__mm float

Radius of the muscle cross-section in millimeters. Default is set to 6.91 mm as determined by Jacobson et al. 1992 [3].

6.91
length__mm float

Length of the muscle in millimeters. Default is 30.0 mm, a nominal value chosen to match the order of magnitude of the FDI muscle; adjust to match the muscle under study.

30.0
fiber_density__fibers_per_mm2 float

Density of muscle fibers per square millimeter. Default is set to 350 fibers/mm² as determined by Bettelho et al. 2019 [7].

350
max_innervation_area_to_total_muscle_area__ratio float

Ratio defining the maximum territory size relative to total muscle area. Default is 0.25 as a pragmatic upper bound for the FDI, with no single published source; revisit for larger muscles. A value of 0.25 means the largest motor unit can innervate up to 25% of the total muscle cross-sectional area. Must be in range (0, 1].

0.25
mean_conduction_velocity__m_per_s float

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.

4.2
mean_fiber_length__mm float

Mean fiber length in mm. Default is set to 31.7 mm as determined by Jacobson et al. 1992 [3] (Table 1).

31.7
var_fiber_length__mm float

Fiber length variance in mm. Default is set to 2.8 mm as determined by Jacobson et al. 1992 [3] (Table 1).

2.8
radius_bone__mm float

Bone radius in mm. Default is set to 1 mm.

1
fat_thickness__mm float

Fat thickness in mm. Default is set to 0.3 mm as determined by Störchle et al. 2018 [5].

0.3
skin_thickness__mm float

Skin thickness in mm. Default is set to the male skin thickness average of 1.29 mm as determined by Brodar 1960 [6].

1.29
muscle_conductivity_radial__S_per_m float

Muscle conductivity in radial direction. Default is set to 0.09 S/m as determined by Botelho et al. 2019 [7] (Table 1).

0.09
muscle_conductivity_longitudinal__S_per_m float

Muscle conductivity in longitudinal direction. Default is set to 0.4 S/m as determined by Botelho et al. 2019 [7] (Table 1).

0.4
fat_conductivity__S_per_m float

Fat conductivity. Default is set to 4.07E-2 S/m as determined by Botelho et al. 2019 [7] (Table 1).

4.07E-2
skin_conductivity__S_per_m float

Skin conductivity. Default is set to 4.88E-4 S/m as determined by Botelho et al. 2019 [7] (Table 1).

4.88E-4
grid_resolution int

Resolution of the computational grid used for innervation the muscle. Higher values provide more accurate spatial distribution but increase computational cost. Default is set to 256.

256
autorun bool

If True, automatically executes the complete muscle simulation pipeline: innervation distribution, muscle fiber generation, and fiber-to-motor unit assignment. If False, these steps must be called manually.

False

Attributes:

Name Type Description
innervation_center_positions__mm ndarray

Motor unit innervation center positions [x, y] in mm. Available after distribute_innervation_centers().

muscle_fiber_centers__mm ndarray

Muscle fiber center positions [x, y] in mm. Available after generate_muscle_fiber_centers().

muscle_fiber_diameters__mm ndarray

Muscle fiber diameters in mm. Available after _generate_fiber_properties().

muscle_fiber_conduction_velocities__mm_per_s ndarray

Muscle fiber conduction velocities in mm/s. Available after _generate_fiber_properties().

assignment ndarray

Motor unit assignment for each muscle fiber. Available after assign_mfs2mns().

number_of_muscle_fibers int

Total number of muscle fibers. Available after generate_muscle_fiber_centers().

muscle_border__mm ndarray

Muscle boundary points for visualization. Available after generate_muscle_fiber_centers().

resulting_number_of_innervated_fibers ndarray

Actual number of fibers per motor unit. Available after assign_mfs2mns().

resulting_innervation_areas__mm2 ndarray

Actual innervation areas per motor unit in mm². Available after assign_mfs2mns().

Raises:

Type Description
ValueError

If max_innervation_area_to_total_muscle_area__ratio is not in (0, 1].

References

[1] Farina, D., Mesin, L., Martina, S., Merletti, R., 2004. A surface EMG generation model with multilayer cylindrical description of the volume conductor. IEEE Transactions on Biomedical Engineering 51, 415–426. https://doi.org/10.1109/TBME.2003.820998
[2] Konstantin, A., Yu, T., Le Carpentier, E., Aoustin, Y., Farina, D., 2020. Simulation of Motor Unit Action Potential Recordings From Intramuscular Multichannel Scanning Electrodes. IEEE Transactions on Biomedical Engineering 67, 2005–2014. https://doi.org/10.1109/TBME.2019.2953680
[3] Jacobson, M.D., Raab, R., Fazeli, B.M., Abrams, R.A., Botte, M.J., Lieber, R.L., 1992. Architectural design of the human intrinsic hand muscles. The Journal of Hand Surgery 17, 804–809. https://doi.org/10.1016/0363-5023(92)90446-V
[4] Nishizono, H., Fujimoto, T., Ohtake, H., Miyashita, M., 1990. Muscle fiber conduction velocity and contractile properties estimated from surface electrode arrays. Electroencephalography and Clinical Neurophysiology 75, 75–81. https://doi.org/10.1016/0013-4694(90)90154-C
[5] Störchle, P., Müller, W., Sengeis, M., Lackner, S., Holasek, S., Fürhapter-Rieger, A., 2018. Measurement of mean subcutaneous fat thickness: eight standardised ultrasound sites compared to 216 randomly selected sites. Sci Rep 8, 16268. https://doi.org/10.1038/s41598-018-34213-0
[6] Brodar, V., 1960. Observations on skin thickness and subcutaneous tissue in man. Zeitschrift für Morphologie und Anthropologie 50, 386–395.
[7] Botelho, D.P., Curran, K., Lowery, M.M., 2019. Anatomically accurate model of EMG during index finger flexion and abduction derived from diffusion tensor imaging. PLOS Computational Biology 15, e1007267. https://doi.org/10.1371/journal.pcbi.1007267

Source code in myogen/simulator/core/muscle/muscle.py
def __init__(
    self,
    recruitment_thresholds: RECRUITMENT_THRESHOLDS__ARRAY,
    radius__mm: Quantity__mm = 6.91 * pq.mm,
    length__mm: Quantity__mm = 30.0 * pq.mm,
    fiber_density__fibers_per_mm2: Quantity__per_mm2 = 350 * pq.mm**-2,
    max_innervation_area_to_total_muscle_area__ratio: float = 1 / 4,
    mean_conduction_velocity__m_per_s: Quantity__m_per_s = 4.2 * pq.m / pq.s,
    mean_fiber_length__mm: Quantity__mm = 31.7 * pq.mm,
    var_fiber_length__mm: Quantity__mm = 2.8 * pq.mm,
    radius_bone__mm: Quantity__mm = 0 * pq.mm,
    fat_thickness__mm: Quantity__mm = 0.3 * pq.mm,
    skin_thickness__mm: Quantity__mm = 1.29 * pq.mm,
    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,
    grid_resolution: int = 256,
    autorun: bool = False,
) -> None:
    # Muscle properties - immutable public access
    self.radius__mm = radius__mm
    self.length__mm = length__mm
    self.fiber_density__fibers_per_mm2 = fiber_density__fibers_per_mm2
    self.max_innervation_area_to_total_muscle_area__ratio = (
        max_innervation_area_to_total_muscle_area__ratio
    )
    self.mean_conduction_velocity__m_s = mean_conduction_velocity__m_per_s
    self.mean_fiber_length__mm = mean_fiber_length__mm
    self.var_fiber_length__mm = var_fiber_length__mm
    self.radius_bone__mm = radius_bone__mm
    self.fat_thickness__mm = fat_thickness__mm
    self.skin_thickness__mm = skin_thickness__mm
    self.muscle_conductivity_radial__S_m = muscle_conductivity_radial__S_per_m
    self.muscle_conductivity_longitudinal__S_m = muscle_conductivity_longitudinal__S_per_m
    self.fat_conductivity__S_m = fat_conductivity__S_per_m
    self.skin_conductivity__S_m = skin_conductivity__S_per_m
    self.grid_resolution = grid_resolution
    self.autorun = autorun
    # Private copies for internal modifications
    self._radius__mm = radius__mm
    self._length__mm = length__mm
    self._fiber_density__fibers_per_mm2 = fiber_density__fibers_per_mm2
    self._max_innervation_area_to_total_muscle_area__ratio = (
        max_innervation_area_to_total_muscle_area__ratio
    )
    self._mean_conduction_velocity__m_s = mean_conduction_velocity__m_per_s
    self._mean_fiber_length__mm = mean_fiber_length__mm
    self._var_fiber_length__mm = var_fiber_length__mm
    self._radius_bone__mm = radius_bone__mm
    self._fat_thickness__mm = fat_thickness__mm
    self._skin_thickness__mm = skin_thickness__mm
    self._muscle_conductivity_radial__S_m = muscle_conductivity_radial__S_per_m
    self._muscle_conductivity_longitudinal__S_m = muscle_conductivity_longitudinal__S_per_m
    self._fat_conductivity__S_m = fat_conductivity__S_per_m
    self._skin_conductivity__S_m = skin_conductivity__S_per_m
    self._grid_resolution = grid_resolution
    self._autorun = autorun
    self._recruitment_thresholds = recruitment_thresholds.copy()

    # Derived properties
    self.muscle_area__mm2 = np.pi * (self._radius__mm**2)
    self.max_innervation_area_scaling_factor = (
        1 / self._max_innervation_area_to_total_muscle_area__ratio
    )
    self._number_of_neurons = len(self._recruitment_thresholds)

    # Simulation results - stored privately, accessed via properties
    self._innervation_center_positions__mm: Optional[Quantity__mm] = None
    self._muscle_fiber_centers__mm: Optional[Quantity__mm] = None
    self._assignment: Optional[np.ndarray] = None
    self._muscle_fiber_diameters__mm: Optional[Quantity__mm] = None
    self._muscle_fiber_conduction_velocities__mm_per_s: Optional[Quantity__mm_per_s] = None
    self._number_of_muscle_fibers: Optional[int] = None
    self._muscle_border__mm: Optional[Quantity__mm] = None

    # Validate the ratio
    if not (0 < max_innervation_area_to_total_muscle_area__ratio <= 1):
        raise ValueError(
            '"max_innervation_area_to_total_muscle_area__ratio" must be in (0, 1]. '
            "This ratio defines how much of the muscle area the largest motor unit can occupy. "
            "For realistic simulations, try values between 0.1 and 0.5."
        )

    self.desired_innervation_areas__mm2 = (
        self._recruitment_thresholds
        / np.max(self._recruitment_thresholds)
        * self.muscle_area__mm2
        / self.max_innervation_area_scaling_factor
    )

    self.desired_number_of_innervated_fibers = np.round(
        (
            self.desired_innervation_areas__mm2
            / np.sum(self.desired_innervation_areas__mm2)
            * self.muscle_area__mm2
            * self._fiber_density__fibers_per_mm2
        ).magnitude
    ).astype(int)

    if autorun:
        self.distribute_innervation_centers()
        self.generate_muscle_fiber_centers()
        self.assign_mfs2mns()
        self._generate_fiber_properties()

resulting_number_of_innervated_fibers property

resulting_number_of_innervated_fibers: ndarray

Calculate the actual number of muscle fibers assigned to each motor unit.

This property returns the final fiber counts after the assignment process, which may differ slightly from the desired counts due to the stochastic assignment algorithm and discrete fiber distribution.

Returns:

Type Description
ndarray

Array of length n_motor_units where each element represents the actual number of muscle fibers assigned to the corresponding motor unit. The sum of all elements equals the total number of muscle fibers.

Raises:

Type Description
ValueError

If muscle fiber assignment has not been completed yet.

Examples:

>>> actual_counts = muscle.resulting_number_of_innervated_fibers
>>> desired_counts = muscle.desired_number_of_innervated_fibers
>>> print(f"Motor unit 0: desired {desired_counts[0]}, actual {actual_counts[0]}")
Notes

This property can be used to assess how well the assignment algorithm achieved the target fiber distribution. Large deviations may indicate the need to adjust assignment parameters or increase grid resolution.

resulting_innervation_areas__mm2 property

resulting_innervation_areas__mm2: Quantity__mm2

Calculate the actual innervation areas for each motor unit based on assigned fibers.

The innervation area is computed as the area of a circle that encompasses all muscle fibers assigned to a motor unit, centered on the motor unit's innervation center. This provides a measure of the spatial extent of each motor unit territory.

Returns:

Type Description
Quantity__mm2

Array of length n_motor_units containing the innervation area (in mm²) for each motor unit. Areas are calculated as π × r², where r is the maximum distance from the innervation center to any assigned fiber.

Raises:

Type Description
ValueError

If innervation_center_positions is None or assignment has not been completed.

Examples:

>>> actual_areas = muscle.resulting_innervation_areas__mm2
>>> desired_areas = muscle.desired_innervation_areas__mm2
>>> for i, (actual, desired) in enumerate(zip(actual_areas, desired_areas)):
...     print(f"MU {i}: desired {desired:.2f} mm², actual {actual:.2f} mm²")
Notes

The resulting areas may differ from desired areas due to the discrete nature of fiber assignment and the constraint of the circular muscle boundary. Motor units near the muscle periphery may have smaller actual areas than desired due to boundary effects.

innervation_center_positions__mm property

innervation_center_positions__mm: Quantity__mm

Motor unit innervation center positions [x, y] in mm.

Returns:

Type Description
Quantity__mm

Array of shape (n_motor_units, 2) containing [x, y] coordinates in mm.

Raises:

Type Description
ValueError

If innervation centers have not been computed yet.

muscle_fiber_centers__mm property

muscle_fiber_centers__mm: Quantity__mm

Muscle fiber center positions [x, y] in mm.

Returns:

Type Description
Quantity__mm

Array of shape (n_fibers, 2) containing [x, y] coordinates in mm.

Raises:

Type Description
ValueError

If muscle fiber centers have not been computed yet.

muscle_fiber_diameters__mm property

muscle_fiber_diameters__mm: Quantity__mm

Muscle fiber diameters in mm.

Returns:

Type Description
Quantity__mm

Array of muscle fiber diameters in mm.

Raises:

Type Description
ValueError

If fiber properties have not been computed yet.

muscle_fiber_conduction_velocities__mm_per_s property

muscle_fiber_conduction_velocities__mm_per_s: Quantity__mm_per_s

Muscle fiber conduction velocities in mm/s.

Returns:

Type Description
Quantity__mm_per_s

Array of muscle fiber conduction velocities in mm/s.

Raises:

Type Description
ValueError

If fiber properties have not been computed yet.

assignment property

assignment: ndarray

Motor unit assignment for each muscle fiber.

Returns:

Type Description
ndarray

Array where each element indicates the motor unit index (0 to n_motor_units-1) assigned to that fiber.

Raises:

Type Description
ValueError

If muscle fiber assignment has not been completed yet.

number_of_muscle_fibers property

number_of_muscle_fibers: int

Total number of muscle fibers.

Returns:

Type Description
int

Total number of muscle fibers.

Raises:

Type Description
ValueError

If muscle fiber centers have not been computed yet.

muscle_border__mm property

muscle_border__mm: Quantity__mm

Muscle boundary points for visualization.

Returns:

Type Description
Quantity__mm

Array of boundary points for the circular muscle cross-section in mm.

Raises:

Type Description
ValueError

If muscle fiber centers have not been computed yet.

recruitment_thresholds property

recruitment_thresholds: RECRUITMENT_THRESHOLDS__ARRAY

Motor unit recruitment thresholds.

Returns:

Type Description
ndarray

Array of recruitment thresholds for each motor unit.

distribute_innervation_centers

distribute_innervation_centers() -> None

Distribute innervation center positions using the fast marching method.

This method implements an optimal packing algorithm to distribute motor unit innervation centers within the circular muscle cross-section. The algorithm uses the Fast Marching Method to ensure that each new innervation center is placed at the location that maximizes the minimum distance to all previously placed centers.

Results are stored in the innervation_center_positions property after execution.

Notes

This method must be called before generate_muscle_fiber_centers() and assign_mfs2mns(). The resulting distribution approximates the optimal packing problem for circles, leading to realistic motor unit territory arrangements.

Source code in myogen/simulator/core/muscle/muscle.py
def distribute_innervation_centers(self) -> None:
    """
    Distribute innervation center positions using the fast marching method.

    This method implements an optimal packing algorithm to distribute motor unit
    innervation centers within the circular muscle cross-section. The algorithm
    uses the Fast Marching Method to ensure that each new innervation center is
    placed at the location that maximizes the minimum distance to all previously
    placed centers.

    Results are stored in the `innervation_center_positions` property after execution.

    Notes
    -----
    This method must be called before generate_muscle_fiber_centers() and
    assign_mfs2mns(). The resulting distribution approximates the optimal
    packing problem for circles, leading to realistic motor unit territory
    arrangements.
    """
    density_map = np.ones((self._grid_resolution, self._grid_resolution))
    X, Y = np.meshgrid(
        np.arange(self._grid_resolution),
        np.arange(self._grid_resolution),
    )
    density_map[
        np.sqrt((X - self._grid_resolution / 2) ** 2 + (Y - self._grid_resolution / 2) ** 2)
        > self._grid_resolution / 2 - 1
    ] = 1e-10

    vertices = np.zeros((2, self._number_of_neurons + 1))
    vertices[:, 0] = [1, 1]

    # MATLAB: for i = 2:(obj.N+1)
    for i in range(1, self._number_of_neurons + 1):
        # Use scikit-fmm for fast marching
        # Create speed map, avoiding division by zero
        ind = np.argmax(_perform_fast_marching(density_map.copy(), vertices[:, :i]))
        x, y = np.unravel_index(ind, (self._grid_resolution, self._grid_resolution))
        vertices[:, i] = [x, y]

    vertices = vertices * pq.mm

    # MATLAB: obj.innervation_center_positions = vertices(:,end:-1:2)';
    # This takes columns from end down to 2 (1-indexed), then transposes
    # In Python: vertices[:, -1:0:-1] gives us columns from end down to 1 (0-indexed)
    self._innervation_center_positions__mm = vertices[:, -1:0:-1].T

    # Only proceed if we have valid innervation_center_positions
    if (
        self._innervation_center_positions__mm.shape[0] > 0
        and self._innervation_center_positions__mm.shape[1] == 2
    ):
        center_offset = (
            self._innervation_center_positions__mm - (self._grid_resolution / 2) * pq.mm
        )
        max_dist = np.max(np.sqrt(center_offset[:, 0] ** 2 + center_offset[:, 1] ** 2))
        if max_dist > 0:  # Avoid division by zero
            self._innervation_center_positions__mm = (
                center_offset / max_dist
            ) * self._radius__mm
        else:
            self._innervation_center_positions__mm = (
                center_offset  # Keep original if max_dist is 0
            )

generate_muscle_fiber_centers

generate_muscle_fiber_centers(verbose: bool = True) -> None

Generate muscle fiber center positions using a pre-computed Voronoi distribution.

This method creates the spatial distribution of muscle fiber centers within the circular muscle cross-section. The distribution is based on a Voronoi tessellation pattern that mimics the natural packing of muscle fibers observed in histological studies.

Parameters:

Name Type Description Default
verbose bool

If True, display status messages. Set to False to disable.

True
Notes

Results are stored in the following properties after execution:

  • mf_centers: Array of shape (n_fibers, 2) with fiber positions [x, y] in mm
  • number_of_muscle_fibers: Total number of muscle fibers
  • muscle_border: Array of border points for visualization

This method should be called after distribute_innervation_centers() and before assign_mfs2mns(). The Voronoi-based distribution provides more realistic fiber spacing compared to regular grids or purely random distributions.

The reference dataset ('voronoi_pi1e5.csv') contains 100,000 pre-computed Voronoi cell centers optimized for circular domains, ensuring efficient and consistent fiber distributions across simulations.

Source code in myogen/simulator/core/muscle/muscle.py
def generate_muscle_fiber_centers(self, verbose: bool = True) -> None:
    """
    Generate muscle fiber center positions using a pre-computed Voronoi distribution.

    This method creates the spatial distribution of muscle fiber centers
    within the circular muscle cross-section. The distribution is based on a
    Voronoi tessellation pattern that mimics the natural packing of muscle fibers
    observed in histological studies.

    Parameters
    ----------
    verbose : bool, default=True
        If True, display status messages. Set to False to disable.

    Notes
    -----
    Results are stored in the following properties after execution:

    - `mf_centers`: Array of shape (n_fibers, 2) with fiber positions [x, y] in mm
    - `number_of_muscle_fibers`: Total number of muscle fibers
    - `muscle_border`: Array of border points for visualization

    This method should be called after distribute_innervation_centers() and
    before assign_mfs2mns(). The Voronoi-based distribution provides more
    realistic fiber spacing compared to regular grids or purely random distributions.

    The reference dataset ('voronoi_pi1e5.csv') contains 100,000 pre-computed
    Voronoi cell centers optimized for circular domains, ensuring efficient
    and consistent fiber distributions across simulations.
    """

    # Expected number of muscle fibers in the muscle
    self._number_of_muscle_fibers = int(
        np.rint(((self._radius__mm**2) * np.pi * self._fiber_density__fibers_per_mm2).magnitude)
    )

    self._muscle_fiber_centers__mm = (
        pd.read_csv(
            Path(inspect.getfile(self.__class__)).parent / "voronoi_pi1e5.csv",
            header=None,
        ).values
        * pq.mm
    )

    # Adjust the loaded innervation_center_positions to the expected number of fibers and muscle radius
    self._muscle_fiber_centers__mm = (self._muscle_fiber_centers__mm - (5 * pq.mm)) / 4
    dists = np.sqrt(
        self._muscle_fiber_centers__mm[:, 0] ** 2 + self._muscle_fiber_centers__mm[:, 1] ** 2
    )
    sorted_indices = np.argsort(dists)

    if len(sorted_indices) >= self._number_of_muscle_fibers + 1:
        self._muscle_fiber_centers__mm = (
            self._muscle_fiber_centers__mm[sorted_indices[: self._number_of_muscle_fibers], :]
            / dists[sorted_indices[self._number_of_muscle_fibers]]
            * self._radius__mm
        )
    else:
        self._muscle_fiber_centers__mm = (
            self._muscle_fiber_centers__mm[sorted_indices, :]
            / dists[sorted_indices[-1]]
            * self._radius__mm
        )
        self._number_of_muscle_fibers = len(self._muscle_fiber_centers__mm)

    # Remove fibers inside the bone boundary
    # Fibers should only exist in the muscle tissue, not in the bone core
    if self._radius_bone__mm.magnitude > 0:
        fiber_radial_dists = np.sqrt(
            self._muscle_fiber_centers__mm[:, 0] ** 2
            + self._muscle_fiber_centers__mm[:, 1] ** 2
        )
        valid_fiber_mask = fiber_radial_dists > self._radius_bone__mm
        n_fibers_removed = np.sum(~valid_fiber_mask)

        if n_fibers_removed > 0:
            self._muscle_fiber_centers__mm = self._muscle_fiber_centers__mm[valid_fiber_mask]
            self._number_of_muscle_fibers = len(self._muscle_fiber_centers__mm)
            if verbose:
                print(
                    f"Removed {n_fibers_removed} fibers inside bone radius "
                    f"(r < {self._radius_bone__mm.magnitude:.3f} mm)"
                )

    # Create muscle border for plotting
    phi_circle = np.linspace(0, 2 * np.pi, 1000)
    phi_circle = phi_circle[:-1]
    self._muscle_border__mm = (
        np.column_stack(
            [
                self._radius__mm.magnitude * np.cos(phi_circle),
                self._radius__mm.magnitude * np.sin(phi_circle),
            ]
        )
        * pq.mm
    )

assign_mfs2mns

assign_mfs2mns(n_neighbours: int = 3, conf: float = 0.999, n_jobs: int = -2, verbose: bool = True) -> None

Assign muscle fibers to motor neurons using biologically realistic principles.

This method implements an assignment algorithm that balances multiple biological constraints:

  1. Proximity: Fibers closer to innervation centers are more likely to be assigned
  2. Territory size: Each motor unit has a target number of fibers based on its size
  3. Self-avoidance: Neighboring fibers avoid belonging to the same motor unit
  4. Gaussian territories: Fiber territories follow roughly Gaussian distributions

The assignment uses a probabilistic approach where each fiber is assigned based on the posterior probability computed from prior probabilities (target fiber numbers) and likelihoods (spatial clustering with Gaussian territories).

Parameters:

Name Type Description Default
n_neighbours int

Number of neighboring fibers to consider for self-avoiding phenomena. Higher values increase intermingling between motor units but may slow computation. Typical range: 2-5.

3
conf float

Confidence interval that defines the relationship between innervation area and Gaussian distribution variance. Higher values create tighter, more compact territories. Should be between 0.9 and 0.999.

0.999
n_jobs int

Number of parallel workers for out-of-circle coefficient computation.

  • n_jobs=-1: Use all CPU cores
  • n_jobs=-2: Use all cores except one (recommended, keeps system responsive)
  • n_jobs=-3: Use all cores except two
  • n_jobs=1: No parallelization
  • n_jobs=N: Use exactly N cores
-2
verbose bool

If True, display progress bars and status messages. Set to False to disable.

True

Raises:

Type Description
ValueError

If innervation_center_positions is None. Call distribute_innervation_centers() first, or if muscle fiber centers are not available.

Notes

Results are stored in the assignment property after execution.

The algorithm compensates for out-of-muscle effects by calculating how much of each motor unit's Gaussian distribution falls outside the circular muscle boundary and adjusting the in-muscle probabilities accordingly.

The self-avoidance mechanism promotes realistic intermingling by reducing the probability of assigning a fiber to a motor unit if its neighbors are already assigned to that unit.

Source code in myogen/simulator/core/muscle/muscle.py
def assign_mfs2mns(self, n_neighbours: int = 3, conf: float = 0.999, n_jobs: int = -2, verbose: bool = True) -> None:
    """
    Assign muscle fibers to motor neurons using biologically realistic principles.

    This method implements an assignment algorithm that balances
    multiple biological constraints:

    1. Proximity: Fibers closer to innervation centers are more likely to be assigned
    2. Territory size: Each motor unit has a target number of fibers based on its size
    3. Self-avoidance: Neighboring fibers avoid belonging to the same motor unit
    4. Gaussian territories: Fiber territories follow roughly Gaussian distributions

    The assignment uses a probabilistic approach where each fiber is assigned
    based on the posterior probability computed from prior probabilities (target
    fiber numbers) and likelihoods (spatial clustering with Gaussian territories).

    Parameters
    ----------
    n_neighbours : int, default 3
        Number of neighboring fibers to consider for self-avoiding phenomena.
        Higher values increase intermingling between motor units but may slow
        computation. Typical range: 2-5.
    conf : float, default 0.999
        Confidence interval that defines the relationship between innervation
        area and Gaussian distribution variance. Higher values create tighter,
        more compact territories. Should be between 0.9 and 0.999.
    n_jobs : int, default -2
        Number of parallel workers for out-of-circle coefficient computation.

        - n_jobs=-1: Use all CPU cores
        - n_jobs=-2: Use all cores except one (recommended, keeps system responsive)
        - n_jobs=-3: Use all cores except two
        - n_jobs=1: No parallelization
        - n_jobs=N: Use exactly N cores
    verbose : bool, default=True
        If True, display progress bars and status messages. Set to False to disable.

    Raises
    ------
    ValueError
        If innervation_center_positions is None. Call distribute_innervation_centers()
        first, or if muscle fiber centers are not available.

    Notes
    -----
    Results are stored in the `assignment` property after execution.

    The algorithm compensates for out-of-muscle effects by calculating how much
    of each motor unit's Gaussian distribution falls outside the circular muscle
    boundary and adjusting the in-muscle probabilities accordingly.

    The self-avoidance mechanism promotes realistic intermingling by reducing
    the probability of assigning a fiber to a motor unit if its neighbors are
    already assigned to that unit.
    """
    # Ensure innervation_center_positions is available
    if self._innervation_center_positions__mm is None:
        raise ValueError(
            "Innervation center positions not computed. "
            "Call distribute_innervation_centers() first."
        )

    if self._muscle_fiber_centers__mm is None:
        raise ValueError(
            "Muscle fiber centers not computed. Call generate_muscle_fiber_centers() first."
        )

    # Out-of-muscle area compensation
    # Calculates how much of the MU's gaussian distribution is outside of the
    # muscle border and inflates the rest of the distribution according to it
    # Work with magnitude to avoid quantity issues in integration
    radius_magnitude = self._radius__mm.magnitude

    c = chi2.ppf(conf, 2)

    def sigma(ia):
        # ia should be in mm^2, extract magnitude if it's a quantity
        ia_magnitude = ia.magnitude if hasattr(ia, "magnitude") else ia
        return np.eye(2) * ia_magnitude / np.pi / c

    # Helper function for parallel computation of out-of-circle coefficients
    def _compute_out_circle_coeff_single_mu(
        mu_index: int,
        radius_mag: float,
        innervation_center_mean: np.ndarray,
        desired_area,
        c_value: float,
    ) -> float:
        """
        Compute out-of-circle coefficient for a single motor unit.

        This function is designed to be called in parallel for each motor unit.

        Parameters
        ----------
        mu_index : int
            Motor unit index (for reference only, not used in computation).
        radius_mag : float
            Muscle radius magnitude in mm.
        innervation_center_mean : np.ndarray
            Innervation center position [x, y] in mm (magnitude only).
        desired_area : float
            Desired innervation area in mm² (magnitude only).
        c_value : float
            Chi-squared value for confidence interval.

        Returns
        -------
        float
            Out-of-circle coefficient for this motor unit.
        """

        def borderfun_pos(x):
            return np.real(np.sqrt(radius_mag**2 - x**2))

        def borderfun_neg(x):
            return np.real(-np.sqrt(radius_mag**2 - x**2))

        # Compute covariance matrix
        ia_magnitude = (
            desired_area.magnitude if hasattr(desired_area, "magnitude") else desired_area
        )
        cov = np.eye(2) * ia_magnitude / np.pi / c_value

        def probfun(y, x):
            points = (
                np.column_stack([x.ravel(), y.ravel()])
                if hasattr(x, "ravel")
                else np.array([[x, y]])
            )
            return multivariate_normal.pdf(
                points, mean=innervation_center_mean, cov=cov
            ).reshape(np.array(x).shape)

        # Use dblquad for integration (equivalent to MATLAB's integral2)
        in_circle_int = dblquad(
            probfun,
            -radius_mag,
            radius_mag,
            borderfun_neg,
            borderfun_pos,
        )[0]  # dblquad returns (integral, error)

        return 1 / in_circle_int

    # Parallel computation of out-of-circle coefficients
    results = []
    with tqdm(
        total=self._number_of_neurons,
        desc="Calculating out-of-circle coefficients",
        unit="MU",
        disable=not verbose,
    ) as pbar:
        for coeff in Parallel(
            n_jobs=n_jobs,
            return_as="generator",
            verbose=0,
            batch_size="auto",
        )(
            delayed(_compute_out_circle_coeff_single_mu)(
                mu,
                radius_magnitude,
                self._innervation_center_positions__mm[mu].magnitude,
                self.desired_innervation_areas__mm2[mu],
                c,
            )
            for mu in range(self._number_of_neurons)
        ):
            results.append(coeff)
            pbar.update(1)

    out_circle_coeff = np.array(results)

    # Find nearest neighbors for suppression (equivalent to MATLAB's knnsearch)
    # Use magnitudes for sklearn compatibility
    if n_neighbours > 0:
        nbrs = NearestNeighbors(n_neighbors=n_neighbours + 1).fit(
            self._muscle_fiber_centers__mm.magnitude
        )
        _, neighbours = nbrs.kneighbors(self._muscle_fiber_centers__mm.magnitude)
        neighbours = neighbours[:, 1:]  # Exclude self (equivalent to neighbours(:,2:end))

    # Pre-compute constant values for vectorized assignment (optimization)
    # A priori probabilities (constant for all fibers)
    apriori_probs = self.desired_number_of_innervated_fibers / self._number_of_muscle_fibers

    # Pre-compute means and covariances for all motor units
    mu_means = np.array(
        [
            self._innervation_center_positions__mm[mu].magnitude
            for mu in range(self._number_of_neurons)
        ]
    )  # Shape: (n_neurons, 2)

    mu_covs = np.array(
        [
            sigma(self.desired_innervation_areas__mm2[mu])
            for mu in range(self._number_of_neurons)
        ]
    )  # Shape: (n_neurons, 2, 2)

    # Pre-compute inverse covariances and determinants for faster PDF computation
    mu_cov_invs = np.array([np.linalg.pinv(cov) for cov in mu_covs])
    mu_cov_dets = np.array([np.linalg.det(cov) for cov in mu_covs])

    # Assignment procedure
    self._assignment = np.full(self._number_of_muscle_fibers, np.nan)
    randomized_mf = get_random_generator().permutation(self._number_of_muscle_fibers)

    for mf in tqdm(randomized_mf, desc="Assigning muscle fibers to motor neurons", unit="MF", disable=not verbose):
        # Vectorized computation of likelihoods for all motor units
        # Compute differences: (n_neurons, 2)
        diffs = self._muscle_fiber_centers__mm[mf, :].magnitude - mu_means

        # Compute Mahalanobis distances efficiently
        # For each MU: (x - mean)^T * inv(cov) * (x - mean)
        mahal_dists = np.sum(
            np.sum(diffs[:, np.newaxis, :] * mu_cov_invs, axis=2) * diffs, axis=1
        )

        # Compute PDF values for all MUs at once
        clust_hoods = 1.0 / np.sqrt((2 * np.pi) ** 2 * mu_cov_dets) * np.exp(-0.5 * mahal_dists)
        clust_hoods *= out_circle_coeff

        # Compute posterior probabilities (vectorized)
        probs = apriori_probs * clust_hoods

        # Apply neighbor suppression
        if n_neighbours > 0:
            neighbor_assignments = self._assignment[neighbours[mf]]
            for mu in range(self._number_of_neurons):
                if np.any(neighbor_assignments == mu):
                    probs[mu] = 0

        # Normalize probabilities
        prob_sum = np.sum(probs)
        if prob_sum > 0:
            probs = probs / prob_sum
        else:
            # Fallback if all probabilities are zero
            probs = np.ones(self._number_of_neurons) / self._number_of_neurons

        # Sample from the probability distribution (equivalent to MATLAB's randsample)
        self._assignment[mf] = get_random_generator().choice(self._number_of_neurons, p=probs)

    if verbose:
        print(f"Assignment completed. {self._number_of_muscle_fibers} muscle fibers assigned.")

resulting_fiber_assignment

resulting_fiber_assignment(mu: int) -> Quantity__mm

Get the muscle fiber positions assigned to a specific motor unit.

Parameters:

Name Type Description Default
mu int

Motor unit index (0-based). Must be less than the total number of motor units.

required

Returns:

Type Description
Quantity__mm

Array of shape (n_assigned_fibers, 2) containing the [x, y] coordinates (in mm) of all muscle fibers assigned to the specified motor unit. If no fibers are assigned to the motor unit, returns an empty array.

Raises:

Type Description
IndexError

If mu is outside the valid range [0, n_motor_units-1].

ValueError

If the muscle fiber assignment has not been completed yet.

Examples:

>>> fiber_positions = muscle.resulting_fiber_assignment(0)
>>> print(f"Motor unit 0 has {len(fiber_positions)} fibers")
>>> print(f"First fiber position: x={fiber_positions[0,0]:.2f}, y={fiber_positions[0,1]:.2f}")
Notes

This method should only be called after assign_mfs2mns() has been executed. The returned coordinates are in the muscle's coordinate system with the origin at the muscle center.

Source code in myogen/simulator/core/muscle/muscle.py
def resulting_fiber_assignment(self, mu: int) -> Quantity__mm:
    """
    Get the muscle fiber positions assigned to a specific motor unit.

    Parameters
    ----------
    mu : int
        Motor unit index (0-based). Must be less than the total number of motor units.

    Returns
    -------
    Quantity__mm
        Array of shape (n_assigned_fibers, 2) containing the [x, y] coordinates
        (in mm) of all muscle fibers assigned to the specified motor unit.
        If no fibers are assigned to the motor unit, returns an empty array.

    Raises
    ------
    IndexError
        If mu is outside the valid range [0, n_motor_units-1].
    ValueError
        If the muscle fiber assignment has not been completed yet.

    Examples
    --------
    >>> fiber_positions = muscle.resulting_fiber_assignment(0)
    >>> print(f"Motor unit 0 has {len(fiber_positions)} fibers")
    >>> print(f"First fiber position: x={fiber_positions[0,0]:.2f}, y={fiber_positions[0,1]:.2f}")

    Notes
    -----
    This method should only be called after assign_mfs2mns() has been executed.
    The returned coordinates are in the muscle's coordinate system with the
    origin at the muscle center.
    """
    if self._assignment is None:
        raise ValueError(
            "Muscle fiber assignment not completed. "
            "Call assign_mfs2mns() first to assign fibers to motor units."
        )

    if self._muscle_fiber_centers__mm is None:
        raise ValueError(
            "Muscle fiber centers not computed. Call generate_muscle_fiber_centers() first."
        )

    if not (0 <= mu < len(self._recruitment_thresholds)):
        raise IndexError(
            f"Motor unit index {mu} is out of range. "
            f"Valid range is [0, {len(self._recruitment_thresholds) - 1}]."
        )

    return self._muscle_fiber_centers__mm[
        np.where(self._assignment == np.arange(len(self._recruitment_thresholds))[mu])[0]
    ]

HillModel

HillModel(simulation_time__ms: Quantity__ms, time_step__ms: Quantity__ms, muscle_parameters: dict[str, Any], n_motor_units_type1: int, n_motor_units_type2: int, initial_joint_angle__deg: float, initial_muscle_length__L0: float = -1.0, muscle_role: Literal['flexor', 'extensor'] = 'flexor')

API wrapper for the Hill muscle model.

This class provides an intuitive interface for creating Hill muscle models with user-friendly parameter names that are internally mapped to the correct format expected by the underlying Hill implementation.

Parameters:

Name Type Description Default
simulation_time__ms float

Total simulation time in milliseconds

required
time_step__ms float

Integration time step in milliseconds

required
muscle_parameters dict[str, Any]

Dictionary containing Hill muscle model parameters

required
n_motor_units_type1 int

Number of type I motor units

required
n_motor_units_type2 int

Number of type II motor units

required
initial_joint_angle__deg float

Initial joint angle in degrees

required
initial_muscle_length__L0 float

Initial muscle length normalized to L0. If -1, automatically calculated from joint angle. Must be between 0.7 and 1.3 if specified.

-1.0
muscle_role str

Muscle role for antagonist pairs ("flexor" or "extensor"), by default "flexor". Used for joint dynamics calculations and result organization.

'flexor'
Source code in myogen/simulator/neuron/muscle.py
@beartowertype
def __init__(
    self,
    simulation_time__ms: Quantity__ms,
    time_step__ms: Quantity__ms,
    muscle_parameters: dict[str, Any],
    n_motor_units_type1: int,
    n_motor_units_type2: int,
    initial_joint_angle__deg: float,
    initial_muscle_length__L0: float = -1.0,
    muscle_role: Literal["flexor", "extensor"] = "flexor",
):
    # Store original parameters (immutable)
    self.simulation_time__ms = simulation_time__ms
    self.time_step__ms = time_step__ms
    self.muscle_parameters = muscle_parameters.copy()
    self.n_motor_units_type1 = n_motor_units_type1
    self.n_motor_units_type2 = n_motor_units_type2
    self.initial_joint_angle__deg = initial_joint_angle__deg
    self.initial_muscle_length__L0 = initial_muscle_length__L0
    self.muscle_role = muscle_role

    # Private working copies for internal use
    self._simulation_time__ms = simulation_time__ms
    self._time_step__ms = time_step__ms
    self._muscle_parameters = muscle_parameters.copy()
    self._n_motor_units_type1 = n_motor_units_type1
    self._n_motor_units_type2 = n_motor_units_type2
    self._initial_joint_angle__deg = initial_joint_angle__deg
    self._initial_muscle_length__L0 = initial_muscle_length__L0
    self._muscle_role = muscle_role

    # Validate inputs
    self._validate_parameters()

    # Create the underlying Hill model
    self._hill_model = self._create_hill_model()

muscle_length property

muscle_length: ndarray

Get muscle length time series (normalized to L0).

muscle_velocity property

muscle_velocity: ndarray

Get muscle velocity time series (L0/s).

muscle_acceleration property

muscle_acceleration: ndarray

Get muscle acceleration time series (L0/s^2).

muscle_force property

muscle_force: ndarray

Get muscle force time series (normalized to F0).

muscle_torque property

muscle_torque: ndarray

Get muscle torque time series (F0*m).

signed_muscle_torque property

signed_muscle_torque: ndarray

Get muscle torque with correct sign for joint dynamics (F0*m).

type1_activation property

type1_activation: ndarray

Get Type I motor unit activation time series.

type2_activation property

type2_activation: ndarray

Get Type II motor unit activation time series.

motor_unit_forces property

motor_unit_forces: ndarray

Get individual motor unit forces matrix (N_units x time_points).

time_vector property

time_vector: ndarray

Get simulation time vector in milliseconds.

F0 property

F0: float

Get maximum isometric force (F0) in Newtons.

L0 property

L0: float

Get optimal muscle length (L0) in meters.

add_spike

add_spike(motor_unit_id: int, delay__ms: float = 0.0) -> None

Add a spike event for a specific motor unit.

Parameters:

Name Type Description Default
motor_unit_id int

ID of the motor unit (0-based index)

required
delay__ms float

Spike delay in milliseconds, by default 0.0

0.0
Source code in myogen/simulator/neuron/muscle.py
def add_spike(self, motor_unit_id: int, delay__ms: float = 0.0) -> None:
    """
    Add a spike event for a specific motor unit.

    Parameters
    ----------
    motor_unit_id : int
        ID of the motor unit (0-based index)
    delay__ms : float, optional
        Spike delay in milliseconds, by default 0.0
    """
    self._hill_model.addSpike(motor_unit_id, delay__ms)

integrate

integrate(joint_angle__deg: float) -> tuple[float, float, float]

Integrate the muscle model for one time step.

Parameters:

Name Type Description Default
joint_angle__deg float

Current joint angle in degrees

required

Returns:

Type Description
tuple[float, float, float]

Muscle length (normalized to L0), velocity (L0/s), acceleration (L0/s^2)

Source code in myogen/simulator/neuron/muscle.py
def integrate(self, joint_angle__deg: float) -> tuple[float, float, float]:
    """
    Integrate the muscle model for one time step.

    Parameters
    ----------
    joint_angle__deg : float
        Current joint angle in degrees

    Returns
    -------
    tuple[float, float, float]
        Muscle length (normalized to L0), velocity (L0/s), acceleration (L0/s^2)
    """
    return self._hill_model.integrate(np.radians(joint_angle__deg))

create_default_muscle_parameters staticmethod

create_default_muscle_parameters(muscle_type: str = 'FDI') -> dict[str, Any]

Create default muscle parameter dictionary.

Parameters:

Name Type Description Default
muscle_type str

Type of muscle model ("FDI", "Sol"), by default "FDI"

'FDI'

Returns:

Type Description
dict[str, Any]

Dictionary of muscle parameters

Raises:

Type Description
ValueError

If muscle_type is not recognized

Source code in myogen/simulator/neuron/muscle.py
@staticmethod
def create_default_muscle_parameters(muscle_type: str = "FDI") -> dict[str, Any]:
    """
    Create default muscle parameter dictionary.

    Parameters
    ----------
    muscle_type : str, optional
        Type of muscle model ("FDI", "Sol"), by default "FDI"

    Returns
    -------
    dict[str, Any]
        Dictionary of muscle parameters

    Raises
    ------
    ValueError
        If muscle_type is not recognized
    """
    if muscle_type == "FDI":
        return {
            # Muscle geometry
            "alfa0": 0.1606,  # Initial pennation angle [rad]
            "F0": 33.75,  # Maximum isometric force [N]
            "L0": 38.9e-3,  # Optimal fascicle length [m]
            "m": 4.67e-3,  # Muscle mass [kg]
            # Passive elements
            "Kpe": 5,  # Passive elastic element stiffness [F0/L0]
            "b": 0.01,  # Muscle fiber viscous element [F0*s/L0]
            "Em_0": 0.5,  # Normalized muscle deformation
            # Tendon parameters
            "LT_0": 49e-3,  # Tendon length for max isometric force [m]
            "Kse": 27.8,  # Tendon elastic element [F0/LT_0]
            "cT": 0.0047,  # Toe region coefficient
            "LT_r": 0.964,  # Linear region start [LT_0]
            # Force-Length curve parameters (Type I fibers)
            # F_L = exp(-|((L^b - 1)/o)|^r) where L is normalized length
            "b1": 2.3,  # Shape parameter for Type I length-force curve
            "o1": 1.12,  # Width parameter for Type I length-force curve
            "r1": 1.62,  # Asymmetry parameter for Type I length-force curve
            # Force-Length curve parameters (Type II fibers)
            "b2": 1.55,  # Shape parameter for Type II length-force curve
            "o2": 0.75,  # Width parameter for Type II length-force curve
            "r2": 2.12,  # Asymmetry parameter for Type II length-force curve
            # Force-Velocity curve parameters (Type I fibers)
            # For concentric: F_V = (bv - V*(av0 + av1*L + av2*L²))/(bv + V)
            # For eccentric: F_V = (Vmax - V)/(Vmax + V*(cv0 + cv1*L))
            "Vmax1": -7.88,  # Maximum shortening velocity for Type I [L0/s]
            "av01": -4.7,  # Concentric velocity coefficient a0 for Type I
            "av11": 8.41,  # Concentric velocity coefficient a1 for Type I (length-dependent)
            "av21": -5.34,  # Concentric velocity coefficient a2 for Type I (length²-dependent)
            "bv1": 0.35,  # Concentric force-velocity scaling for Type I
            "cv01": 5.88,  # Eccentric velocity coefficient c0 for Type I
            "cv11": 0,  # Eccentric velocity coefficient c1 for Type I (length-dependent)
            # Force-Velocity curve parameters (Type II fibers)
            "Vmax2": -9.15,  # Maximum shortening velocity for Type II [L0/s]
            "av02": -1.53,  # Concentric velocity coefficient a0 for Type II
            "av12": 0,  # Concentric velocity coefficient a1 for Type II
            "av22": 0,  # Concentric velocity coefficient a2 for Type II
            "bv2": 0.69,  # Concentric force-velocity scaling for Type II
            "cv02": 5.7,  # Eccentric velocity coefficient c0 for Type II
            "cv12": 9.18,  # Eccentric velocity coefficient c1 for Type II
            # Muscle-tendon length and moment arm coefficients
            "Ak": [
                85.199931e-3,
                -1.184782e-4,
                -4.6264098e-7,
                9.416143e-10,
                4.854117e-12,
            ],
            "Bk": [6.82847e-3, 4.8396e-5, 3.6942e-8, 6.3113e-10, -6.35837e-11],
            # Motor unit parameters
            "RP": 130,  # Range of twitch force amplitude
            "fP": 3,  # First peak twitch force [mN]
            "RT": 3,  # Range of contraction time
            "durType": 1,  # Distribution type (1=exponential)
            "Tl": 90,  # Longest twitch duration [ms]
            "fsatf": 50,  # First MU saturation frequency [Hz]
            "lsatf": 100,  # Last MU saturation frequency [Hz]
            "satType": 1,  # Saturation type (1=linear)
        }

    elif muscle_type == "Sol":
        return {
            # Soleus muscle parameters (larger, stronger muscle)
            "alfa0": 0.494,
            "F0": 3586,
            "L0": 49e-3,
            "m": 0.526,
            "Kpe": 5,
            "b": 0.005,
            "Em_0": 0.5,
            "LT_0": 0.289,
            "Kse": 27.8,
            "cT": 0.0047,
            "LT_r": 0.964,
            "b1": 2.3,
            "o1": 1.12,
            "r1": 1.62,
            "b2": 1.55,
            "o2": 0.75,
            "r2": 2.12,
            "Vmax1": -7.88,
            "av01": -4.7,
            "av11": 8.41,
            "av21": -5.34,
            "bv1": 0.35,
            "cv01": 5.88,
            "cv11": 0,
            "Vmax2": -9.15,
            "av02": -1.53,
            "av12": 0,
            "av22": 0,
            "bv2": 0.69,
            "cv02": 5.7,
            "cv12": 9.18,
            "Ak": [0.323, 7.219e-4, -2.243e-6, -3.148e-8, 9.274e-11],
            "Bk": [-0.041, 2.574e-4, 5.451e-6, -2.219e-8, -5.494e-11],
            "RP": 130,
            "fP": 3,
            "RT": 3,
            "durType": 1,
            "Tl": 90,
            "fsatf": 50,
            "lsatf": 100,
            "satType": 1,
        }

    else:
        raise ValueError(f"Unknown muscle type: {muscle_type}. Use 'FDI' or 'Sol'.")

ForceModel

ForceModel(recruitment_thresholds: RECRUITMENT_THRESHOLDS__ARRAY, recording_frequency__Hz: Quantity__Hz, longest_duration_rise_time__ms: Quantity__ms = 90.0 * ms, contraction_time_range_factor: float = 3.0)

Force model based on Fuglevand et al. (1993) [1].

This class implements the Fuglevand force generation model for motor unit pools, computing individual motor unit twitch responses and their nonlinear gain modulation based on discharge rates. The model generates realistic force outputs from spike trains using physiologically-based parameters.

Parameters:

Name Type Description Default
recruitment_thresholds RECRUITMENT_THRESHOLDS__ARRAY

Recruitment thresholds for each motor unit. Array of values typically ranging from 0 to 1 where larger motor units have higher thresholds.

required
recording_frequency__Hz Quantity__Hz

Recording frequency in Hz. Determines temporal resolution of force calculations. Typical values: 100-1000 Hz.

required
longest_duration_rise_time__ms Quantity__ms

Longest duration of the rise time in milliseconds. This parameter (T_L in _[1]) determines the contraction time of the slowest motor unit. Typical range: 50-150 ms.

90.0 * pq.ms
contraction_time_range_factor float

Contraction time range factor (RT in _[1]). Determines the spread of contraction times across motor units. Generally between 2 and 5. Higher values create larger differences between fast and slow motor units.

3.0

Attributes:

Name Type Description
peak_twitch_forces__unitless ndarray

Peak twitch forces for each motor unit (unitless). Available after initialization. Computed according to Fuglevand equation 13.

contraction_times__samples ndarray

Contraction times for each motor unit in samples. Available after initialization. Computed according to Fuglevand equation 14.

twitch_mat ndarray

Complete twitch matrix for all motor units. Available after initialization. Shape: (max_twitch_length, n_motor_units).

twitch_list list[ndarray]

List of individual twitch responses for each motor unit. Available after initialization. Each element contains the twitch response for one motor unit.

Raises:

Type Description
ValueError

If recruitment_thresholds is empty or contains invalid values. If recording_frequency__Hz is not positive. If longest_duration_rise_time__ms is not positive. If contraction_time_range_factor is not greater than 1.

References

[1] Fuglevand, A. J., Winter, D. A., & Patla, A. E. (1993). Models of recruitment and rate coding in motor-unit pools. Journal of Neurophysiology, 70(2), 782-797.

Examples:

>>> import numpy as np
>>> from myogen.simulator.core.force import ForceModel
>>> thresholds = np.linspace(0.1, 1.0, 10)
>>> force_model = ForceModel(
...     recruitment_thresholds=thresholds,
...     recording_frequency__Hz=2000.0
... )
Source code in myogen/simulator/core/force/force_model.py
def __init__(
    self,
    recruitment_thresholds: RECRUITMENT_THRESHOLDS__ARRAY,
    recording_frequency__Hz: Quantity__Hz,
    longest_duration_rise_time__ms: Quantity__ms = 90.0 * pq.ms,
    contraction_time_range_factor: float = 3.0,
) -> None:
    # Input validation
    if len(recruitment_thresholds) == 0:
        raise ValueError(
            "recruitment_thresholds cannot be empty. "
            "Please provide at least one recruitment threshold value."
        )

    if not np.all(recruitment_thresholds > 0):
        raise ValueError(
            "All recruitment thresholds must be positive. "
            "Found values: min={:.3f}, max={:.3f}. "
            "Recruitment thresholds typically range from 0.01 to 1.0.".format(
                np.min(recruitment_thresholds), np.max(recruitment_thresholds)
            )
        )

    if recording_frequency__Hz <= 0:
        raise ValueError(
            f"recording_frequency__Hz must be positive, got {recording_frequency__Hz}. "
            "Typical values for EMG/force recordings are between 1000-10000 Hz."
        )

    if longest_duration_rise_time__ms <= 0:
        raise ValueError(
            f"longest_duration_rise_time__ms must be positive, got {longest_duration_rise_time__ms}. "
            "Typical values range from 50-150 ms for human motor units."
        )

    if contraction_time_range_factor <= 1.0:
        raise ValueError(
            f"contraction_time_range_factor must be greater than 1.0, got {contraction_time_range_factor}. "
            "This parameter determines the spread of contraction times. Typical values are 2.0-5.0."
        )

    # Immutable public access
    self.recruitment_thresholds = recruitment_thresholds
    self.recording_frequency__Hz = recording_frequency__Hz
    self.longest_duration_rise_time__ms = longest_duration_rise_time__ms
    self.contraction_time_range_factor = contraction_time_range_factor

    # Private copies for internal modifications
    self._recruitment_thresholds = recruitment_thresholds.copy()
    self._recording_frequency__Hz = recording_frequency__Hz
    self._longest_duration_rise_time__ms = longest_duration_rise_time__ms
    self._contraction_time_range_factor = contraction_time_range_factor

    # Derived properties
    self._number_of_neurons = len(self._recruitment_thresholds)
    self._recruitment_ratio = (
        self._recruitment_thresholds[-1] / self._recruitment_thresholds[0]
    )  # referred in [1] as RP

    self._longest_duration_rise_time__samples = float(
        (
            self._longest_duration_rise_time__ms.rescale("s") * self._recording_frequency__Hz
        ).magnitude
    )  # referred in [1] as T_L (see eq. 14)

    # Simulation results - stored privately, accessed via properties
    self._peak_twitch_forces__unitless: Optional[np.ndarray] = None
    self._contraction_times__samples: Optional[np.ndarray] = None
    self._twitch_mat: Optional[np.ndarray] = None
    self._twitch_list: Optional[list[np.ndarray]] = None

    # Initialize model parameters
    self._compute_twitch_parameters()

peak_twitch_forces__unitless property

peak_twitch_forces__unitless: ndarray

Peak twitch forces for each motor unit (unitless).

Returns:

Type Description
ndarray

Array of peak twitch forces according to Fuglevand model equation 13.

Raises:

Type Description
ValueError

If twitch parameters have not been computed yet.

contraction_times__samples property

contraction_times__samples: ndarray

Contraction times for each motor unit in samples.

Returns:

Type Description
ndarray

Array of contraction times according to Fuglevand model equation 14.

Raises:

Type Description
ValueError

If twitch parameters have not been computed yet.

twitch_mat property

twitch_mat: ndarray

Complete twitch matrix for all motor units.

Returns:

Type Description
ndarray

Matrix of shape (max_twitch_length, n_motor_units) containing the twitch responses for each motor unit.

Raises:

Type Description
ValueError

If twitches have not been initialized yet.

twitch_list property

twitch_list: list[ndarray]

List of individual twitch responses for each motor unit.

Returns:

Type Description
list[ndarray]

List where each element is the twitch response array for one motor unit. Each array may have different lengths based on the motor unit's contraction time.

Raises:

Type Description
ValueError

If twitches have not been initialized yet.

generate_force

generate_force(spike_train__Block: SPIKE_TRAIN__Block, verbose: bool = True) -> FORCE__AnalogSignal

Generate force output from motor unit spike trains using the Fuglevand model.

This method simulates muscle force by converting spike trains into force output through individual motor unit twitches with nonlinear gain modulation based on discharge rate. Each motor unit contributes to the total force according to its twitch properties and firing pattern. The output is resampled to match the recording_frequency__Hz parameter.

Parameters:

Name Type Description Default
spike_train__Block SPIKE_TRAIN__Block

Spike train block containing spike train data for multiple motor neuron pools.

required
verbose bool

If True, display progress bars during force generation. Set to False to disable.

True

Returns:

Type Description
FORCE__AnalogSignal

Force output neo.AnalogSignal representing muscle force over time. Each channel corresponds to one motor neuron pool's force response. Sampling rate matches the recording_frequency__Hz parameter.

Raises:

Type Description
ValueError

If spike train matrix dimensions don't match the number of motor units. If twitch parameters have not been computed.

Notes

The force generation follows these steps: 1. Convert spike trains to inter-pulse intervals (IPIs) 2. Calculate nonlinear gain based on discharge rates 3. Sum weighted twitch responses for each spike 4. Apply gain modulation to final force output 5. Resample output to match recording_frequency__Hz

Source code in myogen/simulator/core/force/force_model.py
def generate_force(self, spike_train__Block: SPIKE_TRAIN__Block, verbose: bool = True) -> FORCE__AnalogSignal:
    """
    Generate force output from motor unit spike trains using the Fuglevand model.

    This method simulates muscle force by converting spike trains into force output
    through individual motor unit twitches with nonlinear gain modulation based on
    discharge rate. Each motor unit contributes to the total force according to its
    twitch properties and firing pattern. The output is resampled to match the
    recording_frequency__Hz parameter.

    Parameters
    ----------
    spike_train__Block : SPIKE_TRAIN__Block
        Spike train block containing spike train data for multiple motor neuron pools.
    verbose : bool, default=True
        If True, display progress bars during force generation. Set to False to disable.

    Returns
    -------
    FORCE__AnalogSignal
        Force output neo.AnalogSignal representing muscle force over time.
        Each channel corresponds to one motor neuron pool's force response.
        Sampling rate matches the recording_frequency__Hz parameter.

    Raises
    ------
    ValueError
        If spike train matrix dimensions don't match the number of motor units.
        If twitch parameters have not been computed.

    Notes
    -----
    The force generation follows these steps:
    1. Convert spike trains to inter-pulse intervals (IPIs)
    2. Calculate nonlinear gain based on discharge rates
    3. Sum weighted twitch responses for each spike
    4. Apply gain modulation to final force output
    5. Resample output to match recording_frequency__Hz
    """
    if self._twitch_list is None:
        raise ValueError(
            "Twitch parameters not available. "
            "This should not occur if the model was properly initialized. "
            "Please reinitialize the ForceModel."
        )

    # Extract timing information from spike trains
    spiketrain_timestep__ms = float(
        spike_train__Block.segments[0].spiketrains[0].sampling_period.rescale("ms").magnitude
    )

    forces = []
    for i, segment in enumerate(spike_train__Block.segments):
        if len(segment.spiketrains) != self._number_of_neurons:
            raise ValueError(
                f"MU pool {i} has {len(segment.spiketrains)} neurons, "
                f"but force model was initialized with {self._number_of_neurons} motor units. "
                "The number of neurons in the spike train neo.Block must match the number of recruitment thresholds."
            )

        spike_array = bin_spike_trains(
            segment.spiketrains,
            bin_size=segment.spiketrains[0].sampling_period,
            t_start=segment.t_start,
            t_stop=segment.t_stop,
            sparse=True,
        ).T

        # Generate force with resampling handled internally
        force_output = self._generate_force(
            spike_array, spiketrain_timestep__ms, prefix=f"Pool {i + 1}", verbose=verbose
        )
        forces.append(force_output)

    return AnalogSignal(
        np.stack(forces, axis=-1) * pq.dimensionless,
        t_start=spike_train__Block.segments[0].t_start.rescale("s"),
        sampling_rate=self._recording_frequency__Hz,
    )

ForceModelVectorized

ForceModelVectorized(recruitment_thresholds: RECRUITMENT_THRESHOLDS__ARRAY, recording_frequency__Hz: Quantity__Hz, longest_duration_rise_time__ms: Quantity__ms = 90.0 * ms, contraction_time_range_factor: float = 3.0)

Vectorized force model based on Fuglevand et al. (1993) [1].

This is an optimized version of ForceModel that uses numpy vectorization for significantly better performance, especially for long simulations. It shares the IPI/gain/twitch pipeline with the reference implementation via force_utils so the output is guaranteed to match the reference model bit-for-bit (modulo per-spike accumulation order).

Parameters:

Name Type Description Default
recruitment_thresholds RECRUITMENT_THRESHOLDS__ARRAY

Recruitment thresholds for each motor unit.

required
recording_frequency__Hz Quantity__Hz

Recording frequency in Hz. Determines temporal resolution of force calculations. Typical values: 100-1000 Hz.

required
longest_duration_rise_time__ms Quantity__ms

Longest duration of the rise time in milliseconds.

90.0 * pq.ms
contraction_time_range_factor float

Contraction time range factor. Determines the spread of contraction times across motor units. Generally between 2 and 5.

3.0
References

[1] Fuglevand, A. J., Winter, D. A., & Patla, A. E. (1993). Models of recruitment and rate coding in motor-unit pools. Journal of Neurophysiology, 70(2), 782-797.

Source code in myogen/simulator/core/force/force_model_vectorized.py
def __init__(
    self,
    recruitment_thresholds: RECRUITMENT_THRESHOLDS__ARRAY,
    recording_frequency__Hz: Quantity__Hz,
    longest_duration_rise_time__ms: Quantity__ms = 90.0 * pq.ms,
    contraction_time_range_factor: float = 3.0,
) -> None:
    # Input validation
    if len(recruitment_thresholds) == 0:
        raise ValueError(
            "recruitment_thresholds cannot be empty. "
            "Please provide at least one recruitment threshold value."
        )

    if not np.all(recruitment_thresholds > 0):
        raise ValueError(
            "All recruitment thresholds must be positive. "
            "Found values: min={:.3f}, max={:.3f}. "
            "Recruitment thresholds typically range from 0.01 to 1.0.".format(
                np.min(recruitment_thresholds), np.max(recruitment_thresholds)
            )
        )

    if recording_frequency__Hz <= 0:
        raise ValueError(
            f"recording_frequency__Hz must be positive, got {recording_frequency__Hz}. "
            "Typical values for EMG/force recordings are between 1000-10000 Hz."
        )

    if longest_duration_rise_time__ms <= 0:
        raise ValueError(
            f"longest_duration_rise_time__ms must be positive, got {longest_duration_rise_time__ms}. "
            "Typical values range from 50-150 ms for human motor units."
        )

    if contraction_time_range_factor <= 1.0:
        raise ValueError(
            f"contraction_time_range_factor must be greater than 1.0, got {contraction_time_range_factor}. "
            "This parameter determines the spread of contraction times. Typical values are 2.0-5.0."
        )

    # Immutable public access
    self.recruitment_thresholds = recruitment_thresholds
    self.recording_frequency__Hz = recording_frequency__Hz
    self.longest_duration_rise_time__ms = longest_duration_rise_time__ms
    self.contraction_time_range_factor = contraction_time_range_factor

    # Private copies for internal modifications
    self._recruitment_thresholds = recruitment_thresholds.copy()
    self._recording_frequency__Hz = recording_frequency__Hz
    self._longest_duration_rise_time__ms = longest_duration_rise_time__ms
    self._contraction_time_range_factor = contraction_time_range_factor

    # Derived properties
    self._number_of_neurons = len(self._recruitment_thresholds)
    self._recruitment_ratio = (
        self._recruitment_thresholds[-1] / self._recruitment_thresholds[0]
    )

    # Match ForceModel's quantity-aware sample conversion exactly.
    self._longest_duration_rise_time__samples = float(
        (
            self._longest_duration_rise_time__ms.rescale("s")
            * self._recording_frequency__Hz
        ).magnitude
    )

    # Simulation results
    self._peak_twitch_forces__unitless: Optional[np.ndarray] = None
    self._contraction_times__samples: Optional[np.ndarray] = None
    self._twitch_mat: Optional[np.ndarray] = None
    self._twitch_list: Optional[list[np.ndarray]] = None

    # Initialize model parameters
    self._compute_twitch_parameters()

generate_force

generate_force(spike_train__Block: SPIKE_TRAIN__Block, verbose: bool = True) -> FORCE__AnalogSignal

Generate force output from motor unit spike trains.

The body mirrors ForceModel.generate_force so that the two implementations cannot drift apart silently. Only the per-spike accumulation differs (vectorized vs. per-spike loop).

Source code in myogen/simulator/core/force/force_model_vectorized.py
def generate_force(
    self, spike_train__Block: SPIKE_TRAIN__Block, verbose: bool = True
) -> FORCE__AnalogSignal:
    """
    Generate force output from motor unit spike trains.

    The body mirrors `ForceModel.generate_force` so that the two
    implementations cannot drift apart silently. Only the per-spike
    accumulation differs (vectorized vs. per-spike loop).
    """
    if self._twitch_list is None:
        raise ValueError(
            "Twitch parameters not available. "
            "This should not occur if the model was properly initialized."
        )

    # Extract timing information
    spiketrain_timestep__ms = float(
        spike_train__Block.segments[0]
        .spiketrains[0]
        .sampling_period.rescale("ms")
        .magnitude
    )

    forces = []
    for i, segment in enumerate(spike_train__Block.segments):
        if len(segment.spiketrains) != self._number_of_neurons:
            raise ValueError(
                f"MU pool {i} has {len(segment.spiketrains)} neurons, "
                f"but force model was initialized with {self._number_of_neurons} motor units."
            )

        spike_array = bin_spike_trains(
            segment.spiketrains,
            bin_size=segment.spiketrains[0].sampling_period,
            t_start=segment.t_start,
            t_stop=segment.t_stop,
            sparse=True,
        ).T

        # Generate force with vectorized implementation
        force_output = self._generate_force_vectorized(
            spike_array,
            spiketrain_timestep__ms,
            prefix=f"Pool {i + 1}",
            verbose=verbose,
        )
        forces.append(force_output)

    return AnalogSignal(
        np.stack(forces, axis=-1) * pq.dimensionless,
        t_start=spike_train__Block.segments[0].t_start.rescale("s"),
        sampling_rate=self._recording_frequency__Hz,
    )

EMG

SurfaceEMG

SurfaceEMG(muscle_model: Muscle, electrode_arrays: list[SurfaceElectrodeArray], sampling_frequency__Hz: Quantity__Hz = 2048.0 * Hz, sampling_points_in_t_and_z_domains: int = 256, sampling_points_in_theta_domain: int = 32, MUs_to_simulate: list[int] | None = None, internal_sampling_frequency__Hz: Quantity__Hz | None = None, iap_kernel_length__mm: float | None = None, use_unified: bool = False)

Surface Electromyography (sEMG) Simulation.

This class provides a simulation framework for generating surface electromyography signals from the muscle. It implements the multi-layered cylindrical volume conductor model from Farina et al. 2004 [1].

Parameters:

Name Type Description Default
muscle_model Muscle

Pre-computed muscle model (see myogen.simulator.Muscle).

required
electrode_arrays list[SurfaceElectrodeArray]

List of electrode arrays to use for simulation (see myogen.simulator.SurfaceElectrodeArray).

required
sampling_frequency__Hz float

Sampling frequency in Hz. Default is set to 2048 Hz as used by the Quattrocento (OT Bioelettronica, Turin, Italy) system.

2048.0
sampling_points_in_t_and_z_domains int

Spatial and temporal discretization resolution for numerical integration. Controls the accuracy of the volume conductor calculations but significantly impacts computational cost (scales quadratically). Higher values provide better numerical accuracy at the expense of simulation time. Default is set to 256 samples.

256
sampling_points_in_theta_domain int

Angular discretization for cylindrical coordinate system in degrees. Higher values provide better spatial resolution but cause numerical overflow in Bessel functions. Default is set to 32 points to avoid numerical instability. WARNING: Values >64 cause extreme Bessel function overflow leading to incorrect results. This is suitable for most EMG studies.

32
MUs_to_simulate list[int]

Indices of motor units to simulate. If None, all motor units are simulated. Default is None. For computational efficiency, consider simulating subsets for initial analysis. Indices correspond to the recruitment order (0 is recruited first).

None
internal_sampling_frequency__Hz Quantity__Hz

Internal sampling frequency for MUAP computation before downsampling. If None, defaults to 10 kHz. Higher values provide better temporal resolution but increase computation time. Default is 10 kHz.

None
iap_kernel_length__mm float

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.

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.

None

Attributes:

Name Type Description
muaps__Block SURFACE_MUAP__Block

Motor Unit Action Potential (MUAP) templates for each electrode array as a neo.Block. Available after simulate_muaps().

surface_emg__Block SURFACE_EMG__Block

Surface EMG signals for each electrode array as a neo.Block. Available after simulate_surface_emg().

noisy_surface_emg__Block SURFACE_EMG__Block

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().

References

[1] Farina, D., Mesin, L., Martina, S., Merletti, R., 2004. A surface EMG generation model with multilayer cylindrical description of the volume conductor. IEEE Transactions on Biomedical Engineering 51, 415–426. https://doi.org/10.1109/TBME.2003.820998

Source code in myogen/simulator/core/emg/surface/surface_emg.py
def __init__(
    self,
    muscle_model: Muscle,
    electrode_arrays: list[SurfaceElectrodeArray],
    sampling_frequency__Hz: Quantity__Hz = 2048.0 * pq.Hz,
    sampling_points_in_t_and_z_domains: int = 256,
    sampling_points_in_theta_domain: int = 32,
    MUs_to_simulate: list[int] | None = None,
    internal_sampling_frequency__Hz: Quantity__Hz | None = None,
    iap_kernel_length__mm: float | None = None,
    use_unified: bool = False,
):
    # Immutable public arguments - never modify these
    self.muscle_model = muscle_model
    self.electrode_arrays = electrode_arrays
    self.sampling_frequency__Hz = sampling_frequency__Hz
    self.sampling_points_in_t_and_z_domains = sampling_points_in_t_and_z_domains
    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._use_unified = use_unified

    # Internal sampling frequency for higher resolution MUAP computation
    # If not specified, defaults to 10 kHz for better MUAP resolution
    if internal_sampling_frequency__Hz is None:
        internal_sampling_frequency__Hz = 10000.0 * pq.Hz
    self.internal_sampling_frequency__Hz = internal_sampling_frequency__Hz

    # Private copies for internal modifications (extract magnitudes)
    self._muscle_model = muscle_model
    self._electrode_arrays = electrode_arrays
    self._sampling_frequency__Hz = float(sampling_frequency__Hz.rescale(pq.Hz).magnitude)
    self._internal_sampling_frequency__Hz = float(
        internal_sampling_frequency__Hz.rescale(pq.Hz).magnitude
    )

    # Calculate upsampling factor and internal sample count
    self._upsampling_factor = (
        self._internal_sampling_frequency__Hz / self._sampling_frequency__Hz
    )
    self._internal_sampling_points = int(
        np.round(sampling_points_in_t_and_z_domains * self._upsampling_factor)
    )

    self._sampling_points_in_t_and_z_domains = sampling_points_in_t_and_z_domains
    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

    # Derived properties from muscle model - immutable public access
    self.mean_conduction_velocity__m_s = self._muscle_model.mean_conduction_velocity__m_s
    self.mean_fiber_length__mm = self._muscle_model.mean_fiber_length__mm
    self.var_fiber_length__mm = self._muscle_model.var_fiber_length__mm
    self.radius_bone__mm = self._muscle_model.radius_bone__mm
    self.fat_thickness__mm = self._muscle_model.fat_thickness__mm
    self.skin_thickness__mm = self._muscle_model.skin_thickness__mm
    self.muscle_conductivity_radial__S_m = self._muscle_model.muscle_conductivity_radial__S_m
    self.muscle_conductivity_longitudinal__S_m = (
        self._muscle_model.muscle_conductivity_longitudinal__S_m
    )
    self.fat_conductivity__S_m = self._muscle_model.fat_conductivity__S_m
    self.skin_conductivity__S_m = self._muscle_model.skin_conductivity__S_m

    # Private copies for internal modifications (extract magnitudes if Quantity objects)
    def _extract_value(val):
        """Helper to extract magnitude from Quantity or return float directly."""
        if hasattr(val, "magnitude"):
            return float(val.magnitude)
        return float(val)

    self._mean_conduction_velocity__m_s = _extract_value(
        self._muscle_model.mean_conduction_velocity__m_s
    )
    self._mean_fiber_length__mm = _extract_value(self._muscle_model.mean_fiber_length__mm)
    self._var_fiber_length__mm = _extract_value(self._muscle_model.var_fiber_length__mm)
    self._radius_bone__mm = _extract_value(self._muscle_model.radius_bone__mm)
    self._fat_thickness__mm = _extract_value(self._muscle_model.fat_thickness__mm)
    self._skin_thickness__mm = _extract_value(self._muscle_model.skin_thickness__mm)
    self._muscle_conductivity_radial__S_m = _extract_value(
        self._muscle_model.muscle_conductivity_radial__S_m
    )
    self._muscle_conductivity_longitudinal__S_m = _extract_value(
        self._muscle_model.muscle_conductivity_longitudinal__S_m
    )
    self._fat_conductivity__S_m = _extract_value(self._muscle_model.fat_conductivity__S_m)
    self._skin_conductivity__S_m = _extract_value(self._muscle_model.skin_conductivity__S_m)

    # Calculate total radius - immutable and private
    self._radius_muscle__mm = _extract_value(self._muscle_model.radius__mm)
    self.radius_total = (
        self._radius_muscle__mm + self._fat_thickness__mm + self._skin_thickness__mm
    )
    self._radius_total = self.radius_total

    # Simulation results - stored privately, accessed via properties
    self._muaps__Block: Optional[SURFACE_MUAP__Block] = None
    self._surface_emg__Block: Optional[SURFACE_EMG__Block] = None
    self._noisy_surface_emg__Block: Optional[SURFACE_EMG__Block] = None
    self._spike_train__Block: Optional[SPIKE_TRAIN__Block] = None

muaps__Block property

muaps__Block: SURFACE_MUAP__Block

Motor Unit Action Potential (MUAP) templates for each electrode array.

Returns:

Type Description
list[SURFACE_MUAP_SHAPE__TENSOR]

List of MUAP templates for each electrode array.

Raises:

Type Description
ValueError

If MUAP templates have not been computed yet.

surface_emg__Block property

surface_emg__Block: SURFACE_EMG__Block

Surface EMG signals for each electrode array stored in a neo.Block.

Returns:

Type Description
SURFACE_EMG__Block

Surface EMG signals for each electrode array stored in a neo.Block.

Raises:

Type Description
ValueError

If surface EMG has not been computed yet.

noisy_surface_emg__Block property

noisy_surface_emg__Block: SURFACE_EMG__Block

Noisy surface EMG signals for each electrode array.

Returns:

Type Description
SURFACE_EMG__Block

Noisy surface EMG signals for each electrode array.

Raises:

Type Description
ValueError

If noisy surface EMG has not been computed yet.

spike_train__Block property

spike_train__Block: SPIKE_TRAIN__Block

Spike train block used for EMG generation.

Returns:

Type Description
SPIKE_TRAIN__Block

The spike train block used in the simulation.

Raises:

Type Description
ValueError

If spike train block has not been set yet.

simulate_muaps

simulate_muaps(n_jobs: int = -2, verbose: bool = True, use_gpu: Optional[bool] = None) -> SURFACE_MUAP__Block

Simulate MUAPs for all electrode arrays using the provided muscle model.

This method generates Motor Unit Action Potential (MUAP) templates that represent the electrical signature of each motor unit as recorded by the surface electrodes. The simulation uses the multi-layered cylindrical volume conductor model with parallel processing for improved performance.

Parameters:

Name Type Description Default
n_jobs int

Number of parallel workers for motor unit processing. Default is -2. - n_jobs=-1: Use all CPU cores - n_jobs=-2: Use all cores except one (recommended, keeps system responsive) - n_jobs=-3: Use all cores except two - n_jobs=1: No parallelization - n_jobs=N: Use exactly N cores

-2
verbose bool

If True, display progress bars. Set to False to disable.

True
use_gpu bool or None

GPU acceleration control for the per-fiber volume-conductor solve (mirrors myogen.simulator.MotorUnitSim.calc_sfaps):

  • None → auto: use GPU if CuPy is available and MYOGEN_DISABLE_GPU is not set in the environment; CPU otherwise.
  • True → require GPU; raises RuntimeError if unavailable.
  • False → force CPU execution.

Note: CuPy only supports NVIDIA GPUs (CUDA). AMD/ROCm is not supported. For typical surface-EMG problem sizes, CPU is often faster than GPU due to host↔device transfer overhead.

None

Returns:

Type Description
SURFACE_MUAP__Block

neo.Block of generated MUAP templates for each electrode array. Results are stored in the muaps property after execution.

Notes

This method must be called before simulate_surface_emg(). The generated MUAP templates are used as basis functions for EMG signal synthesis.

The motor units are processed in parallel using joblib, with each motor unit's fibers processed sequentially to maintain optimization efficiency.

Source code in myogen/simulator/core/emg/surface/surface_emg.py
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
def simulate_muaps(
    self,
    n_jobs: int = -2,
    verbose: bool = True,
    use_gpu: Optional[bool] = None,
) -> SURFACE_MUAP__Block:
    """
    Simulate MUAPs for all electrode arrays using the provided muscle model.

    This method generates Motor Unit Action Potential (MUAP) templates that represent
    the electrical signature of each motor unit as recorded by the surface electrodes.
    The simulation uses the multi-layered cylindrical volume conductor model with
    parallel processing for improved performance.

    Parameters
    ----------
    n_jobs : int, optional
        Number of parallel workers for motor unit processing. Default is -2.
        - n_jobs=-1: Use all CPU cores
        - n_jobs=-2: Use all cores except one (recommended, keeps system responsive)
        - n_jobs=-3: Use all cores except two
        - n_jobs=1: No parallelization
        - n_jobs=N: Use exactly N cores
    verbose : bool, default=True
        If True, display progress bars. Set to False to disable.
    use_gpu : bool or None, default=None
        GPU acceleration control for the per-fiber volume-conductor solve
        (mirrors `myogen.simulator.MotorUnitSim.calc_sfaps`):

        - ``None``  → auto: use GPU if CuPy is available and
          ``MYOGEN_DISABLE_GPU`` is not set in the environment; CPU otherwise.
        - ``True``  → require GPU; raises ``RuntimeError`` if unavailable.
        - ``False`` → force CPU execution.

        Note: CuPy only supports NVIDIA GPUs (CUDA). AMD/ROCm is not
        supported. For typical surface-EMG problem sizes, CPU is often
        faster than GPU due to host↔device transfer overhead.

    Returns
    -------
    SURFACE_MUAP__Block
        neo.Block of generated MUAP templates for each electrode array.
        Results are stored in the `muaps` property after execution.

    Notes
    -----
    This method must be called before simulate_surface_emg(). The generated MUAP
    templates are used as basis functions for EMG signal synthesis.

    The motor units are processed in parallel using joblib, with each motor unit's
    fibers processed sequentially to maintain optimization efficiency.
    """
    # Set default MUs to simulate
    if self._MUs_to_simulate is None:
        self._MUs_to_simulate = list(
            range(len(self._muscle_model.resulting_number_of_innervated_fibers))
        )

    # Calculate innervation zone variance
    innervation_zone_variance = (
        self._mean_fiber_length__mm * 0.1
    )  # 10% of the mean fiber length (see Botelho et al. 2019 [6])

    # Extract fiber counts
    number_of_fibers_per_MUs = self._muscle_model.resulting_number_of_innervated_fibers

    # Create time array at INTERNAL sampling frequency for higher resolution
    t_internal = np.linspace(
        0,
        (self._internal_sampling_points - 1) / self._internal_sampling_frequency__Hz * 1e-3,
        self._internal_sampling_points,
    )

    # Get total number of motor units
    n_motor_units = len(number_of_fibers_per_MUs)

    # Pre-calculate innervation zones for all MUs
    innervation_zones = get_random_generator().uniform(
        low=-innervation_zone_variance / 2,
        high=innervation_zone_variance / 2,
        size=n_motor_units,
    )

    # Pre-calculate per-MU fiber-length variations in the PARENT process.
    # Drawing inside `_process_single_mu` (a joblib worker) is unsafe under
    # the default loky backend: each worker inherits a copy of the parent's
    # RNG state via fork/pickle and then advances independently, so two MUs
    # processed in two different workers would draw from identical state.
    # Mirroring the `innervation_zones` pre-compute pattern guarantees each
    # MU gets a distinct, reproducible draw regardless of worker layout.
    fiber_length_variations_per_MU = [
        get_random_generator().uniform(
            low=-self._var_fiber_length__mm,
            high=self._var_fiber_length__mm,
            size=int(number_of_fibers_per_MUs[mu_idx]),
        )
        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),
    )

    # 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,
    )

    # Helper function to process a single motor unit
    def _process_single_mu(
        MU_index: int,
        electrode_array_original: SurfaceElectrodeArray,
    ) -> tuple[np.ndarray, str]:
        """
        Process a single motor unit (all its fibers) in parallel.

        Parameters
        ----------
        MU_index : int
            Index of the motor unit to process.
        electrode_array_original : SurfaceElectrodeArray
            Original electrode array (will be deep-copied to avoid threading issues).

        Returns
        -------
        tuple[np.ndarray, str]
            Tuple of (array_result, segment_name) where array_result is the accumulated
            MUAP signal for this MU and segment_name is the name for the segment.
        """
        try:
            # Deep copy electrode array only when running in parallel (thread safety)
            # For n_jobs=1 sequential mode, skip the expensive deepcopy
            if n_jobs != 1:
                electrode_array = deepcopy(electrode_array_original)
            else:
                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)

            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)
                return array_result_downsampled, f"MUAP_{MU_index}"

            # Get fiber positions
            position_of_fibers_raw = self._muscle_model.resulting_fiber_assignment(MU_index)
            # Extract magnitude if Quantity, otherwise use as-is
            if hasattr(position_of_fibers_raw, "magnitude"):
                position_of_fibers = position_of_fibers_raw.magnitude
            else:
                position_of_fibers = position_of_fibers_raw

            innervation_zone = innervation_zones[MU_index]

            # Use the parent-precomputed fiber-length variations to keep
            # per-MU RNG state independent of the joblib worker layout
            # (see comment at the `fiber_length_variations_per_MU` site).
            fiber_length_variations = fiber_length_variations_per_MU[MU_index]

            # Pre-compute geometric values for all fibers (optimization: vectorized)
            R_values = np.sqrt(position_of_fibers[:, 0] ** 2 + position_of_fibers[:, 1] ** 2)
            theta_values = np.arctan2(position_of_fibers[:, 1], position_of_fibers[:, 0])
            fiber_lengths = self._mean_fiber_length__mm + fiber_length_variations

            # Matrix optimization variables (local to this MU)
            A_matrix = None
            B_incomplete = None

            # Unified path cache variables
            A_matrix_unified = None
            b_z_cached = None

            # Pre-compute base electrode positions ONCE (avoid per-fiber grid recomputation)
            import quantities as pq
            base_pos_z = electrode_array.pos_z.rescale(pq.mm).magnitude.copy()
            base_pos_theta = electrode_array.pos_theta.rescale(pq.rad).magnitude.copy()
            base_rele = float(electrode_array.electrode_radius__mm.rescale(pq.mm).magnitude)

            # Pre-extract scalar values from self (avoid repeated quantities access in fiber loop)
            Fs_internal = float(self._internal_sampling_frequency__Hz * 1e-3)
            v_conduction = float(self._mean_conduction_velocity__m_s)
            N_internal = int(self._internal_sampling_points)
            M_theta = int(self._sampling_points_in_theta_domain)
            r_total = float(self._radius_total)
            r_bone = float(self._radius_bone__mm)
            th_fat = float(self._fat_thickness__mm)
            th_skin = float(self._skin_thickness__mm)
            sig_rho = float(self._muscle_conductivity_radial__S_m)
            sig_z = float(self._muscle_conductivity_longitudinal__S_m)
            sig_skin_val = float(self._skin_conductivity__S_m)
            sig_fat_val = float(self._fat_conductivity__S_m)

            # Determine IAP kernel length once (same for all fibers)
            if self._iap_kernel_length__mm is not None:
                kernel_length = self._iap_kernel_length__mm
            else:
                IAP_SCALE_FACTOR = 2.5
                kernel_length = self._mean_fiber_length__mm * IAP_SCALE_FACTOR

            # Process each fiber (inner loop - must remain sequential)
            fiber_iter = tqdm(
                range(number_of_fibers),
                desc=f"  MU {MU_index} fibers",
                leave=False,
                disable=not verbose,
            )
            for fiber_number in fiber_iter:
                # Use pre-computed values (optimization: vectorized calculations)
                R = R_values[fiber_number]
                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)

                if self._use_unified:
                    from myogen.simulator.core.emg.fiber_simulation import (
                        simulate_fiber_hybrid,
                    )

                    # Hybrid: time-domain Rosenfalck + frequency-domain volume conductor
                    # Same signature as _simulate_fiber_v2_python (drop-in replacement)
                    phi_temp, A_matrix, B_incomplete = simulate_fiber_hybrid(
                        Fs=Fs_internal,
                        v=v_conduction,
                        N=N_internal,
                        M=M_theta,
                        r=r_total,
                        r_bone=r_bone,
                        th_fat=th_fat,
                        th_skin=th_skin,
                        R=R,
                        L1=L1,
                        L2=L2,
                        zi=innervation_zone,
                        electrode_array=electrode_array,
                        sig_muscle_rho=sig_rho,
                        sig_muscle_z=sig_z,
                        sig_fat=sig_fat_val,
                        sig_skin=sig_skin_val,
                        fiber_length__mm=kernel_length,
                        A_matrix=None if fiber_number == 0 else A_matrix,
                        B_incomplete=None if fiber_number == 0 else B_incomplete,
                        # `simulate_fiber_hybrid` takes a plain bool; map
                        # the tri-state surface flag (None → CPU here, since
                        # hybrid is experimental) before forwarding.
                        use_gpu=bool(use_gpu) if use_gpu is not None else False,
                        theta_offset=-theta,
                        pos_z_precomputed=base_pos_z,
                        pos_theta_precomputed=base_pos_theta,
                        rele_precomputed=base_rele,
                        D1=96.0,
                    )
                else:
                    # Existing frequency-domain path (unchanged)
                    phi_temp, A_matrix, B_incomplete = _simulate_fiber_v2_python(
                        Fs=Fs_internal,
                        v=v_conduction,
                        N=N_internal,
                        M=M_theta,
                        r=r_total,
                        r_bone=r_bone,
                        th_fat=th_fat,
                        th_skin=th_skin,
                        R=R,
                        L1=L1,
                        L2=L2,
                        zi=innervation_zone,
                        electrode_array=electrode_array,
                        sig_muscle_rho=sig_rho,
                        sig_muscle_z=sig_z,
                        sig_fat=sig_fat_val,
                        sig_skin=sig_skin_val,
                        fiber_length__mm=kernel_length,
                        A_matrix=None if fiber_number == 0 else A_matrix,
                        B_incomplete=None if fiber_number == 0 else B_incomplete,
                        # Tri-state forwarded directly: None → auto, True →
                        # require GPU, False → force CPU (resolved inside
                        # `_simulate_fiber_v2_python`).
                        use_gpu=use_gpu,
                        theta_offset=-theta,
                        pos_z_precomputed=base_pos_z,
                        pos_theta_precomputed=base_pos_theta,
                        rele_precomputed=base_rele,
                    )

                array_result_internal += phi_temp

            # Downsample from internal resolution to output resolution
            # resample operates on the last axis (time), which is axis=2
            array_result_downsampled = resample(
                array_result_internal, self._sampling_points_in_t_and_z_domains, axis=2
            )

            return array_result_downsampled, f"MUAP_{MU_index}"

        except Exception as e:
            # Log error and return empty result to avoid crashing entire parallel job
            logging.error(
                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)
            return empty_result, f"MUAP_{MU_index}_FAILED"

    block = Block()
    for array_idx, electrode_array in enumerate(self._electrode_arrays):
        group = Group(name=f"ElectrodeArray_{array_idx}")
        block.groups.append(group)

        # Process only specified motor units in parallel
        n_mus_to_compute = len(self._MUs_to_simulate)
        logging.info(
            f"Processing {n_mus_to_compute}/{n_motor_units} motor units for electrode array {array_idx + 1}/{len(self._electrode_arrays)} using parallel processing..."
        )

        # Process motor units with progress bar
        results = {}
        if n_jobs == 1:
            # Sequential: call directly, no joblib overhead
            for MU_index in tqdm(
                self._MUs_to_simulate,
                desc=f"Electrode Array {array_idx + 1}/{len(self._electrode_arrays)}",
                disable=not verbose,
            ):
                array_result, segment_name = _process_single_mu(MU_index, electrode_array)
                mu_idx = int(segment_name.split("_")[1].split("_")[0])
                results[mu_idx] = (array_result, segment_name)
        else:
            # Parallel: use joblib for multi-core
            with tqdm(
                total=n_mus_to_compute,
                desc=f"Electrode Array {array_idx + 1}/{len(self._electrode_arrays)}",
                disable=not verbose,
            ) as pbar:
                for array_result, segment_name in Parallel(
                    n_jobs=n_jobs,
                    return_as="generator",
                    verbose=0,
                    batch_size="auto",
                )(
                    delayed(_process_single_mu)(MU_index, electrode_array)
                    for MU_index in self._MUs_to_simulate
                ):
                    mu_idx = int(segment_name.split("_")[1].split("_")[0])
                    results[mu_idx] = (array_result, segment_name)
                    pbar.update(1)

        # Calculate actual MUAP duration based on fiber lengths
        # Use iap_kernel_length__mm if specified, otherwise use scaled fiber length from muscle model
        if self._iap_kernel_length__mm is not None:
            kernel_length_mm = self._iap_kernel_length__mm
        else:
            # Scale fiber length to avoid boundary truncation of IAP kernel
            # The IAP kernel needs ~2.5x the fiber length to fully develop and decay
            # This prevents edge artifacts while maintaining proportionality to fiber length
            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

        # Create segments for ALL MUs (maintaining index order)
        # Non-computed MUs get empty signals at output resolution
        for MU_index in range(n_motor_units):
            if MU_index in results:
                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)
                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])
            segment.analogsignals.append(
                create_grid_signal(
                    signal=np.transpose(array_result, (2, 0, 1)) * pq.mV,
                    grid_shape=grid_shape,
                    sampling_rate=effective_sampling_rate__Hz * pq.Hz,
                    t_start=0 * pq.s,
                )
            )

    # Store results privately
    self._muaps__Block = block

    return block

simulate_surface_emg

simulate_surface_emg(spike_train__Block: SPIKE_TRAIN__Block, verbose: bool = True) -> SURFACE_EMG__Block

Generate surface EMG signals for all electrode arrays using the provided spike train block.

This method convolves the pre-computed MUAP templates with the spike trains to synthesize realistic surface EMG signals. The process includes temporal resampling to match the spike train timestep and supports both CPU and GPU acceleration.

Parameters:

Name Type Description Default
spike_train__Block SPIKE_TRAIN__Block

Block containing spike trains organized as segments (pools) with spiketrains.

required
verbose bool

If True, display progress bars. Set to False to disable.

True

Returns:

Type Description
SURFACE_EMG__Block

Surface EMG signals for each electrode array stored in a neo.Block. Results are stored in the surface_emg__tensors property after execution.

Raises:

Type Description
ValueError

If MUAP templates have not been generated. Call simulate_muaps() first.

Source code in myogen/simulator/core/emg/surface/surface_emg.py
def simulate_surface_emg(self, spike_train__Block: SPIKE_TRAIN__Block, verbose: bool = True) -> SURFACE_EMG__Block:
    """
    Generate surface EMG signals for all electrode arrays using the provided spike train block.

    This method convolves the pre-computed MUAP templates with the spike trains
    to synthesize realistic surface EMG signals. The process includes temporal resampling
    to match the spike train timestep and supports both CPU and GPU acceleration.

    Parameters
    ----------
    spike_train__Block : SPIKE_TRAIN__Block
        Block containing spike trains organized as segments (pools) with spiketrains.
    verbose : bool, default=True
        If True, display progress bars. Set to False to disable.

    Returns
    -------
    SURFACE_EMG__Block
        Surface EMG signals for each electrode array stored in a neo.Block.
        Results are stored in the `surface_emg__tensors` property after execution.

    Raises
    ------
    ValueError
        If MUAP templates have not been generated. Call simulate_muaps() first.
    """
    if self._muaps__Block is None:
        raise ValueError("MUAP templates have not been generated. Call simulate_muaps() first.")

    # Store spike train data privately
    self._spike_train__Block = spike_train__Block

    # Extract timestep from the MUAP native sampling rate (not output rate)
    muap_native_rate = float(self._muaps__Block.groups[0].segments[0].analogsignals[0].sampling_rate)
    muap_timestep__ms = float((1 / muap_native_rate) * 1000) * pq.ms

    # Convert spike train block to numpy arrays
    n_pools = len(spike_train__Block.segments)
    n_neurons = len(spike_train__Block.segments[0].spiketrains)

    # Extract spike train durations to determine time length
    first_spiketrain = spike_train__Block.segments[0].spiketrains[0]
    spiketrain_timestep__ms = first_spiketrain.sampling_period.rescale("ms")

    # Bin each pool's spike trains into a boolean occupancy array.
    spike_trains = np.array(
        [
            bin_spike_trains(segment.spiketrains, bin_size=spiketrain_timestep__ms)
            for segment in spike_train__Block.segments
        ]
    )

    # Handle MUs to simulate
    if self._MUs_to_simulate is None:
        MUs_to_simulate = set(range(n_neurons))
    else:
        MUs_to_simulate = set(self._MUs_to_simulate)

    # Create active neuron indices (all neurons are active in each pool for spike train block)
    active_neuron_indices = [list(range(n_neurons)) for _ in range(n_pools)]

    block = Block()

    muap_data_list = [
        np.array([signal_to_grid(seg.analogsignals[0]) for seg in group.segments])
        for group in self._muaps__Block.groups
    ]

    for array_idx, muap_array in enumerate(muap_data_list):
        emg_group = Group(name=f"ElectrodeArray_{array_idx}")
        block.groups.append(emg_group)

        muap_array = np.transpose(muap_array, (0, 2, 3, 1))

        # Temporal resampling: MUAP samples are at native rate, resample to spike train rate
        muap_duration__s = muap_array.shape[3] * muap_timestep__ms.rescale(pq.s).magnitude
        new_muap_time_length = max(
            1,
            np.round(
                muap_duration__s
                / spiketrain_timestep__ms.rescale("s").magnitude
            ).astype(int),
        )

        muap_shapes = np.zeros(
            (
                muap_array.shape[0],
                muap_array.shape[1],
                muap_array.shape[2],
                new_muap_time_length,
            )
        )

        # Time axes for interpolation
        xp_native = np.arange(muap_array.shape[3]) * muap_timestep__ms.rescale(pq.s).magnitude
        x_target = np.arange(new_muap_time_length) * spiketrain_timestep__ms.rescale(pq.s).magnitude

        for muap_nr in range(muap_shapes.shape[0]):
            for row in range(muap_shapes.shape[1]):
                for col in range(muap_shapes.shape[2]):
                    muap_shapes[muap_nr, row, col] = np.interp(
                        x=x_target,
                        xp=xp_native,
                        fp=muap_array[muap_nr, row, col],
                    )

        # n_pools already defined above from spike_train_block
        n_rows = muap_shapes.shape[1]
        n_cols = muap_shapes.shape[2]

        # Initialize result array
        sample_conv = np.convolve(spike_trains[0, 0], muap_shapes[0, 0, 0], mode="same")

        surface_emg = np.zeros((n_pools, n_rows, n_cols, len(sample_conv)))

        # No normalization needed - MUAPs are in absolute units (mV) from biophysical model

        # Perform convolution for each pool using GPU acceleration if available
        if HAS_CUPY:
            # Use GPU acceleration with CuPy
            spike_gpu = cp.asarray(spike_trains)
            muap_gpu = cp.asarray(muap_shapes)
            surface_emg_gpu = cp.zeros((n_pools, n_rows, n_cols, len(sample_conv)))

            for pool_idx in tqdm(
                range(n_pools),
                desc=f"Electrode Array {array_idx + 1}/{len(self._muaps__Block.groups)} Surface EMG (GPU)",
                unit="pools",
                disable=not verbose,
            ):
                pool_active_neurons = set(active_neuron_indices[pool_idx])

                for row_idx in range(n_rows):
                    for col_idx in range(n_cols):
                        # Process all active MUs on GPU
                        convolutions = cp.array(
                            [
                                cp.convolve(
                                    spike_gpu[pool_idx, mu_idx],
                                    muap_gpu[mu_idx, row_idx, col_idx],
                                    mode="same",
                                )
                                for mu_idx in MUs_to_simulate.intersection(pool_active_neurons)
                            ]
                        )
                        # Sum across MUAPs on GPU
                        if len(convolutions) > 0:
                            surface_emg_gpu[pool_idx, row_idx, col_idx] = cp.sum(
                                convolutions, axis=0
                            )

            # Transfer results back to CPU
            surface_emg = cp.asnumpy(surface_emg_gpu)
        else:
            # Fallback to CPU computation with NumPy
            for pool_idx in tqdm(
                range(n_pools),
                desc=f"Electrode Array {array_idx + 1}/{len(self._muaps__Block.groups)} Surface EMG (CPU)",
                unit="pools",
                disable=not verbose,
            ):
                pool_active_neurons = set(active_neuron_indices[pool_idx])

                for row_idx in range(n_rows):
                    for col_idx in range(n_cols):
                        # Process all active MUs
                        convolutions = []
                        for mu_idx in MUs_to_simulate.intersection(pool_active_neurons):
                            conv = np.convolve(
                                spike_trains[pool_idx, mu_idx],
                                muap_shapes[mu_idx, row_idx, col_idx],
                                mode="same",
                            )
                            convolutions.append(conv)

                        if convolutions:
                            surface_emg[pool_idx, row_idx, col_idx] = np.sum(
                                convolutions, axis=0
                            )

        # Temporal resampling
        surface_emg_resampled = np.zeros(
            (
                n_pools,
                n_rows,
                n_cols,
                int(
                    surface_emg.shape[-1]
                    * spiketrain_timestep__ms.rescale(pq.s).magnitude
                    * self._sampling_frequency__Hz
                ),
            )
        )
        # Native (spike-train-rate) and target (output-rate) time axes,
        # built from integer sample counts so each grid matches its data
        # length exactly while preserving the original sampling periods
        # (dt for the native axis, 1 / fs for the target axis). np.arange
        # with a float step can land one element off when N * timestep is
        # not exactly representable in IEEE 754, which makes np.interp
        # raise on an xp/fp length mismatch (issue #12).
        resample_timestep__s = spiketrain_timestep__ms.rescale(pq.s).magnitude
        n_native = surface_emg.shape[-1]
        n_resampled = surface_emg_resampled.shape[-1]
        xp_time = np.arange(n_native) * resample_timestep__s
        x_time = np.arange(n_resampled) / self._sampling_frequency__Hz
        for pool_idx in range(n_pools):
            for row_idx in range(n_rows):
                for col_idx in range(n_cols):
                    surface_emg_resampled[pool_idx, row_idx, col_idx] = np.interp(
                        x=x_time,
                        xp=xp_time,
                        fp=surface_emg[pool_idx, row_idx, col_idx],
                    )

        # Create segments for each motor unit pool within this electrode array group
        for pool_idx in range(n_pools):
            segment = Segment(name=f"Pool_{pool_idx}")
            emg_group.segments.append(segment)

            # Create grid-annotated AnalogSignal for this pool's EMG data
            segment.analogsignals.append(
                create_grid_signal(
                    signal=np.transpose(surface_emg_resampled[pool_idx], (2, 0, 1)) * pq.mV,
                    grid_shape=(n_rows, n_cols),
                    sampling_rate=self._sampling_frequency__Hz * pq.Hz,
                )
            )

    # Store results privately
    self._surface_emg__Block = block
    return block

add_noise

add_noise(snr__dB: float, noise_type: str = 'gaussian') -> SURFACE_EMG__Block

Add noise to all electrode arrays.

This method adds realistic noise to the simulated surface EMG signals based on a specified signal-to-noise ratio. The noise is calculated and applied independently for each electrode channel to ensure that channels with different signal amplitudes maintain the specified SNR.

Parameters:

Name Type Description Default
snr__dB float

Signal-to-noise ratio in dB. Higher values result in cleaner signals. Typical physiological EMG has SNR ranging from 10-40 dB. The SNR is applied independently to each electrode channel.

required
noise_type str

Type of noise to add. Currently supports "gaussian" for white noise.

"gaussian"

Returns:

Type Description
SURFACE_EMG__Block

Noisy EMG signals for each electrode array as a neo.Block. Results are stored in the noisy_surface_emg__Block property after execution.

Raises:

Type Description
ValueError

If surface EMG has not been simulated. Call simulate_surface_emg() first.

Notes

The noise is computed per-channel (per electrode) to maintain the specified SNR independently across all channels. This ensures that electrodes with different signal amplitudes receive appropriately scaled noise.

Source code in myogen/simulator/core/emg/surface/surface_emg.py
def add_noise(self, snr__dB: float, noise_type: str = "gaussian") -> SURFACE_EMG__Block:
    """
    Add noise to all electrode arrays.

    This method adds realistic noise to the simulated surface EMG signals
    based on a specified signal-to-noise ratio. The noise is calculated
    and applied independently for each electrode channel to ensure that
    channels with different signal amplitudes maintain the specified SNR.

    Parameters
    ----------
    snr__dB : float
        Signal-to-noise ratio in dB. Higher values result in cleaner signals.
        Typical physiological EMG has SNR ranging from 10-40 dB.
        The SNR is applied independently to each electrode channel.
    noise_type : str, default="gaussian"
        Type of noise to add. Currently supports "gaussian" for white noise.

    Returns
    -------
    SURFACE_EMG__Block
        Noisy EMG signals for each electrode array as a neo.Block.
        Results are stored in the `noisy_surface_emg__Block` property after execution.

    Raises
    ------
    ValueError
        If surface EMG has not been simulated. Call simulate_surface_emg() first.

    Notes
    -----
    The noise is computed per-channel (per electrode) to maintain the specified
    SNR independently across all channels. This ensures that electrodes with
    different signal amplitudes receive appropriately scaled noise.
    """
    if self._surface_emg__Block is None:
        raise ValueError(
            "Surface EMG has not been simulated. Call simulate_surface_emg() first."
        )

    noisy_block = Block()

    for array_idx, emg_group in enumerate(self._surface_emg__Block.groups):
        noisy_group = Group(name=f"ElectrodeArray_{array_idx}")
        noisy_block.groups.append(noisy_group)

        for pool_idx, segment in enumerate(emg_group.segments):
            noisy_segment = Segment(name=f"Pool_{pool_idx}")
            noisy_group.segments.append(noisy_segment)

            # Get the EMG signal data
            emg_signal = segment.analogsignals[0]
            grid_shape = emg_signal.annotations["grid_shape"]
            emg_array = signal_to_grid(emg_signal)  # Shape: (time, rows, cols)

            # Calculate signal power PER CHANNEL (per electrode)
            # Mean along time axis (axis=0) gives power per spatial location
            signal_power_per_channel = np.mean(emg_array**2, axis=0)  # Shape: (rows, cols)

            # Calculate noise power per channel
            snr_linear = 10 ** (snr__dB / 10)
            noise_power_per_channel = signal_power_per_channel / snr_linear
            noise_std_per_channel = np.sqrt(noise_power_per_channel)  # Shape: (rows, cols)

            # Generate noise
            if noise_type.lower() == "gaussian":
                # Generate standard normal noise, then scale per channel
                noise = get_random_generator().normal(loc=0.0, scale=1.0, size=emg_array.shape)
                # Broadcast noise_std_per_channel along time axis
                # noise shape: (time, rows, cols)
                # noise_std_per_channel shape: (rows, cols)
                # Broadcasting: (time, rows, cols) * (1, rows, cols)
                noise = noise * noise_std_per_channel[np.newaxis, :, :]
            else:
                raise ValueError(f"Unsupported noise type: {noise_type}")

            # Add noise
            noisy_emg = emg_array + noise

            # Create new grid-annotated AnalogSignal with noise
            noisy_segment.analogsignals.append(
                create_grid_signal(
                    signal=noisy_emg * emg_signal.units,
                    grid_shape=grid_shape,
                    t_start=emg_signal.t_start,
                    sampling_rate=emg_signal.sampling_rate,
                )
            )

    # Store results privately
    self._noisy_surface_emg__Block = noisy_block
    return noisy_block

IntramuscularEMG

IntramuscularEMG(muscle_model: Muscle, electrode_array: IntramuscularElectrodeArray, sampling_frequency__Hz: Quantity__Hz = 10240.0 * Hz, spatial_resolution__mm: Quantity__mm = 0.01 * mm, endplate_center__percent: float = 50, nmj_jitter__s: Quantity__s = 3.5e-05 * s, branch_cvs__m_per_s: tuple[Quantity__m_per_s, Quantity__m_per_s] = (5.0 * m / s, 2.0 * m / s), MUs_to_simulate: list[int] | None = None)

Intramuscular Electromyography (iEMG) Simulation.

This class provides a comprehensive simulation framework for generating intramuscular EMG signals detected by needle electrodes.

Parameters:

Name Type Description Default
muscle_model Muscle

Pre-computed muscle model (see myogen.simulator.Muscle).

required
electrode_array IntramuscularElectrodeArray

Intramuscular electrode array configuration to use for simulation (see myogen.simulator.IntramuscularElectrodeArray).

required
sampling_frequency__Hz Quantity__Hz

Sampling frequency in Hz for EMG simulation. Default is set to 10240 Hz as used by the Quattrocento (OT Bioelettronica, Turin, Italy) system.

10240.0 * pq.Hz
spatial_resolution__mm Quantity__mm

Spatial resolution for fiber action potential calculation in mm. Default is set to 0.01 mm.

0.01 * pq.mm
endplate_center__percent float

Percentage of muscle length where the endplate is located. By default, the endplate is located at the center of the muscle (50% of the muscle length).

50
nmj_jitter__s Quantity__s

Standard deviation of neuromuscular junction jitter in seconds. Default is set to 35e-6 s as determined by Konstantin et al. 2020 [1].

35e-6 * pq.s
branch_cvs__m_per_s tuple[Quantity__m_per_s, Quantity__m_per_s]

Conduction velocities for the two-layer model of the neuromuscular junction in m/s. Default is set to (5.0, 2.0) m/s as determined by Konstantin et al. 2020 [1].

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].

(5.0 * pq.m / pq.s, 2.0 * pq.m / pq.s)
MUs_to_simulate list[int]

Indices of motor units to simulate. If None, all motor units are simulated. Default is None. For computational efficiency, consider simulating subsets for initial analysis. Indices correspond to the recruitment order (0 is recruited first).

None

Attributes:

Name Type Description
muaps__Block INTRAMUSCULAR_MUAP__Block

Intramuscular MUAP shapes for the electrode array as a neo.Block. Available after simulate_muaps().

intramuscular_emg__Block INTRAMUSCULAR_EMG__Block

Intramuscular EMG signals for the electrode array as a neo.Block. Available after simulate_intramuscular_emg().

noisy_intramuscular_emg__Block INTRAMUSCULAR_EMG__Block

Noisy intramuscular EMG signals for the electrode array as a neo.Block. Available after add_noise().

spike_train__Block SPIKE_TRAIN__Block

Spike train block used for EMG generation. Available after simulate_intramuscular_emg().

References

[1] Konstantin, A., Yu, T., Le Carpentier, E., Aoustin, Y., Farina, D., 2020. Simulation of Motor Unit Action Potential Recordings From Intramuscular Multichannel Scanning Electrodes. IEEE Transactions on Biomedical Engineering 67, 2005–2014. https://doi.org/10.1109/TBME.2019.2953680

Source code in myogen/simulator/core/emg/intramuscular/intramuscular_emg.py
def __init__(
    self,
    muscle_model: Muscle,
    electrode_array: IntramuscularElectrodeArray,
    sampling_frequency__Hz: Quantity__Hz = 10240.0 * pq.Hz,
    spatial_resolution__mm: Quantity__mm = 0.01 * pq.mm,
    endplate_center__percent: float = 50,
    nmj_jitter__s: Quantity__s = 35e-6 * pq.s,
    branch_cvs__m_per_s: tuple[Quantity__m_per_s, Quantity__m_per_s] = (
        5.0 * pq.m / pq.s,
        2.0 * pq.m / pq.s,
    ),
    MUs_to_simulate: list[int] | None = None,
):
    # Immutable public arguments - never modify these
    self.muscle_model = muscle_model
    self.electrode_array = electrode_array
    self.sampling_frequency__Hz = sampling_frequency__Hz
    self.spatial_resolution__mm = spatial_resolution__mm
    self.endplate_center__percent = endplate_center__percent
    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

    # Private copies for internal modifications (extract magnitudes)
    self._muscle_model = muscle_model
    self._electrode_array = electrode_array
    self._sampling_frequency__Hz = float(sampling_frequency__Hz.rescale(pq.Hz).magnitude)
    self._spatial_resolution__mm = float(spatial_resolution__mm.rescale(pq.mm).magnitude)
    self._endplate_center__percent = endplate_center__percent
    self._nmj_jitter__s = float(nmj_jitter__s.rescale(pq.s).magnitude)
    self._branch_cvs__m_per_s = (
        float(branch_cvs__m_per_s[0].rescale(pq.m / pq.s).magnitude),
        float(branch_cvs__m_per_s[1].rescale(pq.m / pq.s).magnitude),
    )
    self._MUs_to_simulate = MUs_to_simulate

    # Derived parameters - immutable public access
    self.endplate_center__mm = self._muscle_model.length__mm * (
        self._endplate_center__percent / 100.0
    )

    # Private copies for internal modifications
    self._branch_cvs__mm_per_s: tuple[float, float] = (
        self._branch_cvs__m_per_s[0] * 1000.0,
        self._branch_cvs__m_per_s[1] * 1000.0,
    )
    # Extract magnitude for internal use (calculations expect floats)
    length_mm = (
        float(self._muscle_model.length__mm.rescale(pq.mm).magnitude)
        if hasattr(self._muscle_model.length__mm, "magnitude")
        else float(self._muscle_model.length__mm)
    )
    self._endplate_center__mm = length_mm * (self._endplate_center__percent / 100.0)

    # Derived parameters - private for internal use
    self._dt = 1.0 / self._sampling_frequency__Hz
    self._dz = self._spatial_resolution__mm
    self._n_motor_units = len(self._muscle_model.recruitment_thresholds)

    # Motor unit selection
    if self._MUs_to_simulate is None:
        self._MUs_to_simulate = list(range(self._n_motor_units))
    else:
        self._MUs_to_simulate = self._MUs_to_simulate

    # Motor unit simulations - private storage
    self._motor_units: list[MotorUnitSim] = []  # List of motor unit simulators
    self._muaps__Block: Optional[INTRAMUSCULAR_MUAP__Block] = None
    self._max_muap_length: int = 0

    # Simulation results - stored privately, accessed via properties
    self._intramuscular_emg__Block: Optional[INTRAMUSCULAR_EMG__Block] = None
    self._noisy_intramuscular_emg__Block: Optional[INTRAMUSCULAR_EMG__Block] = None
    self._spike_train__Block: Optional[SPIKE_TRAIN__Block] = None

muaps__Block property

Intramuscular MUAP shapes for the electrode array.

Returns:

Type Description
INTRAMUSCULAR_MUAP__Block

Intramuscular MUAP templates for the electrode array as a neo.Block.

Raises:

Type Description
ValueError

If MUAP templates have not been computed yet.

intramuscular_emg__Block property

intramuscular_emg__Block: INTRAMUSCULAR_EMG__Block

Intramuscular EMG signals for the electrode array.

Returns:

Type Description
INTRAMUSCULAR_EMG__Block

Intramuscular EMG signals for the electrode array as a neo.Block.

Raises:

Type Description
ValueError

If intramuscular EMG has not been computed yet.

noisy_intramuscular_emg__Block property

noisy_intramuscular_emg__Block: INTRAMUSCULAR_EMG__Block

Noisy intramuscular EMG signals for the electrode array.

Returns:

Type Description
INTRAMUSCULAR_EMG__Block

Noisy intramuscular EMG signals for the electrode array as a neo.Block.

Raises:

Type Description
ValueError

If noisy intramuscular EMG has not been computed yet.

spike_train__Block property

spike_train__Block: SPIKE_TRAIN__Block

Spike train block used for EMG generation.

Returns:

Type Description
SPIKE_TRAIN__Block

The spike train block used in the simulation.

Raises:

Type Description
ValueError

If spike train block has not been set yet.

simulate_muaps

simulate_muaps(n_jobs: int = -2, verbose: bool = True) -> INTRAMUSCULAR_MUAP__Block

Simulate MUAPs for all electrode arrays using the provided muscle model.

This method generates intramuscular Motor Unit Action Potential (MUAP) templates by simulating individual motor units with realistic neuromuscular junction distributions and fiber action potential propagation.

Parameters:

Name Type Description Default
n_jobs int

Number of parallel workers for motor unit processing. Default is -2. - n_jobs=-1: Use all CPU cores - n_jobs=-2: Use all cores except one (recommended, keeps system responsive) - n_jobs=-3: Use all cores except two - n_jobs=1: No parallelization - n_jobs=N: Use exactly N cores

-2
verbose bool

If True, display progress bars. Set to False to disable.

True

Returns:

Type Description
INTRAMUSCULAR_MUAP__Block

Intramuscular MUAP shapes for all motor units stored in a neo.Block. Results are stored in the muaps__Block property after execution.

Notes

This method must be called before simulate_intramuscular_emg(). The process includes: (1) motor unit initialization, (2) neuromuscular junction simulation, and (3) MUAP calculation with spatial filtering.

Source code in myogen/simulator/core/emg/intramuscular/intramuscular_emg.py
def simulate_muaps(self, n_jobs: int = -2, verbose: bool = True) -> INTRAMUSCULAR_MUAP__Block:
    """
    Simulate MUAPs for all electrode arrays using the provided muscle model.

    This method generates intramuscular Motor Unit Action Potential (MUAP) templates
    by simulating individual motor units with realistic neuromuscular junction
    distributions and fiber action potential propagation.

    Parameters
    ----------
    n_jobs : int, optional
        Number of parallel workers for motor unit processing. Default is -2.
        - n_jobs=-1: Use all CPU cores
        - n_jobs=-2: Use all cores except one (recommended, keeps system responsive)
        - n_jobs=-3: Use all cores except two
        - n_jobs=1: No parallelization
        - n_jobs=N: Use exactly N cores
    verbose : bool, default=True
        If True, display progress bars. Set to False to disable.

    Returns
    -------
    INTRAMUSCULAR_MUAP__Block
        Intramuscular MUAP shapes for all motor units stored in a neo.Block.
        Results are stored in the `muaps__Block` property after execution.

    Notes
    -----
    This method must be called before simulate_intramuscular_emg(). The process
    includes: (1) motor unit initialization, (2) neuromuscular junction simulation,
    and (3) MUAP calculation with spatial filtering.
    """
    self._initialize_motor_units(verbose=verbose)
    self._simulate_neuromuscular_junctions(verbose=verbose)
    return self._calculate_muaps(n_jobs=n_jobs, verbose=verbose)

simulate_intramuscular_emg

simulate_intramuscular_emg(spike_train__Block: SPIKE_TRAIN__Block, verbose: bool = True) -> INTRAMUSCULAR_EMG__Block

Generate intramuscular EMG signals using the provided spike train block.

This method convolves the pre-computed MUAP templates with spike trains to synthesize realistic intramuscular EMG signals. The process includes temporal resampling and supports both CPU and GPU acceleration for efficient computation.

Parameters:

Name Type Description Default
spike_train__Block SPIKE_TRAIN__Block

Block containing spike trains organized as segments (pools) with spiketrains.

required
verbose bool

If True, display progress bars. Set to False to disable.

True

Returns:

Type Description
INTRAMUSCULAR_EMG__Block

Intramuscular EMG signals for the electrode array stored in a neo.Block. Results are stored in the intramuscular_emg__Block property after execution.

Raises:

Type Description
ValueError

If MUAP templates have not been generated. Call simulate_muaps() first.

Source code in myogen/simulator/core/emg/intramuscular/intramuscular_emg.py
def simulate_intramuscular_emg(
    self,
    spike_train__Block: SPIKE_TRAIN__Block,
    verbose: bool = True,
) -> INTRAMUSCULAR_EMG__Block:
    """
    Generate intramuscular EMG signals using the provided spike train block.

    This method convolves the pre-computed MUAP templates with spike trains
    to synthesize realistic intramuscular EMG signals. The process includes temporal
    resampling and supports both CPU and GPU acceleration for efficient computation.

    Parameters
    ----------
    spike_train__Block : SPIKE_TRAIN__Block
        Block containing spike trains organized as segments (pools) with spiketrains.
    verbose : bool, default=True
        If True, display progress bars. Set to False to disable.

    Returns
    -------
    INTRAMUSCULAR_EMG__Block
        Intramuscular EMG signals for the electrode array stored in a neo.Block.
        Results are stored in the `intramuscular_emg__Block` property after execution.

    Raises
    ------
    ValueError
        If MUAP templates have not been generated. Call simulate_muaps() first.
    """
    if self._muaps__Block is None:
        raise ValueError("MUAP templates have not been generated. Call simulate_muaps() first.")

    # Store spike train data privately
    self._spike_train__Block = spike_train__Block

    # Handle MUs to simulate
    if self._MUs_to_simulate is None:
        MUs_to_simulate = set(
            range(len(self._muscle_model.resulting_number_of_innervated_fibers))
        )
    else:
        MUs_to_simulate = set(self._MUs_to_simulate)

    # Extract MUAP data from Block and pad to same length
    muap_data_list = [seg.analogsignals[0].magnitude for seg in self._muaps__Block.segments]

    # Find the maximum length among all MUAPs
    max_length = max(muap.shape[0] for muap in muap_data_list)
    n_electrodes = muap_data_list[0].shape[1]

    # Pad all MUAPs to the same length (centered, pad both sides)
    padded_muaps = []
    for muap in muap_data_list:
        pad_total = max_length - muap.shape[0]
        if pad_total > 0:
            pad_left = pad_total // 2
            pad_right = pad_total - pad_left
            pad_width = ((pad_left, pad_right), (0, 0))
            padded_muap = np.pad(muap, pad_width, mode="constant", constant_values=0)
        else:
            padded_muap = muap
        padded_muaps.append(padded_muap)

    muap_array = np.array(padded_muaps)

    # Extract timestep from the first spike train
    first_spiketrain = spike_train__Block.segments[0].spiketrains[0]
    spiketrain_timestep__ms = first_spiketrain.sampling_period.rescale("ms")

    target_length = int(
        np.round(
            muap_array.shape[1]
            / self._sampling_frequency__Hz
            * 1
            / (spiketrain_timestep__ms.rescale("s").magnitude)
        )
    )
    muap_shapes = np.zeros((muap_array.shape[0], muap_array.shape[2], target_length))
    for muap_nr in range(muap_shapes.shape[0]):
        for electrode_nr in range(muap_shapes.shape[1]):
            muap_shapes[muap_nr, electrode_nr] = np.interp(
                np.linspace(
                    0,
                    muap_array.shape[1] / self._sampling_frequency__Hz,
                    target_length,
                    endpoint=False,
                ),
                np.arange(
                    0,
                    muap_array.shape[1] / self._sampling_frequency__Hz,
                    1 / self._sampling_frequency__Hz,
                ),
                muap_array[muap_nr, :, electrode_nr],
            )

    # Convert spike train block to numpy arrays
    n_pools = len(spike_train__Block.segments)
    n_neurons = len(spike_train__Block.segments[0].spiketrains)
    n_electrodes = muap_shapes.shape[1]

    # Bin each pool's spike trains into a boolean occupancy array.
    spike_trains = np.array(
        [
            bin_spike_trains(segment.spiketrains, bin_size=spiketrain_timestep__ms)
            for segment in spike_train__Block.segments
        ]
    )

    # Create active neuron indices (all neurons are active in each pool for spike train block)
    active_neuron_indices = [list(range(n_neurons)) for _ in range(n_pools)]

    # Initialize result array
    sample_conv = np.convolve(
        spike_trains[0, 0],
        muap_shapes[0, 0],
        mode="same",
    )
    intramuscular_emg = np.zeros((n_pools, n_electrodes, len(sample_conv)))

    # Normalize MUAP shapes before convolution
    # Note: Unlike surface EMG, intramuscular MUAP amplitudes from biophysical
    # calculations are in arbitrary units and must be normalized to prevent
    # numerical overflow during convolution. Final EMG amplitudes are determined
    # by the spike train convolution, not the raw MUAP amplitudes.
    max_muap_amplitude = np.max(np.abs(muap_shapes))
    if max_muap_amplitude > 0:
        muap_shapes /= max_muap_amplitude

    # Perform convolution for each pool using GPU acceleration if available
    if HAS_CUPY:
        # Use GPU acceleration with CuPy
        spike_gpu = cp.asarray(spike_trains)
        muap_gpu = cp.asarray(muap_shapes)
        intramuscular_emg_gpu = cp.zeros((n_pools, n_electrodes, len(sample_conv)))

        for pool_idx in tqdm(
            range(n_pools),
            desc="Intramuscular EMG (GPU)",
            unit="pool",
            disable=not verbose,
        ):
            pool_active_neurons = set(active_neuron_indices[pool_idx])

            for e_idx in range(n_electrodes):
                # Process all active MUs on GPU
                convolutions = []
                for mu_idx in MUs_to_simulate.intersection(pool_active_neurons):
                    # Use mu_idx directly since muap_shapes now contains all MUs
                    if mu_idx < muap_gpu.shape[0]:
                        conv = cp.convolve(
                            spike_gpu[pool_idx, mu_idx],
                            muap_gpu[mu_idx, e_idx],
                            mode="same",
                        )
                        convolutions.append(conv)

                convolutions = cp.array(convolutions) if convolutions else cp.array([])
                # Sum across MUAPs on GPU
                if len(convolutions) > 0:
                    intramuscular_emg_gpu[pool_idx, e_idx] = cp.sum(convolutions, axis=0)

        # Transfer results back to CPU
        intramuscular_emg = cp.asnumpy(intramuscular_emg_gpu)
    else:
        # Fallback to CPU computation with NumPy
        for pool_idx in tqdm(
            range(n_pools),
            desc="Intramuscular EMG (CPU)",
            unit="pool",
            disable=not verbose,
        ):
            pool_active_neurons = set(active_neuron_indices[pool_idx])

            for e_idx in range(n_electrodes):
                # Process all active MUs
                convolutions = []
                for mu_idx in MUs_to_simulate.intersection(pool_active_neurons):
                    # Use mu_idx directly since muap_shapes now contains all MUs
                    if mu_idx < muap_shapes.shape[0]:
                        conv = np.convolve(
                            spike_trains[pool_idx, mu_idx],
                            muap_shapes[mu_idx, e_idx],
                            mode="same",
                        )
                        convolutions.append(conv)

                if convolutions:
                    intramuscular_emg[pool_idx, e_idx] = np.sum(convolutions, axis=0)

    # Temporal resampling
    intramuscular_emg_resampled = np.zeros(
        (
            n_pools,
            n_electrodes,
            int(
                intramuscular_emg.shape[-1]
                * spiketrain_timestep__ms.rescale("s").magnitude
                * self._sampling_frequency__Hz
            ),
        )
    )
    for pool_idx in range(n_pools):
        for e_idx in range(n_electrodes):
            intramuscular_emg_resampled[pool_idx, e_idx] = np.interp(
                x=np.arange(
                    start=0,
                    stop=intramuscular_emg.shape[-1]
                    * spiketrain_timestep__ms.rescale("s").magnitude,
                    step=1 / self._sampling_frequency__Hz,
                ),
                xp=np.arange(
                    start=0,
                    stop=intramuscular_emg.shape[-1]
                    * spiketrain_timestep__ms.rescale("s").magnitude,
                    step=spiketrain_timestep__ms.rescale("s").magnitude,
                ),
                fp=intramuscular_emg[pool_idx, e_idx],
            )

    # Create neo Block structure
    block = Block()

    # Create segments for each motor unit pool
    for pool_idx in range(n_pools):
        segment = Segment(name=f"Pool_{pool_idx}")
        block.segments.append(segment)

        # Create AnalogSignal for this pool's EMG data
        segment.analogsignals.append(
            AnalogSignal(
                intramuscular_emg_resampled[pool_idx].T * pq.dimensionless,
                t_start=0 * pq.ms,
                sampling_rate=self._sampling_frequency__Hz * pq.Hz,
            )
        )

    # Store results privately
    self._intramuscular_emg__Block = block
    return block

add_noise

add_noise(snr__dB: float, noise_type: str = 'gaussian', *, spectral_slope: float = -0.5, excess_kurtosis: float = 3.0, powerline_hz: float = 50.0, powerline_amplitude: float = 0.1, powerline_harmonic_ratios: Optional[list[float]] = None, powerline_frequency_drift_hz: float = 0.3, powerline_amplitude_modulation_depth: float = 0.15, peak_hz: float = 1000.0, analog_hpf_hz: float = 10.0, baseline_drift_rms_uv: float = 0.0, baseline_drift_alpha: float = 1.75, baseline_drift_low_hz: float | None = None, baseline_drift_high_hz: float = 1.0) -> INTRAMUSCULAR_EMG__Block

Add noise to the electrode array.

Two noise models are available:

  • "gaussian" – white Gaussian noise (legacy default). Each electrode channel gets independent normal noise scaled to hit the requested per-channel SNR.
  • "realistic" – spectrally colored noise calibrated against real bipolar fine-wire iEMG recordings (TEP / Synergy studies). 1/f-like base, mid-band spectral emphasis from electrode–amplifier bandwidth, heavy tails from cross-talk artifacts, and additive 50/60 Hz powerline interference with harmonics. See myogen.utils.emg_noise for the math.

Per-channel SNR is preserved across both modes: each electrode's noise RMS is computed from that channel's own signal RMS so electrodes with different amplitudes get appropriately scaled noise.

Parameters:

Name Type Description Default
snr__dB float

Signal-to-noise ratio in dB. Higher values result in cleaner signals. Typical intramuscular EMG has SNR ranging from 15–50 dB. Applied independently to each electrode channel.

required
noise_type (gaussian, realistic)

Noise model to use. "realistic" activates the colored noise pipeline; the remaining keyword-only parameters then apply.

"gaussian"
spectral_slope float

PSD slope in log–log space for the colored-noise base. 0 = white, -1 = pink. Real iEMG: -0.4 to -0.8. Ignored when noise_type="gaussian".

-0.5
excess_kurtosis float

Target excess kurtosis (0 = Gaussian). Controlled isometric iEMG: 1–6. Dynamic tasks: higher. Ignored when noise_type="gaussian".

3.0
powerline_hz float

Powerline interference frequency. Use 60.0 for North America, 0.0 to disable. Ignored when noise_type="gaussian".

50.0
powerline_amplitude float

Powerline fundamental amplitude as a fraction of noise RMS. Set to 0 (or powerline_hz=0) to disable. Ignored when noise_type="gaussian".

0.1
powerline_harmonic_ratios list of float

Per-harmonic amplitude ratios relative to the fundamental. None uses the default [1.0, 0.5, 0.3, 0.15, 0.08] (fundamental + 4 harmonics). Ignored when noise_type="gaussian".

None
powerline_frequency_drift_hz float

Standard deviation (Hz) of the slow random walk of the mains instantaneous frequency. Broadens each line peak from a delta into a ~2-5 Hz FWHM bump (typical of real recordings). Set to 0 for a pure-tone powerline. Ignored when noise_type="gaussian".

0.3
powerline_amplitude_modulation_depth float

Fractional AM depth (±15% by default) applied to each powerline harmonic — adds narrow sidebands within ±2 Hz of every line. Set to 0 to disable. Ignored when noise_type="gaussian".

0.15
peak_hz float

Center frequency of the mid-band spectral emphasis from electrode–amplifier bandwidth interaction. Ignored when noise_type="gaussian".

1000.0
baseline_drift_rms_uv float

Target RMS (post-HPF) of a band-limited 1/f^α baseline drift representing electrode/interface noise (Huigen 2002, Gondran 1996). Same units as the EMG block. Set to 0 (the default) to disable, preserving legacy behaviour. The spectral form is paper-constrained; the amplitude defaults are not validated for intramuscular EMG — calibrate against real recordings via myogen.utils.calibrate_baseline_drift_profile. Broadband movement artifacts (0–20 Hz, De Luca 2010) are a separate phenomenon out of scope here. Ignored when noise_type="gaussian".

0.0
baseline_drift_alpha float

Drift PSD slope α (PSD ∝ 1/f^α). Midpoint of the [1.5, 2.0] electrode-noise regime. Must be > 0. Ignored when noise_type="gaussian".

1.75
baseline_drift_low_hz float or None

Lower edge of the drift band, in Hz. None resolves to the lowest nonzero FFT bin available for the simulation length. Ignored when noise_type="gaussian".

None
baseline_drift_high_hz float

Upper edge of the drift band, in Hz. Default keeps the knob scoped to sub-1 Hz baseline wander. Ignored when noise_type="gaussian".

1.0

Returns:

Type Description
INTRAMUSCULAR_EMG__Block

Noisy intramuscular EMG signals for the electrode array as a neo.Block. Results are also stored on noisy_intramuscular_emg__Block.

Raises:

Type Description
ValueError

If intramuscular EMG has not been simulated (call simulate_intramuscular_emg first) or noise_type is unrecognized.

Source code in myogen/simulator/core/emg/intramuscular/intramuscular_emg.py
def add_noise(
    self,
    snr__dB: float,
    noise_type: str = "gaussian",
    *,
    spectral_slope: float = -0.5,
    excess_kurtosis: float = 3.0,
    powerline_hz: float = 50.0,
    powerline_amplitude: float = 0.1,
    powerline_harmonic_ratios: Optional[list[float]] = None,
    powerline_frequency_drift_hz: float = 0.3,
    powerline_amplitude_modulation_depth: float = 0.15,
    peak_hz: float = 1000.0,
    analog_hpf_hz: float = 10.0,
    baseline_drift_rms_uv: float = 0.0,
    baseline_drift_alpha: float = 1.75,
    baseline_drift_low_hz: float | None = None,
    baseline_drift_high_hz: float = 1.0,
) -> INTRAMUSCULAR_EMG__Block:
    """
    Add noise to the electrode array.

    Two noise models are available:

    * ``"gaussian"`` – white Gaussian noise (legacy default). Each
      electrode channel gets independent normal noise scaled to hit
      the requested per-channel SNR.
    * ``"realistic"`` – spectrally colored noise calibrated against
      real bipolar fine-wire iEMG recordings (TEP / Synergy studies).
      1/f-like base, mid-band spectral emphasis from
      electrode–amplifier bandwidth, heavy tails from cross-talk
      artifacts, and additive 50/60 Hz powerline interference with
      harmonics. See `myogen.utils.emg_noise` for the math.

    Per-channel SNR is preserved across both modes: each electrode's
    noise RMS is computed from that channel's own signal RMS so
    electrodes with different amplitudes get appropriately scaled
    noise.

    Parameters
    ----------
    snr__dB : float
        Signal-to-noise ratio in dB. Higher values result in cleaner
        signals. Typical intramuscular EMG has SNR ranging from
        15–50 dB. Applied independently to each electrode channel.
    noise_type : {"gaussian", "realistic"}, default="gaussian"
        Noise model to use. ``"realistic"`` activates the colored
        noise pipeline; the remaining keyword-only parameters then
        apply.
    spectral_slope : float, default=-0.5
        PSD slope in log–log space for the colored-noise base.
        ``0`` = white, ``-1`` = pink. Real iEMG: -0.4 to -0.8.
        Ignored when ``noise_type="gaussian"``.
    excess_kurtosis : float, default=3.0
        Target excess kurtosis (``0`` = Gaussian). Controlled
        isometric iEMG: 1–6. Dynamic tasks: higher.
        Ignored when ``noise_type="gaussian"``.
    powerline_hz : float, default=50.0
        Powerline interference frequency. Use ``60.0`` for North
        America, ``0.0`` to disable.
        Ignored when ``noise_type="gaussian"``.
    powerline_amplitude : float, default=0.1
        Powerline fundamental amplitude as a fraction of noise RMS.
        Set to ``0`` (or ``powerline_hz=0``) to disable.
        Ignored when ``noise_type="gaussian"``.
    powerline_harmonic_ratios : list of float, optional
        Per-harmonic amplitude ratios relative to the fundamental.
        ``None`` uses the default ``[1.0, 0.5, 0.3, 0.15, 0.08]``
        (fundamental + 4 harmonics).
        Ignored when ``noise_type="gaussian"``.
    powerline_frequency_drift_hz : float, default=0.3
        Standard deviation (Hz) of the slow random walk of the
        mains instantaneous frequency. Broadens each line peak
        from a delta into a ~2-5 Hz FWHM bump (typical of real
        recordings). Set to 0 for a pure-tone powerline.
        Ignored when ``noise_type="gaussian"``.
    powerline_amplitude_modulation_depth : float, default=0.15
        Fractional AM depth (±15% by default) applied to each
        powerline harmonic — adds narrow sidebands within ±2 Hz
        of every line. Set to 0 to disable.
        Ignored when ``noise_type="gaussian"``.
    peak_hz : float, default=1000.0
        Center frequency of the mid-band spectral emphasis from
        electrode–amplifier bandwidth interaction.
        Ignored when ``noise_type="gaussian"``.
    baseline_drift_rms_uv : float, default=0.0
        Target RMS (post-HPF) of a band-limited 1/f^α baseline
        drift representing electrode/interface noise (Huigen 2002,
        Gondran 1996). Same units as the EMG block. Set to 0 (the
        default) to disable, preserving legacy behaviour.
        The spectral form is paper-constrained; the amplitude
        defaults are *not* validated for intramuscular EMG —
        calibrate against real recordings via
        `myogen.utils.calibrate_baseline_drift_profile`.
        Broadband movement artifacts (0–20 Hz, De Luca 2010) are
        a separate phenomenon out of scope here.
        Ignored when ``noise_type="gaussian"``.
    baseline_drift_alpha : float, default=1.75
        Drift PSD slope α (PSD ∝ 1/f^α). Midpoint of the [1.5, 2.0]
        electrode-noise regime. Must be > 0.
        Ignored when ``noise_type="gaussian"``.
    baseline_drift_low_hz : float or None, default=None
        Lower edge of the drift band, in Hz. ``None`` resolves
        to the lowest nonzero FFT bin available for the
        simulation length.
        Ignored when ``noise_type="gaussian"``.
    baseline_drift_high_hz : float, default=1.0
        Upper edge of the drift band, in Hz. Default keeps the
        knob scoped to sub-1 Hz baseline wander.
        Ignored when ``noise_type="gaussian"``.

    Returns
    -------
    INTRAMUSCULAR_EMG__Block
        Noisy intramuscular EMG signals for the electrode array as a
        ``neo.Block``. Results are also stored on
        ``noisy_intramuscular_emg__Block``.

    Raises
    ------
    ValueError
        If intramuscular EMG has not been simulated (call
        `simulate_intramuscular_emg` first) or ``noise_type``
        is unrecognized.
    """
    if self._intramuscular_emg__Block is None:
        raise ValueError(
            "Intramuscular EMG has not been simulated. Call simulate_intramuscular_emg() first."
        )

    noise_kind = noise_type.lower()
    if noise_kind not in ("gaussian", "realistic"):
        raise ValueError(f"Unsupported noise type: {noise_type}")

    noisy_block = Block()
    rng_global = get_random_generator()

    for pool_idx, segment in enumerate(self._intramuscular_emg__Block.segments):
        noisy_segment = Segment(name=f"Pool_{pool_idx}")
        noisy_block.segments.append(noisy_segment)

        # Get the EMG signal data
        emg_signal = segment.analogsignals[0]
        emg_array = emg_signal.magnitude  # Shape: (time, n_electrodes)

        # Calculate signal power PER CHANNEL (per electrode)
        # Mean along time axis (axis=0) gives power per electrode
        signal_power_per_channel = np.mean(emg_array**2, axis=0)  # Shape: (n_electrodes,)

        # Calculate noise power per channel
        snr_linear = 10 ** (snr__dB / 10)
        noise_power_per_channel = signal_power_per_channel / snr_linear
        noise_std_per_channel = np.sqrt(noise_power_per_channel)  # Shape: (n_electrodes,)

        if noise_kind == "gaussian":
            # Generate standard normal noise, then scale per channel
            noise = rng_global.normal(loc=0.0, scale=1.0, size=emg_array.shape)
            noise = noise * noise_std_per_channel[np.newaxis, :]
        else:
            # Colored noise (signal-shape preserving SNR semantics).
            # Per-channel call so each electrode hits the requested SNR
            # individually, matching the legacy gaussian path.
            from myogen.utils.emg_noise import generate_realistic_noise

            n_samples, n_channels = emg_array.shape
            fs_hz = float(self._sampling_frequency__Hz)
            noise = np.empty_like(emg_array)
            for ch in range(n_channels):
                ch_rng = np.random.default_rng(rng_global.integers(0, 2**31 - 1))
                noise[:, ch] = generate_realistic_noise(
                    n_samples,
                    fs_hz,
                    noise_rms=float(noise_std_per_channel[ch]),
                    spectral_slope=spectral_slope,
                    excess_kurtosis=excess_kurtosis,
                    powerline_hz=powerline_hz,
                    powerline_amplitude=powerline_amplitude,
                    powerline_harmonic_ratios=powerline_harmonic_ratios,
                    powerline_frequency_drift_hz=powerline_frequency_drift_hz,
                    powerline_amplitude_modulation_depth=powerline_amplitude_modulation_depth,
                    peak_hz=peak_hz,
                    analog_hpf_hz=analog_hpf_hz,
                    baseline_drift_rms_uv=baseline_drift_rms_uv,
                    baseline_drift_alpha=baseline_drift_alpha,
                    baseline_drift_low_hz=baseline_drift_low_hz,
                    baseline_drift_high_hz=baseline_drift_high_hz,
                    rng=ch_rng,
                )

        # Add noise
        noisy_emg = emg_array + noise

        # Create new AnalogSignal with noise
        noisy_segment.analogsignals.append(
            AnalogSignal(
                noisy_emg * emg_signal.units,
                t_start=emg_signal.t_start,
                sampling_rate=emg_signal.sampling_rate,
            )
        )

    # Store results privately
    self._noisy_intramuscular_emg__Block = noisy_block
    return noisy_block

SurfaceElectrodeArray

SurfaceElectrodeArray(num_rows: int, num_cols: int, inter_electrode_distance__mm: Quantity__mm, electrode_radius__mm: Quantity__mm, center_point__mm_deg: tuple[Quantity__mm, Quantity__deg] = (0.0 * mm, 0.0 * deg), bending_radius__mm: Quantity__mm = 0.0 * mm, rotation_angle__deg: Quantity__deg = 0.0 * deg, differentiation_mode: Literal['monopolar', 'bipolar_longitudinal', 'bipolar_transversal', 'laplacian'] = 'monopolar')

Surface electrode array for EMG recording.

Represents a grid of surface electrodes with configurable spacing, size, and differentiation modes.

Parameters:

Name Type Description Default
num_rows int

Number of rows in the electrode array

required
num_cols int

Number of columns in the electrode array

required
inter_electrode_distance__mm float

Inter-electrode distance in mm.

required
electrode_radius__mm float

Radius of the electrodes in mm

required
center_point__mm_deg tuple[float, float]

Position along z in mm and rotation around the muscle theta in degrees.

(0.0 * mm, 0.0 * deg)
bending_radius__mm float

Bending radius around which the electrode grid is bent. Usually this is equal to the radius of the muscle.

0.0 * mm
rotation_angle__deg float

Rotation angle of the electrodes in degrees. This is the angle between the electrode grid and the muscle surface.

0.0 * deg
differentiation_mode (monopolar, bipolar_longitudinal, bipolar_transversal, laplacian)

Differentiation mode. Default is monopolar.

"monopolar"

Attributes:

Name Type Description
pos_z ndarray

Longitudinal electrode positions in mm, shape (num_rows, num_cols). Available after class initialization via _create_electrode_grid().

pos_theta ndarray

Angular electrode positions in radians, shape (num_rows, num_cols). Available after class initialization via _create_electrode_grid().

electrode_positions tuple[ndarray, ndarray]

Complete electrode position arrays (pos_z, pos_theta). Available after class initialization via _create_electrode_grid().

num_electrodes int

Total number of electrodes (num_rows * num_cols).

num_channels int

Number of recording channels based on differentiation mode.

Source code in myogen/simulator/core/emg/electrodes.py
def __init__(
    self,
    num_rows: int,
    num_cols: int,
    inter_electrode_distance__mm: Quantity__mm,
    electrode_radius__mm: Quantity__mm,
    center_point__mm_deg: tuple[Quantity__mm, Quantity__deg] = (0.0 * pq.mm, 0.0 * pq.deg),
    bending_radius__mm: Quantity__mm = 0.0 * pq.mm,
    rotation_angle__deg: Quantity__deg = 0.0 * pq.deg,
    differentiation_mode: Literal[
        "monopolar", "bipolar_longitudinal", "bipolar_transversal", "laplacian"
    ] = "monopolar",
):
    # Immutable public arguments - never modify these
    self.num_rows = num_rows
    self.num_cols = num_cols
    self.center_point__mm_deg = center_point__mm_deg
    self.bending_radius__mm = bending_radius__mm
    self.rotation_angle__deg = rotation_angle__deg
    self.inter_electrode_distance__mm = inter_electrode_distance__mm
    self.electrode_radius__mm = electrode_radius__mm
    self.differentiation_mode = differentiation_mode

    # Private copies for internal modifications (extract magnitudes)
    self._num_rows = num_rows
    self._num_cols = num_cols
    self._center_point__mm_deg = (
        float(center_point__mm_deg[0].rescale(pq.mm).magnitude),
        float(center_point__mm_deg[1].rescale(pq.deg).magnitude),
    )
    self._bending_radius__mm = float(bending_radius__mm.rescale(pq.mm).magnitude)
    self._rotation_angle__deg = float(rotation_angle__deg.rescale(pq.deg).magnitude)
    self._inter_electrode_distance__mm = float(
        inter_electrode_distance__mm.rescale(pq.mm).magnitude
    )
    self._electrode_radius__mm = float(electrode_radius__mm.rescale(pq.mm).magnitude)
    self._differentiation_mode = differentiation_mode

    self.num_electrodes = num_rows * num_cols

    # Handle zero bending radius
    if self._bending_radius__mm == 0:
        self._bending_radius__mm = np.finfo(np.float32).eps

    # Set up channel configuration based on differentiation mode
    if differentiation_mode == "monopolar":
        self._num_channels = self.num_electrodes
    elif differentiation_mode in ["bipolar_longitudinal", "bipolar_transversal"]:
        # For bipolar, we lose one channel per dimension
        if differentiation_mode == "bipolar_longitudinal":
            self._num_channels = max(1, num_rows - 1) * num_cols
        else:  # bipolar_transversal
            self._num_channels = num_rows * max(1, num_cols - 1)
    elif differentiation_mode == "laplacian":
        # For Laplacian, we lose border electrodes
        self._num_channels = max(1, num_rows - 2) * max(1, num_cols - 2)
    else:
        self._num_channels = self.num_electrodes

    # Create electrode grid in local coordinate system
    self._create_electrode_grid()

pos_z property

pos_z: Quantity__mm

Longitudinal positions of electrodes in mm.

Returns:

Type Description
Quantity__mm

Array of shape (num_rows, num_cols) containing z-coordinates of each electrode position in mm.

Raises:

Type Description
AttributeError

If electrode grid has not been created. Run constructor first.

pos_theta property

pos_theta: Quantity__rad

Angular positions of electrodes in radians.

Returns:

Type Description
Quantity__rad

Array of shape (num_rows, num_cols) containing angular coordinates of each electrode position in radians.

Raises:

Type Description
AttributeError

If electrode grid has not been created. Run constructor first.

electrode_positions property

electrode_positions: tuple[Quantity__mm, Quantity__rad]

Complete electrode position arrays (z, theta) in physical coordinates.

Returns:

Type Description
tuple[Quantity__mm, Quantity__rad]

Tuple containing: - pos_z: Longitudinal positions in mm, shape (num_rows, num_cols) - pos_theta: Angular positions in radians, shape (num_rows, num_cols)

Raises:

Type Description
AttributeError

If electrode grid has not been created. Run constructor first.

num_channels property

num_channels: int

Number of recording channels based on differentiation mode.

Returns:

Type Description
int

Number of recording channels. Depends on differentiation_mode: - "monopolar": num_rows * num_cols - "bipolar_longitudinal": (num_rows - 1) * num_cols - "bipolar_transversal": num_rows * (num_cols - 1) - "laplacian": (num_rows - 2) * (num_cols - 2)

Raises:

Type Description
AttributeError

If channel count has not been calculated. Run constructor first.

get_H_sf

get_H_sf(ktheta_mesh_kzktheta: ndarray, kz_mesh_kzktheta: ndarray) -> ndarray | float

Get the spatial filter for the electrode array.

Parameters:

Name Type Description Default
ktheta_mesh_kzktheta ndarray

Angular spatial frequency mesh

required
kz_mesh_kzktheta ndarray

Longitudinal spatial frequency mesh

required

Returns:

Name Type Description
H_sf ndarray or float

Spatial filter for the specified differentiation mode

Source code in myogen/simulator/core/emg/electrodes.py
def get_H_sf(
    self, ktheta_mesh_kzktheta: np.ndarray, kz_mesh_kzktheta: np.ndarray
) -> np.ndarray | float:
    """
    Get the spatial filter for the electrode array.

    Parameters
    ----------
    ktheta_mesh_kzktheta : np.ndarray
        Angular spatial frequency mesh
    kz_mesh_kzktheta : np.ndarray
        Longitudinal spatial frequency mesh

    Returns
    -------
    H_sf : np.ndarray or float
        Spatial filter for the specified differentiation mode
    """
    if self.differentiation_mode == "monopolar":
        H_sf = 1.0

    elif self.differentiation_mode == "bipolar_longitudinal":
        # Differential along muscle fiber direction (z-axis)
        # Apply coordinate transformation for rotation (Farina 2004, eq 38)
        alpha_rad = self._rotation_angle__deg * np.pi / 180
        kz_new = ktheta_mesh_kzktheta / self._bending_radius__mm * np.sin(
            alpha_rad
        ) + kz_mesh_kzktheta * np.cos(alpha_rad)
        # Spatial filter for longitudinal differential (Farina 2004, eq 30)
        # Two electrodes at ±d_z/2: H_sf = exp(j·k'_z·d_z/2) - exp(-j·k'_z·d_z/2)
        half_ied = self._inter_electrode_distance__mm / 2
        H_sf = np.exp(1j * kz_new * half_ied) - np.exp(-1j * kz_new * half_ied)

    elif self.differentiation_mode == "bipolar_transversal":
        # Differential around muscle circumference (theta-axis)
        # Apply coordinate transformation for rotation (Farina 2004, eq 38)
        alpha_rad = self._rotation_angle__deg * np.pi / 180
        ktheta_new = ktheta_mesh_kzktheta * np.cos(
            alpha_rad
        ) - kz_mesh_kzktheta * self._bending_radius__mm * np.sin(alpha_rad)
        # Spatial filter for transversal differential (Farina 2004, eq 30-31)
        # d_theta = IED / R_ele in radians
        half_d_theta = self._inter_electrode_distance__mm / (2 * self._bending_radius__mm)
        H_sf = np.exp(1j * ktheta_new * half_d_theta) - np.exp(
            -1j * ktheta_new * half_d_theta
        )

    elif self.differentiation_mode == "laplacian":
        # Laplacian (second-order spatial differential)
        # Combination of longitudinal and transversal second derivatives
        alpha_rad = self._rotation_angle__deg * np.pi / 180
        kz_new = ktheta_mesh_kzktheta / self._bending_radius__mm * np.sin(
            alpha_rad
        ) + kz_mesh_kzktheta * np.cos(alpha_rad)
        ktheta_new = ktheta_mesh_kzktheta * np.cos(
            alpha_rad
        ) - kz_mesh_kzktheta * self._bending_radius__mm * np.sin(alpha_rad)

        # Laplacian approximation: second differences along both axes
        half_ied = self._inter_electrode_distance__mm / 2
        half_d_theta = self._inter_electrode_distance__mm / (2 * self._bending_radius__mm)
        k_total_sq = (kz_new * half_ied) ** 2 + (ktheta_new * half_d_theta) ** 2
        H_sf = -k_total_sq

    return H_sf

IntramuscularElectrodeArray

IntramuscularElectrodeArray(num_electrodes: int, inter_electrode_distance__mm: Quantity__mm = 0.5 * mm, position__mm: tuple[Quantity__mm, Quantity__mm, Quantity__mm] = (0.0 * mm, 0.0 * mm, 0.0 * mm), orientation__rad: tuple[Quantity__rad, Quantity__rad, Quantity__rad] = (0.0 * rad, 0.0 * rad, 0.0 * rad), differentiation_mode: Literal['consecutive', 'reference'] = 'consecutive', trajectory_distance__mm: Quantity__mm = 0.0 * mm, trajectory_steps: int = 1)

Intramuscular electrode array for EMG recording.

Represents a linear array of intramuscular electrodes (needle electrodes) with configurable spacing and differentiation modes.

Parameters:

Name Type Description Default
num_electrodes int

Number of electrodes in the array

required
inter_electrode_distance__mm float

Inter-electrode distance in mm

0.5
position__mm tuple[float, float, float]

Position of the electrode array center in mm (x, y, z coordinates)

(0.0, 0.0, 0.0)
orientation__rad tuple[float, float, float]

Orientation of the electrode array in radians (roll, pitch, yaw)

(0.0, 0.0, 0.0)
differentiation_mode Literal['consecutive', 'reference']

Differentiation mode for recording

"consecutive"
trajectory_distance__mm float

Distance for trajectory movement in mm

0.0
trajectory_steps int

Number of steps in the trajectory

1

Attributes:

Name Type Description
electrode_positions ndarray

Current electrode positions in 3D space, shape (n_nodes * num_electrodes, 3). Available after set_linear_trajectory() execution.

differential_matrix ndarray

Differential matrix for signal processing based on differentiation mode. Available after class initialization.

trajectory_transforms ndarray

Transformation matrices for trajectory movement, shape (n_nodes, 6). Available after set_linear_trajectory() execution.

initial_positions ndarray

Initial electrode positions after position/orientation setup, shape (num_electrodes, 3). Available after set_position() execution.

num_channels int

Number of recording channels based on differentiation mode. Available after class initialization.

num_points int

Alias for num_electrodes (compatibility).

n_nodes int

Number of trajectory nodes.

Source code in myogen/simulator/core/emg/electrodes.py
def __init__(
    self,
    num_electrodes: int,
    inter_electrode_distance__mm: Quantity__mm = 0.5 * pq.mm,
    position__mm: tuple[Quantity__mm, Quantity__mm, Quantity__mm] = (
        0.0 * pq.mm,
        0.0 * pq.mm,
        0.0 * pq.mm,
    ),
    orientation__rad: tuple[Quantity__rad, Quantity__rad, Quantity__rad] = (
        0.0 * pq.rad,
        0.0 * pq.rad,
        0.0 * pq.rad,
    ),
    differentiation_mode: Literal["consecutive", "reference"] = "consecutive",
    trajectory_distance__mm: Quantity__mm = 0.0 * pq.mm,
    trajectory_steps: int = 1,
):
    # Immutable public arguments - never modify these
    self.num_electrodes = num_electrodes
    self.inter_electrode_distance__mm = inter_electrode_distance__mm
    self.position__mm = position__mm
    self.orientation__rad = orientation__rad
    self.differentiation_mode = differentiation_mode
    self.trajectory_distance__mm = trajectory_distance__mm
    self.trajectory_steps = trajectory_steps

    # Private copies for internal modifications (extract magnitudes)
    self._num_electrodes = num_electrodes
    self._inter_electrode_distance__mm = float(
        inter_electrode_distance__mm.rescale(pq.mm).magnitude
    )
    self._position__mm = (
        float(position__mm[0].rescale(pq.mm).magnitude),
        float(position__mm[1].rescale(pq.mm).magnitude),
        float(position__mm[2].rescale(pq.mm).magnitude),
    )
    self._orientation__rad = (
        float(orientation__rad[0].rescale(pq.rad).magnitude),
        float(orientation__rad[1].rescale(pq.rad).magnitude),
        float(orientation__rad[2].rescale(pq.rad).magnitude),
    )
    self._differentiation_mode = differentiation_mode
    self._trajectory_distance__mm = float(trajectory_distance__mm.rescale(pq.mm).magnitude)
    self._trajectory_steps = trajectory_steps

    self.num_points = num_electrodes  # Alias for compatibility
    self.n_nodes = trajectory_steps

    self._pts_origin = np.concatenate(
        [
            np.zeros((self._num_electrodes, 2)),
            np.arange(self._num_electrodes)[..., None] * self._inter_electrode_distance__mm,
        ],
        axis=-1,
    )

    self._normal_origin = []
    self._normals_init = []
    self._normals = []

    match differentiation_mode:
        case "consecutive":
            eye_mat = np.eye(self._pts_origin.shape[0] - 1, self._pts_origin.shape[0])
            self._diff_mat = eye_mat - np.roll(eye_mat, shift=1, axis=1)
        case "reference":
            self._diff_mat = np.roll(
                np.eye(self._pts_origin.shape[0] - 1, self._pts_origin.shape[0]),
                shift=1,
                axis=1,
            )
            self._diff_mat[:, 0] = -1

    self._n_channels = self._diff_mat.shape[0]

    # Use public (Quantity) parameters for method calls
    self.set_position(position__mm=position__mm, orientation__rad=orientation__rad)
    self.set_linear_trajectory(
        distance__mm=trajectory_distance__mm, n_nodes=self._trajectory_steps
    )

electrode_positions property

electrode_positions: Quantity__mm

Current electrode positions in 3D space (mm).

Returns:

Type Description
Quantity__mm

Array of shape (n_nodes * num_electrodes, 3) containing x, y, z coordinates of each electrode position for all trajectory nodes, in mm.

Raises:

Type Description
AttributeError

If trajectory has not been calculated. Run set_linear_trajectory() first.

differential_matrix property

differential_matrix: ndarray

Differential matrix for signal processing based on differentiation mode.

Returns:

Type Description
ndarray

Differential matrix for applying spatial differentiation to recorded signals. Shape depends on differentiation mode and number of trajectory nodes.

Raises:

Type Description
AttributeError

If differential matrix has not been created. Run constructor first.

trajectory_transforms property

trajectory_transforms: ndarray

Transformation matrices for trajectory movement.

Returns:

Type Description
ndarray

Array of shape (n_nodes, 6) containing translation and rotation parameters for each trajectory node. First 3 columns are translations (x, y, z), last 3 columns are rotations (roll, pitch, yaw).

Raises:

Type Description
AttributeError

If trajectory has not been set. Run set_linear_trajectory() first.

initial_positions property

initial_positions: Quantity__mm

Initial electrode positions after position/orientation setup.

Returns:

Type Description
Quantity__mm

Array of shape (num_electrodes, 3) containing initial x, y, z coordinates of electrodes before trajectory movement is applied, in mm.

Raises:

Type Description
AttributeError

If initial positions have not been set. Run set_position() first.

num_channels property

num_channels: int

Number of recording channels based on differentiation mode.

Returns:

Type Description
int

Number of differential recording channels available from this electrode array.

Raises:

Type Description
AttributeError

If channel count has not been calculated. Run constructor first.

set_position

set_position(position__mm: tuple[Quantity__mm, Quantity__mm, Quantity__mm], orientation__rad: tuple[Quantity__rad, Quantity__rad, Quantity__rad]) -> None

Set the position and orientation of the intramuscular electrode array.

This method defines the spatial placement and angular orientation of the electrode array within the muscle volume. The array is first oriented according to the specified rotations and then translated to the target position.

Coordinate System: - x-axis: radial direction (outward from muscle center) - y-axis: circumferential direction (around muscle) - z-axis: longitudinal direction (along muscle fibers)

Rotation Order: Applied as: Roll (x) → Pitch (y) → Yaw (z) using Rodrigues rotation

Parameters:

Name Type Description Default
position__mm tuple[float, float, float]

Center position of the electrode array in mm (x, y, z coordinates). This defines where the array center is placed within the muscle.

required
orientation__rad tuple[float, float, float]

Orientation angles in radians (roll, pitch, yaw). - Roll: rotation around x-axis (radial tilt) - Pitch: rotation around y-axis (circumferential tilt) - Yaw: rotation around z-axis (longitudinal rotation)

required
Notes

Position and orientation changes affect all subsequent trajectory calculations. The electrode positions are recalculated based on the new transformation.

Examples:

>>> # Place array at muscle center with 45° yaw rotation
>>> array.set_position(
...     position__mm=(0.0, 0.0, 10.0),
...     orientation__rad=(0.0, 0.0, np.pi/4)
... )
See Also

set_linear_trajectory : Define trajectory movement parameters rodrigues_rot : Rodrigues rotation implementation

Source code in myogen/simulator/core/emg/electrodes.py
def set_position(
    self,
    position__mm: tuple[Quantity__mm, Quantity__mm, Quantity__mm],
    orientation__rad: tuple[Quantity__rad, Quantity__rad, Quantity__rad],
) -> None:
    """
    Set the position and orientation of the intramuscular electrode array.

    This method defines the spatial placement and angular orientation of the
    electrode array within the muscle volume. The array is first oriented
    according to the specified rotations and then translated to the target position.

    **Coordinate System:**
    - x-axis: radial direction (outward from muscle center)
    - y-axis: circumferential direction (around muscle)
    - z-axis: longitudinal direction (along muscle fibers)

    **Rotation Order:**
    Applied as: Roll (x) → Pitch (y) → Yaw (z) using Rodrigues rotation

    Parameters
    ----------
    position__mm : tuple[float, float, float]
        Center position of the electrode array in mm (x, y, z coordinates).
        This defines where the array center is placed within the muscle.
    orientation__rad : tuple[float, float, float]
        Orientation angles in radians (roll, pitch, yaw).
        - Roll: rotation around x-axis (radial tilt)
        - Pitch: rotation around y-axis (circumferential tilt)
        - Yaw: rotation around z-axis (longitudinal rotation)

    Notes
    -----
    Position and orientation changes affect all subsequent trajectory calculations.
    The electrode positions are recalculated based on the new transformation.

    Examples
    --------
    >>> # Place array at muscle center with 45° yaw rotation
    >>> array.set_position(
    ...     position__mm=(0.0, 0.0, 10.0),
    ...     orientation__rad=(0.0, 0.0, np.pi/4)
    ... )

    See Also
    --------
    set_linear_trajectory : Define trajectory movement parameters
    rodrigues_rot : Rodrigues rotation implementation
    """
    self._pts_init = np.copy(self._pts_origin)

    # Extract magnitude values for internal calculations
    orientation__rad_values = (
        float(orientation__rad[0].rescale(pq.rad).magnitude),
        float(orientation__rad[1].rescale(pq.rad).magnitude),
        float(orientation__rad[2].rescale(pq.rad).magnitude),
    )
    position__mm_values = np.array(
        [
            float(position__mm[0].rescale(pq.mm).magnitude),
            float(position__mm[1].rescale(pq.mm).magnitude),
            float(position__mm[2].rescale(pq.mm).magnitude),
        ]
    )

    self._pts_init = self.rodrigues_rot(self._pts_init, [1, 0, 0], orientation__rad_values[0])
    self._pts_init = self.rodrigues_rot(self._pts_init, [0, 1, 0], orientation__rad_values[1])
    self._pts_init = self.rodrigues_rot(self._pts_init, [0, 0, 1], orientation__rad_values[2])

    self._pts_init += np.matlib.repmat(position__mm_values[None], self._pts_init.shape[0], 1)
    self._pts = np.copy(self._pts_init)

rodrigues_rot

rodrigues_rot(v, k, theta)

Apply Rodrigues rotation to vectors around an arbitrary axis.

This method implements 3D rotation of points or vectors around an arbitrary axis using the Rodrigues rotation formula. It is used internally for electrode array positioning and trajectory calculations.

Mathematical Foundation: Based on Rodrigues' rotation formula for rotating a vector v around axis k by angle theta: v_rot = vcos(θ) + (k×v)sin(θ) + k(k·v)(1-cos(θ))

Parameters:

Name Type Description Default
v array_like

Vector(s) to rotate. Can be single vector (3,) or array of vectors (N, 3).

required
k array_like

Rotation axis vector (3,). Will be normalized internally.

required
theta float

Rotation angle in radians. Positive angles follow right-hand rule.

required

Returns:

Type Description
ndarray

Rotated vector(s) with same shape as input v.

Notes

Uses scipy.spatial.transform.Rotation for numerical stability and efficiency. The rotation axis k is automatically normalized to unit length.

Examples:

>>> # Rotate point 90° around z-axis
>>> point = np.array([1.0, 0.0, 0.0])
>>> rotated = array.rodrigues_rot(point, [0, 0, 1], np.pi/2)
>>> # Result: approximately [0, 1, 0]
Source code in myogen/simulator/core/emg/electrodes.py
def rodrigues_rot(self, v, k, theta):
    """
    Apply Rodrigues rotation to vectors around an arbitrary axis.

    This method implements 3D rotation of points or vectors around an arbitrary
    axis using the Rodrigues rotation formula. It is used internally for
    electrode array positioning and trajectory calculations.

    **Mathematical Foundation:**
    Based on Rodrigues' rotation formula for rotating a vector v around
    axis k by angle theta: v_rot = v*cos(θ) + (k×v)*sin(θ) + k*(k·v)*(1-cos(θ))

    Parameters
    ----------
    v : array_like
        Vector(s) to rotate. Can be single vector (3,) or array of vectors (N, 3).
    k : array_like
        Rotation axis vector (3,). Will be normalized internally.
    theta : float
        Rotation angle in radians. Positive angles follow right-hand rule.

    Returns
    -------
    np.ndarray
        Rotated vector(s) with same shape as input v.

    Notes
    -----
    Uses scipy.spatial.transform.Rotation for numerical stability and efficiency.
    The rotation axis k is automatically normalized to unit length.

    Examples
    --------
    >>> # Rotate point 90° around z-axis
    >>> point = np.array([1.0, 0.0, 0.0])
    >>> rotated = array.rodrigues_rot(point, [0, 0, 1], np.pi/2)
    >>> # Result: approximately [0, 1, 0]
    """
    v = np.array(v.copy(), dtype=float)
    k = np.array(k.copy(), dtype=float)
    k = k / np.linalg.norm(k)  # normalize axis

    r = R.from_rotvec(k * theta)  # Create rotation from axis-angle
    return r.apply(v)  # Rotates v (works with (3,), (N, 3))

set_linear_trajectory

set_linear_trajectory(distance__mm: Quantity__mm, n_nodes: int | None = None) -> None

Configure linear trajectory movement for the electrode array.

This method sets up a linear movement path for the electrode array, simulating needle insertion or withdrawal. The trajectory is discretized into nodes for temporal interpolation during EMG simulation.

Trajectory Properties: - Direction: Along the array's longitudinal axis (z-direction in local coordinates) - Movement: Linear progression from start to end position - Discretization: Evenly spaced nodes for smooth interpolation - Default step size: 0.5mm if n_nodes not specified

Parameters:

Name Type Description Default
distance__mm float

Total trajectory distance in mm. Positive values move in the positive z-direction of the array's local coordinate system.

required
n_nodes int

Number of discrete trajectory nodes. If None, automatically calculated as max(ceil(distance/0.5), 1) for 0.5mm steps.

None
Notes

The trajectory is applied after position and orientation transformations. All trajectory transforms are calculated in the array's oriented coordinate system.

Examples:

>>> # Set up 10mm insertion with default step size (~0.5mm)
>>> array.set_linear_trajectory(distance__mm=10.0)
>>> # Set up 5mm trajectory with specific number of nodes
>>> array.set_linear_trajectory(distance__mm=5.0, n_nodes=20)
See Also

calc_observation_points : Calculate electrode positions along trajectory traj_mixing_mat : Generate mixing matrices for trajectory interpolation

Source code in myogen/simulator/core/emg/electrodes.py
def set_linear_trajectory(self, distance__mm: Quantity__mm, n_nodes: int | None = None) -> None:
    """
    Configure linear trajectory movement for the electrode array.

    This method sets up a linear movement path for the electrode array,
    simulating needle insertion or withdrawal. The trajectory is discretized
    into nodes for temporal interpolation during EMG simulation.

    **Trajectory Properties:**
    - Direction: Along the array's longitudinal axis (z-direction in local coordinates)
    - Movement: Linear progression from start to end position
    - Discretization: Evenly spaced nodes for smooth interpolation
    - Default step size: 0.5mm if n_nodes not specified

    Parameters
    ----------
    distance__mm : float
        Total trajectory distance in mm. Positive values move in the
        positive z-direction of the array's local coordinate system.
    n_nodes : int, optional
        Number of discrete trajectory nodes. If None, automatically
        calculated as max(ceil(distance/0.5), 1) for 0.5mm steps.

    Notes
    -----
    The trajectory is applied after position and orientation transformations.
    All trajectory transforms are calculated in the array's oriented coordinate system.

    Examples
    --------
    >>> # Set up 10mm insertion with default step size (~0.5mm)
    >>> array.set_linear_trajectory(distance__mm=10.0)

    >>> # Set up 5mm trajectory with specific number of nodes
    >>> array.set_linear_trajectory(distance__mm=5.0, n_nodes=20)

    See Also
    --------
    calc_observation_points : Calculate electrode positions along trajectory
    traj_mixing_mat : Generate mixing matrices for trajectory interpolation
    """
    # Extract magnitude value for internal calculations
    distance__mm_value = float(distance__mm.rescale(pq.mm).magnitude)

    if n_nodes is None:
        n_nodes = max(np.ceil(distance__mm_value / 0.5), 1)

    self.n_nodes = n_nodes
    self._trajectory_step = distance__mm_value / self.n_nodes

    self.traj_transforms = np.linspace(start=0, stop=distance__mm_value, num=self.n_nodes)
    self.traj_transforms = np.hstack(
        [
            np.zeros((max(self.traj_transforms.shape), 2)),
            self.traj_transforms.reshape(-1, 1),
            np.zeros((max(self.traj_transforms.shape), 3)),
        ]
    )

    # Use private orientation values (already extracted as floats)
    self.traj_transforms[:, :3] = self.rodrigues_rot(
        self.traj_transforms[:, :3], [1, 0, 0], self._orientation__rad[0]
    )
    self.traj_transforms[:, :3] = self.rodrigues_rot(
        self.traj_transforms[:, :3], [0, 1, 0], self._orientation__rad[1]
    )
    self.traj_transforms[:, :3] = self.rodrigues_rot(
        self.traj_transforms[:, :3], [0, 0, 1], self._orientation__rad[2]
    )

    self.calc_observation_points()

traj_mixing_fun

traj_mixing_fun(t, n_nodes, node) -> float

Compute mixing weight for a specific trajectory node at given time.

This function calculates the interpolation weight for a trajectory node based on the current time/position along the trajectory. Uses triangular weighting where nodes closer to the current time get higher weights.

Parameters:

Name Type Description Default
t float

Current normalized time or position in trajectory (0.0 to 1.0).

required
n_nodes int

Total number of nodes in the trajectory.

required
node int or array_like

Node index(es) for which to calculate mixing weights. Can be scalar or array of node indices.

required

Returns:

Type Description
float or ndarray

Mixing weight(s) for the specified node(s) at time t. Returns 0 for distant nodes, max weight 1 for closest node.

Source code in myogen/simulator/core/emg/electrodes.py
def traj_mixing_fun(self, t, n_nodes, node) -> float:
    """
    Compute mixing weight for a specific trajectory node at given time.

    This function calculates the interpolation weight for a trajectory node based
    on the current time/position along the trajectory. Uses triangular weighting
    where nodes closer to the current time get higher weights.

    Parameters
    ----------
    t : float
        Current normalized time or position in trajectory (0.0 to 1.0).
    n_nodes : int
        Total number of nodes in the trajectory.
    node : int or array_like
        Node index(es) for which to calculate mixing weights.
        Can be scalar or array of node indices.

    Returns
    -------
    float or np.ndarray
        Mixing weight(s) for the specified node(s) at time t.
        Returns 0 for distant nodes, max weight 1 for closest node.
    """
    eps = np.finfo(float).eps
    return np.maximum(
        0,
        1 - (n_nodes - 1) * np.abs(t - (node - 1) / (n_nodes - 1 + eps)),
    )

traj_mixing_mat

traj_mixing_mat(t, n_nodes, n_channels) -> ndarray

Generate mixing matrix for trajectory interpolation during EMG simulation.

This method creates a diagonal mixing matrix that weights the contribution of different trajectory nodes during temporal interpolation. The matrix enables smooth transitions between electrode positions as the array moves along its trajectory during needle insertion or withdrawal.

Interpolation Strategy: - Linear interpolation between adjacent trajectory nodes - Weights based on distance from current time/position to node positions - Diagonal matrix structure for efficient computation - Smooth transitions avoid discontinuities in EMG signals

Parameters:

Name Type Description Default
t float

Current normalized time or position in trajectory (0.0 to 1.0). 0.0 corresponds to trajectory start, 1.0 to trajectory end.

required
n_nodes int

Total number of trajectory nodes for interpolation.

required
n_channels int

Number of recording channels in the electrode array. Depends on differentiation mode and electrode count.

required

Returns:

Type Description
ndarray

Diagonal mixing matrix with shape (n_nodes * n_channels, n_nodes * n_channels). Diagonal elements contain repeated mixing weights for each trajectory node, with each node's weight repeated n_channels times.

Notes

The mixing matrix enables temporal interpolation of EMG signals recorded at different trajectory positions. Higher weights are given to nodes closer to the current time/position parameter.

Examples:

>>> # Get mixing weights for mid-trajectory position
>>> mix_mat = array.traj_mixing_mat(t=0.5, n_nodes=10, n_channels=4)
>>> # Matrix will weight middle nodes more heavily
See Also

traj_mixing_fun : Individual node mixing function set_linear_trajectory : Configure trajectory parameters

Source code in myogen/simulator/core/emg/electrodes.py
def traj_mixing_mat(self, t, n_nodes, n_channels) -> np.ndarray:
    """
    Generate mixing matrix for trajectory interpolation during EMG simulation.

    This method creates a diagonal mixing matrix that weights the contribution of
    different trajectory nodes during temporal interpolation. The matrix enables
    smooth transitions between electrode positions as the array moves along its
    trajectory during needle insertion or withdrawal.

    **Interpolation Strategy:**
    - Linear interpolation between adjacent trajectory nodes
    - Weights based on distance from current time/position to node positions
    - Diagonal matrix structure for efficient computation
    - Smooth transitions avoid discontinuities in EMG signals

    Parameters
    ----------
    t : float
        Current normalized time or position in trajectory (0.0 to 1.0).
        0.0 corresponds to trajectory start, 1.0 to trajectory end.
    n_nodes : int
        Total number of trajectory nodes for interpolation.
    n_channels : int
        Number of recording channels in the electrode array.
        Depends on differentiation mode and electrode count.

    Returns
    -------
    np.ndarray
        Diagonal mixing matrix with shape (n_nodes * n_channels, n_nodes * n_channels).
        Diagonal elements contain repeated mixing weights for each trajectory node,
        with each node's weight repeated n_channels times.

    Notes
    -----
    The mixing matrix enables temporal interpolation of EMG signals recorded
    at different trajectory positions. Higher weights are given to nodes
    closer to the current time/position parameter.

    Examples
    --------
    >>> # Get mixing weights for mid-trajectory position
    >>> mix_mat = array.traj_mixing_mat(t=0.5, n_nodes=10, n_channels=4)
    >>> # Matrix will weight middle nodes more heavily

    See Also
    --------
    traj_mixing_fun : Individual node mixing function
    set_linear_trajectory : Configure trajectory parameters
    """

    return np.diag(
        np.repeat(
            self.traj_mixing_fun(t, n_nodes, np.arange(1, n_nodes + 1)),
            n_channels,
        )
    )

Proprioception

SpindleModel

SpindleModel(simulation_time__ms: Quantity__ms, time_step__ms: Quantity__ms, spindle_parameters: dict[str, Any])

API wrapper for the muscle spindle model.

This class provides an intuitive interface for creating muscle spindle models with user-friendly parameter names that are internally mapped to the correct format expected by the underlying Spindle implementation.

The muscle spindle is a proprioceptive sensory organ that detects changes in muscle length and velocity, providing feedback for motor control.

Parameters:

Name Type Description Default
simulation_time__ms Quantity__ms

Total simulation time in milliseconds

required
time_step__ms Quantity__ms

Integration time step in milliseconds

required
spindle_parameters dict[str, Any]

Dictionary containing spindle model parameters

required
Source code in myogen/simulator/neuron/proprioception/spindle.py
def __init__(
    self,
    simulation_time__ms: Quantity__ms,
    time_step__ms: Quantity__ms,
    spindle_parameters: dict[str, Any],
):
    self.simulation_time__ms = simulation_time__ms
    self.time_step__ms = time_step__ms
    self.spindle_parameters = spindle_parameters.copy()

    # Private working copies for internal use
    self._simulation_time__ms = simulation_time__ms
    self._time_step__ms = time_step__ms
    self._spindle_parameters = spindle_parameters.copy()

    # Validate inputs
    self._validate_parameters()

    # Create the underlying Spindle model
    self._spindle_model = self._create_spindle_model()

primary_afferent_firing__Hz property

primary_afferent_firing__Hz: ndarray

Get primary afferent (Ia) firing rate time series in Hz.

secondary_afferent_firing__Hz property

secondary_afferent_firing__Hz: ndarray

Get secondary afferent (II) firing rate time series in Hz.

bag1_activation property

bag1_activation: ndarray

Get Bag1 fiber activation time series.

bag2_activation property

bag2_activation: ndarray

Get Bag2 fiber activation time series.

chain_activation property

chain_activation: ndarray

Get Chain fiber activation time series.

intrafusal_tensions property

intrafusal_tensions: ndarray

Get intrafusal fiber tensions matrix (3 × time_points) [Bag1, Bag2, Chain].

time_vector property

time_vector: ndarray

Get simulation time vector in milliseconds.

integrate

integrate(muscle_length__L0: float, muscle_velocity__L0_per_s: float, muscle_acceleration__L0_per_s2: float, gamma_dynamic_drive__Hz: float, gamma_static_drive__Hz: float) -> tuple[float, float]

Integrate the spindle model for one time step.

Parameters:

Name Type Description Default
muscle_length__L0 float

Current muscle length normalized to L0

required
muscle_velocity__L0_per_s float

Current muscle velocity in L0/s

required
muscle_acceleration__L0_per_s2 float

Current muscle acceleration in L0/s²

required
gamma_dynamic_drive__Hz float

Gamma dynamic motor neuron drive frequency in Hz

required
gamma_static_drive__Hz float

Gamma static motor neuron drive frequency in Hz

required

Returns:

Type Description
tuple[float, float]

Primary afferent (Ia) and secondary afferent (II) firing rates in Hz

Source code in myogen/simulator/neuron/proprioception/spindle.py
def integrate(
    self,
    muscle_length__L0: float,
    muscle_velocity__L0_per_s: float,
    muscle_acceleration__L0_per_s2: float,
    gamma_dynamic_drive__Hz: float,
    gamma_static_drive__Hz: float,
) -> tuple[float, float]:
    """
    Integrate the spindle model for one time step.

    Parameters
    ----------
    muscle_length__L0 : float
        Current muscle length normalized to L0
    muscle_velocity__L0_per_s : float
        Current muscle velocity in L0/s
    muscle_acceleration__L0_per_s2 : float
        Current muscle acceleration in L0/s²
    gamma_dynamic_drive__Hz : float
        Gamma dynamic motor neuron drive frequency in Hz
    gamma_static_drive__Hz : float
        Gamma static motor neuron drive frequency in Hz

    Returns
    -------
    tuple[float, float]
        Primary afferent (Ia) and secondary afferent (II) firing rates in Hz
    """
    return self._spindle_model.integrate(
        muscle_length__L0,
        muscle_velocity__L0_per_s,
        muscle_acceleration__L0_per_s2,
        gamma_dynamic_drive__Hz,
        gamma_static_drive__Hz,
    )

create_default_spindle_parameters staticmethod

create_default_spindle_parameters(species: str = 'human', deafferent_ia: bool = False, deafferent_ii: bool = False) -> dict[str, Any]

Create default spindle parameter dictionary.

Parameters:

Name Type Description Default
species str

Species type ("human" or "cat"), by default "human"

'human'
deafferent_ia bool

Whether to simulate Ia afferent deafferentation, by default False

False
deafferent_ii bool

Whether to simulate II afferent deafferentation, by default False

False

Returns:

Type Description
dict[str, Any]

Dictionary of spindle parameters with detailed explanations

Raises:

Type Description
ValueError

If species is not recognized

Source code in myogen/simulator/neuron/proprioception/spindle.py
@staticmethod
def create_default_spindle_parameters(
    species: str = "human", deafferent_ia: bool = False, deafferent_ii: bool = False
) -> dict[str, Any]:
    """
    Create default spindle parameter dictionary.

    Parameters
    ----------
    species : str, optional
        Species type ("human" or "cat"), by default "human"
    deafferent_ia : bool, optional
        Whether to simulate Ia afferent deafferentation, by default False
    deafferent_ii : bool, optional
        Whether to simulate II afferent deafferentation, by default False

    Returns
    -------
    dict[str, Any]
        Dictionary of spindle parameters with detailed explanations

    Raises
    ------
    ValueError
        If species is not recognized
    """
    # Base spindle parameters (Mileusnic et al., 2006)
    spindle_params = {
        # Fusimotor activation parameters
        "fBag1": 60,  # Fusimotor frequency to activation constant for Bag1 [Hz]
        "fBag2": 60,  # Fusimotor frequency to activation constant for Bag2 [Hz]
        "fChain": 90,  # Fusimotor frequency to activation constant for Chain [Hz]
        "P": 2,  # Fusimotor frequency to activation power constant
        # Force generation coefficients
        "G1": 0.0289,  # Dynamic fusimotor input force generation coef [FU]
        "G2": 0.0636,  # Static fusimotor input force generation coef [FU]
        "G2Chain": 0.0954,  # Static fusimotor input force gen coef for Chain [FU]
        # Sensory Region (SR) mechanical parameters
        "K_SR": 10.4649,  # SR spring constant [FU/L0] - detects length changes
        "L0_SR": 0.04,  # SR rest length [L0] - baseline length
        "LN_SR": 0.0423,  # SR threshold length [L0] - minimum for activation
        # Polar Region (PR) mechanical parameters
        "K_PR": 0.15,  # PR spring constant [FU/L0] - contractile region
        "L0_PR": 0.76,  # PR rest length [L0] - baseline contractile length
        "LN_PR": 0.89,  # PR threshold length [L0] - minimum for activation
        # Intrafusal fiber mechanical properties
        "M": 0.0002,  # Intrafusal fiber mass [FU/(L0/s²)] - inertial component
        # Passive damping coefficients [FU/(L0/s)] - baseline viscosity
        "b0Bag1": 0.0605,  # Bag1 passive damping
        "b0Bag2": 0.0822,  # Bag2 passive damping
        "b0Chain": 0.0822,  # Chain passive damping
        # Fusimotor-dependent damping coefficients [FU/(L0/s)]
        "b1Bag1": 0.2592,  # Dynamic fusimotor damping for Bag1
        "b2Bag2": -0.0460,  # Static fusimotor damping for Bag2
        "b2Chain": -0.0690,  # Static fusimotor damping for Chain
        # Force-velocity relationship parameters
        "a": 0.3,  # Nonlinear velocity dependence power constant
        "C_L": 1,  # Lengthening coefficient of asymmetry in F-V curve
        "C_S": 0.42,  # Shortening coefficient of asymmetry in F-V curve
        "R": 0.46,  # Fascicle length where force production is zero [L0]
        # Afferent firing properties
        "X": 0.7,  # Secondary afferent percentage on sensory region [0-1]
        "Lsec": 0.04,  # Secondary afferent rest length [L0]
        "S": 0.156,  # Occlusion factor for primary afferent interactions
        # Temporal dynamics (low-pass filtering)
        "tau1": 0.149,  # Bag1 activation time constant [s] - fast dynamics
        "tau2": 0.205,  # Bag2 activation time constant [s] - slow dynamics
        # Afferent sensitivity gains [Hz/L0] - firing rate per unit stretch
        "gBag1": 6500,  # Bag1 contribution to primary afferent (Ia)
        "gBag2A1": 3250,  # Bag2 contribution to primary afferent (Ia)
        "gChainA1": 3250,  # Chain contribution to primary afferent (Ia)
        "gBag2A2": 3500,  # Bag2 contribution to secondary afferent (II)
        "gChainA2": 3500,  # Chain contribution to secondary afferent (II)
    }

    # Species-specific and deafferentation modifications
    if species == "human":
        if not deafferent_ii and not deafferent_ia:
            # Normal human spindle (Case 1, Elias thesis pg 66)
            pass  # Use default values above

        elif deafferent_ii and not deafferent_ia:
            # Human with Type II deafferentation (Case 2, Elias thesis pg 66)
            spindle_params.update(
                {
                    "gBag1": 7000,  # Enhanced Bag1 sensitivity
                    "gBag2A1": 3750,  # Enhanced Bag2 primary sensitivity
                    "gChainA1": 3750,  # Enhanced Chain primary sensitivity
                    "gBag2A2": 0,  # No Bag2 secondary afferents
                    "gChainA2": 0,  # No Chain secondary afferents
                }
            )

        elif not deafferent_ii and deafferent_ia:
            # Human with Ia deafferentation (Case 3, Elias thesis pg 66)
            spindle_params.update(
                {
                    "gBag1": 0,  # No Bag1 primary afferents
                    "gBag2A1": 0,  # No Bag2 primary afferents
                    "gChainA1": 0,  # No Chain primary afferents
                    "gBag2A2": 4500,  # Enhanced Bag2 secondary sensitivity
                    "gChainA2": 4500,  # Enhanced Chain secondary sensitivity
                }
            )

    elif species == "cat":
        # Cat spindle parameters (original Mileusnic values)
        spindle_params.update(
            {
                "gBag1": 20000,  # Higher sensitivity in cat
                "gBag2A1": 10000,  # Higher Bag2 primary sensitivity
                "gChainA1": 10000,  # Higher Chain primary sensitivity
                "gBag2A2": 7250,  # Higher Bag2 secondary sensitivity
                "gChainA2": 7250,  # Higher Chain secondary sensitivity
            }
        )

    else:
        raise ValueError(f"Unknown species: {species}. Use 'human' or 'cat'.")

    return spindle_params

GolgiTendonOrganModel

GolgiTendonOrganModel(simulation_time__ms: Quantity__ms, time_step__ms: Quantity__ms, gto_parameters: dict[str, Any])

API wrapper for the Golgi Tendon Organ (GTO) model.

This class provides an intuitive interface for creating GTO models with user-friendly parameter names that are internally mapped to the correct format expected by the underlying GTO implementation.

The Golgi Tendon Organ is a proprioceptive sensory organ located at the muscle-tendon junction that detects muscle force/tension and provides feedback for motor control and protection against excessive forces.

The model is based on Lin & Crago (2002) and implements a logarithmic force-to-firing relationship with digital filtering for realistic afferent discharge patterns.

Parameters:

Name Type Description Default
simulation_time__ms Quantity__ms

Total simulation time in milliseconds

required
time_step__ms Quantity__ms

Integration time step in milliseconds

required
gto_parameters dict[str, Any]

Dictionary containing GTO model parameters

required
Source code in myogen/simulator/neuron/proprioception/golgi.py
def __init__(
    self,
    simulation_time__ms: Quantity__ms,
    time_step__ms: Quantity__ms,
    gto_parameters: dict[str, Any],
):
    # Store original parameters (immutable)
    self.simulation_time__ms = simulation_time__ms
    self.time_step__ms = time_step__ms
    self.gto_parameters = gto_parameters.copy()

    # Private working copies for internal use
    self._simulation_time__ms = simulation_time__ms
    self._time_step__ms = time_step__ms
    self._gto_parameters = gto_parameters.copy()

    # Validate inputs
    self._validate_parameters()

    # Create the underlying GTO model
    self._gto_model = self._create_gto_model()

ib_afferent_firing__Hz property

ib_afferent_firing__Hz: ndarray

Get Ib afferent firing rate time series in Hz.

integrate

integrate(muscle_force__N: float) -> float

Integrate the GTO model for one time step.

Parameters:

Name Type Description Default
muscle_force__N float

Current muscle force in Newtons

required

Returns:

Type Description
float

Ib afferent firing rate in Hz

Source code in myogen/simulator/neuron/proprioception/golgi.py
def integrate(self, muscle_force__N: float) -> float:
    """
    Integrate the GTO model for one time step.

    Parameters
    ----------
    muscle_force__N : float
        Current muscle force in Newtons

    Returns
    -------
    float
        Ib afferent firing rate in Hz
    """
    return self._gto_model.integrate(muscle_force__N)

create_default_gto_parameters staticmethod

create_default_gto_parameters() -> dict[str, Any]

Create default Golgi Tendon Organ parameter dictionary.

The GTO model uses a logarithmic force-to-firing relationship: firing_rate = G1 * log(force/G2 + 1)

This is followed by digital filtering to create realistic temporal dynamics in the afferent discharge pattern.

Returns:

Type Description
dict[str, Any]

Dictionary of GTO parameters with detailed explanations

Notes

Model based on: - Lin & Crago (2002): Mathematical model framework - Aniss et al. (1990b): Human GTO physiological data - Elias PhD thesis (pg 83): Implementation details

The logarithmic relationship captures the GTO's ability to encode force over a wide dynamic range, from threshold detection of small forces to saturation at high forces, providing force feedback for motor control and protective reflexes.

Source code in myogen/simulator/neuron/proprioception/golgi.py
@staticmethod
def create_default_gto_parameters() -> dict[str, Any]:
    """
    Create default Golgi Tendon Organ parameter dictionary.

    The GTO model uses a logarithmic force-to-firing relationship:
    firing_rate = G1 * log(force/G2 + 1)

    This is followed by digital filtering to create realistic temporal
    dynamics in the afferent discharge pattern.

    Returns
    -------
    dict[str, Any]
        Dictionary of GTO parameters with detailed explanations

    Notes
    -----
    Model based on:
    - Lin & Crago (2002): Mathematical model framework
    - Aniss et al. (1990b): Human GTO physiological data
    - Elias PhD thesis (pg 83): Implementation details

    The logarithmic relationship captures the GTO's ability to encode
    force over a wide dynamic range, from threshold detection of small
    forces to saturation at high forces, providing force feedback for
    motor control and protective reflexes.
    """
    return {
        # Force-to-firing relationship parameters
        # Firing rate [Hz] = G1 * log(force[N]/G2 + 1)
        "G1": 40,  # Logarithmic gain coefficient [Hz]
        # Controls the sensitivity of force-to-firing conversion
        # Higher G1 = more sensitive to force changes
        # Typical range: 30-60 Hz for different muscles
        # Source: Lin & Crago (2002), Elias thesis uses 40 Hz
        "G2": 4,  # Force scaling coefficient [N]
        # Sets the force level for logarithmic scaling
        # Lower G2 = more sensitive to low forces
        # Higher G2 = requires higher forces for activation
        # Typical range: 2-8 N depending on muscle strength
        # Source: Calibrated for human muscle force ranges
    }

create_gto_parameters_for_muscle staticmethod

create_gto_parameters_for_muscle(muscle_type: str = 'FDI') -> dict[str, Any]

Create GTO parameters optimized for specific muscle types.

Parameters:

Name Type Description Default
muscle_type str

Type of muscle ("FDI", "Sol", "generic"), by default "FDI"

'FDI'

Returns:

Type Description
dict[str, Any]

Dictionary of muscle-specific GTO parameters

Notes

Different muscles have different force production capabilities and thus require different GTO sensitivity parameters:

  • FDI (First Dorsal Interosseous): Small hand muscle, low forces
  • Sol (Soleus): Large calf muscle, high forces
  • Generic: General-purpose parameters
Source code in myogen/simulator/neuron/proprioception/golgi.py
@staticmethod
def create_gto_parameters_for_muscle(muscle_type: str = "FDI") -> dict[str, Any]:
    """
    Create GTO parameters optimized for specific muscle types.

    Parameters
    ----------
    muscle_type : str, optional
        Type of muscle ("FDI", "Sol", "generic"), by default "FDI"

    Returns
    -------
    dict[str, Any]
        Dictionary of muscle-specific GTO parameters

    Notes
    -----
    Different muscles have different force production capabilities
    and thus require different GTO sensitivity parameters:

    - FDI (First Dorsal Interosseous): Small hand muscle, low forces
    - Sol (Soleus): Large calf muscle, high forces
    - Generic: General-purpose parameters
    """
    if muscle_type == "FDI":
        # Small hand muscle - more sensitive to small forces
        return {
            "G1": 45,  # Higher sensitivity for small force detection
            "G2": 2,  # Lower threshold for activation at small forces
        }

    elif muscle_type == "Sol":
        # Large calf muscle - less sensitive, handles high forces
        return {
            "G1": 35,  # Lower sensitivity appropriate for large forces
            "G2": 8,  # Higher threshold matching muscle's force capacity
        }

    elif muscle_type == "generic":
        # General-purpose parameters
        return create_default_gto_parameters()

    else:
        raise ValueError(f"Unknown muscle type: {muscle_type}. Use 'FDI', 'Sol', or 'generic'.")

JointDynamics

JointDynamics(inertia__kg_m2: float, damping__Nm_s_per_rad: float, stiffness__Nm_per_rad: float = 0.0, initial_angle__deg: float = 0.0, initial_velocity__deg_per_s: float = 0.0)

Joint dynamics integrator for closed-loop neuromechanical control.

This class implements second-order joint dynamics using the equation: I⋅α = τ - B⋅ω - K⋅θ where α = angular acceleration, τ = torque, ω = angular velocity, θ = joint angle, I = inertia, B = damping, K = stiffness.

Parameters:

Name Type Description Default
inertia__kg_m2 float

Joint rotational inertia in kg⋅m².

required
damping__Nm_s_per_rad float

Joint viscous damping coefficient in N⋅m⋅s/rad. Controls velocity-dependent resistance.

required
stiffness__Nm_per_rad float

Joint elastic stiffness in N⋅m/rad. Set to 0 for passive joints, >0 for spring-loaded joints.

0.0
initial_angle__deg float

Initial joint angle in degrees.

0.0
initial_velocity__deg_per_s float

Initial angular velocity in degrees per second.

0.0

Attributes:

Name Type Description
angle__rad float

Current joint angle in radians.

velocity__rad_per_s float

Current angular velocity in rad/s.

angle__deg float

Current joint angle in degrees (computed property).

Notes

JointDynamics uses plain float parameters (not quantities.Quantity) for the inertial/damping coefficients because no canonical Quantity__kg_m2, Quantity__Nm_s_per_rad, or Quantity__deg_per_s type aliases exist in myogen.utils.types. Values are interpreted as:

  • inertia__kg_m2 — rotational inertia in kg·m²
  • damping__Nm_s_per_rad — viscous damping in N·m·s/rad
  • initial_angle__deg — starting joint angle in degrees
  • initial_velocity__deg_per_s — starting angular velocity in degrees/second
Source code in myogen/simulator/neuron/joint_dynamics.py
def __init__(
    self,
    inertia__kg_m2: float,
    damping__Nm_s_per_rad: float,
    stiffness__Nm_per_rad: float = 0.0,
    initial_angle__deg: float = 0.0,
    initial_velocity__deg_per_s: float = 0.0,
) -> None:
    # Input validation
    if inertia__kg_m2 <= 0:
        raise ValueError(
            f"inertia__kg_m2 must be positive, got {inertia__kg_m2}. "
            "Typical values: finger=0.001-0.01, elbow=0.1-0.5, knee=1.0-5.0 kg⋅m²"
        )

    if damping__Nm_s_per_rad < 0:
        raise ValueError(
            f"damping__Nm_s_per_rad must be non-negative, got {damping__Nm_s_per_rad}. "
            "Typical values: 0.001-0.1 N⋅m⋅s/rad"
        )

    if stiffness__Nm_per_rad < 0:
        raise ValueError(
            f"stiffness__Nm_per_rad must be non-negative, got {stiffness__Nm_per_rad}. "
            "Use 0 for passive joints, >0 for spring-loaded joints"
        )

    # Immutable public parameters
    self.inertia__kg_m2 = inertia__kg_m2
    self.damping__Nm_s_per_rad = damping__Nm_s_per_rad
    self.stiffness__Nm_per_rad = stiffness__Nm_per_rad
    self.initial_angle__deg = initial_angle__deg
    self.initial_velocity__deg_per_s = initial_velocity__deg_per_s

    # Private working copies
    self._inertia = inertia__kg_m2
    self._damping = damping__Nm_s_per_rad
    self._stiffness = stiffness__Nm_per_rad

    # Joint state variables
    self.angle__rad = np.radians(initial_angle__deg)
    self.velocity__rad_per_s = np.radians(initial_velocity__deg_per_s)

angle__deg property

angle__deg: float

Current joint angle in degrees.

integrate

integrate(torque__Nm: float, dt__s: float) -> tuple[float, float]

Integrate joint dynamics for one time step.

Parameters:

Name Type Description Default
torque__Nm float

Applied muscle torque in N⋅m.

required
dt__s float

Integration time step in seconds.

required

Returns:

Type Description
tuple[float, float]

Updated (angle__deg, velocity__deg_per_s).

Source code in myogen/simulator/neuron/joint_dynamics.py
def integrate(self, torque__Nm: float, dt__s: float) -> tuple[float, float]:
    """
    Integrate joint dynamics for one time step.

    Parameters
    ----------
    torque__Nm : float
        Applied muscle torque in N⋅m.
    dt__s : float
        Integration time step in seconds.

    Returns
    -------
    tuple[float, float]
        Updated (angle__deg, velocity__deg_per_s).
    """
    # Second-order dynamics: I⋅α = τ - B⋅ω - K⋅θ
    spring_torque = -self._stiffness * self.angle__rad
    damping_torque = -self._damping * self.velocity__rad_per_s

    angular_acceleration = (
        torque__Nm + spring_torque + damping_torque
    ) / self._inertia

    # Euler integration (could use RK4 for higher accuracy)
    self.velocity__rad_per_s += angular_acceleration * dt__s
    self.angle__rad += self.velocity__rad_per_s * dt__s

    return self.angle__deg, np.degrees(self.velocity__rad_per_s)

reset

reset() -> None

Reset joint to initial conditions.

Source code in myogen/simulator/neuron/joint_dynamics.py
def reset(self) -> None:
    """Reset joint to initial conditions."""
    self.angle__rad = np.radians(self.initial_angle__deg)
    self.velocity__rad_per_s = np.radians(self.initial_velocity__deg_per_s)

get_state

get_state() -> dict

Get current joint state.

Returns:

Type Description
dict

Dictionary containing current angle, velocity, and parameters.

Source code in myogen/simulator/neuron/joint_dynamics.py
def get_state(self) -> dict:
    """
    Get current joint state.

    Returns
    -------
    dict
        Dictionary containing current angle, velocity, and parameters.
    """
    return {
        "angle__deg": self.angle__deg,
        "angle__rad": self.angle__rad,
        "velocity__deg_per_s": np.degrees(self.velocity__rad_per_s),
        "velocity__rad_per_s": self.velocity__rad_per_s,
        "inertia__kg_m2": self._inertia,
        "damping__Nm_s_per_rad": self._damping,
        "stiffness__Nm_per_rad": self._stiffness,
    }

Geometry & biomechanics

MuscleGeometry dataclass

MuscleGeometry(origin_coords: tuple[float, float], insertion_coords: tuple[float, float], optimal_length__cm: float, max_length_change__cm: float)

Geometric parameters for a muscle crossing a joint.

JointGeometry dataclass

JointGeometry(center_coords: tuple[float, float], radius__cm: float, range_of_motion__degrees: tuple[float, float])

Geometric parameters for a joint.

JointBiomechanics

JointBiomechanics(joint_type: Literal['hinge', 'ball_socket'], joint_geometry: JointGeometry, muscle_geometries: list[MuscleGeometry], moment_arm_data__cm: Optional[MOMENT_ARM__MATRIX] = None, use_simplified_model: bool = True)

Joint biomechanics with muscle moment arms and length calculations.

This class computes the relationship between joint angles and muscle lengths, moment arms, and resulting joint torques. It supports both simple geometric models and more complex biomechanical relationships.

Parameters:

Name Type Description Default
joint_type (hinge, ball_socket)

Type of joint for biomechanical calculations. "hinge" = single degree of freedom (e.g., elbow, knee) "ball_socket" = multi-degree of freedom (e.g., shoulder, hip)

"hinge"
joint_geometry JointGeometry

Geometric parameters of the joint.

required
muscle_geometries list[MuscleGeometry]

List of muscle geometry parameters for muscles crossing this joint.

required
moment_arm_data__cm MOMENT_ARM__MATRIX

Pre-computed moment arm data as function of joint angle. If None, computed from geometry. Shape: (n_angles, n_muscles).

None
use_simplified_model bool

Whether to use simplified geometric model or detailed biomechanics.

True

Attributes:

Name Type Description
n_muscles int

Number of muscles crossing this joint.

joint_angles__deg ndarray

Array of joint angles for moment arm calculations.

moment_arms__cm ndarray

Moment arms for each muscle at each joint angle.

Source code in myogen/simulator/core/force/biomechanics.py
def __init__(
    self,
    joint_type: Literal["hinge", "ball_socket"],
    joint_geometry: JointGeometry,
    muscle_geometries: list[MuscleGeometry],
    moment_arm_data__cm: Optional[MOMENT_ARM__MATRIX] = None,
    use_simplified_model: bool = True,
) -> None:
    # Validate inputs
    if joint_type not in ["hinge", "ball_socket"]:
        raise ValueError(f"joint_type must be 'hinge' or 'ball_socket', got '{joint_type}'")

    if len(muscle_geometries) == 0:
        raise ValueError("At least one muscle geometry must be provided")

    if joint_geometry.radius__cm <= 0:
        raise ValueError(
            f"joint_geometry.radius__cm must be positive, got {joint_geometry.radius__cm}"
        )

    min_angle, max_angle = joint_geometry.range_of_motion__degrees
    if min_angle >= max_angle:
        raise ValueError(
            f"Joint range of motion invalid: min={min_angle}, max={max_angle}. "
            "Maximum angle must be greater than minimum angle."
        )

    # Store parameters (immutable public access)
    self.joint_type = joint_type
    self.joint_geometry = joint_geometry
    self.muscle_geometries = muscle_geometries
    self.use_simplified_model = use_simplified_model

    # Private copies for internal use
    self._joint_type = joint_type
    self._joint_geometry = joint_geometry
    self._muscle_geometries = muscle_geometries
    self._use_simplified = use_simplified_model

    # Derived properties
    self.n_muscles = len(muscle_geometries)

    # Set up angle array for calculations
    min_angle, max_angle = joint_geometry.range_of_motion__degrees
    self.joint_angles__deg = np.linspace(min_angle, max_angle, 181)  # 1 degree resolution

    # Initialize or validate moment arm data
    if moment_arm_data__cm is not None:
        if moment_arm_data__cm.shape != (
            len(self.joint_angles__deg),
            self.n_muscles,
        ):
            raise ValueError(
                f"moment_arm_data__cm shape {moment_arm_data__cm.shape} does not match "
                f"expected ({len(self.joint_angles__deg)}, {self.n_muscles})"
            )
        self._moment_arms__cm = moment_arm_data__cm
    else:
        self._moment_arms__cm = self._compute_moment_arms()

    # Public access to moment arms
    self.moment_arms__cm = self._moment_arms__cm

compute_muscle_length

compute_muscle_length(joint_angle__deg: float, muscle_index: int) -> float

Compute muscle length from joint angle.

Parameters:

Name Type Description Default
joint_angle__deg float

Joint angle in degrees.

required
muscle_index int

Index of the muscle (0 to n_muscles-1).

required

Returns:

Type Description
float

Muscle length in cm.

Raises:

Type Description
ValueError

If muscle_index is out of range.

Source code in myogen/simulator/core/force/biomechanics.py
def compute_muscle_length(self, joint_angle__deg: float, muscle_index: int) -> float:
    """
    Compute muscle length from joint angle.

    Parameters
    ----------
    joint_angle__deg : float
        Joint angle in degrees.
    muscle_index : int
        Index of the muscle (0 to n_muscles-1).

    Returns
    -------
    float
        Muscle length in cm.

    Raises
    ------
    ValueError
        If muscle_index is out of range.
    """
    if not 0 <= muscle_index < self.n_muscles:
        raise ValueError(
            f"muscle_index must be between 0 and {self.n_muscles - 1}, got {muscle_index}"
        )

    muscle_geom = self._muscle_geometries[muscle_index]

    if self._use_simplified:
        # Simplified model: length varies with joint angle
        angle_rad = np.radians(joint_angle__deg)

        # Length change is proportional to joint rotation and moment arm
        moment_arm = self.get_moment_arm(joint_angle__deg, muscle_index)

        # Reference angle (middle of range of motion)
        min_angle, max_angle = self._joint_geometry.range_of_motion__degrees
        ref_angle = (min_angle + max_angle) / 2

        angle_change = np.radians(joint_angle__deg - ref_angle)
        length_change = moment_arm * angle_change

        muscle_length = muscle_geom.optimal_length__cm + length_change

        # Constrain to reasonable range
        min_length = muscle_geom.optimal_length__cm - muscle_geom.max_length_change__cm
        max_length = muscle_geom.optimal_length__cm + muscle_geom.max_length_change__cm

        return np.clip(muscle_length, min_length, max_length)
    else:
        # Detailed geometric calculation
        return self._compute_geometric_muscle_length(joint_angle__deg, muscle_index)

get_moment_arm

get_moment_arm(joint_angle__deg: float, muscle_index: int) -> float

Get moment arm for a muscle at a specific joint angle.

Parameters:

Name Type Description Default
joint_angle__deg float

Joint angle in degrees.

required
muscle_index int

Index of the muscle.

required

Returns:

Type Description
float

Moment arm in cm.

Source code in myogen/simulator/core/force/biomechanics.py
def get_moment_arm(self, joint_angle__deg: float, muscle_index: int) -> float:
    """
    Get moment arm for a muscle at a specific joint angle.

    Parameters
    ----------
    joint_angle__deg : float
        Joint angle in degrees.
    muscle_index : int
        Index of the muscle.

    Returns
    -------
    float
        Moment arm in cm.
    """
    if not 0 <= muscle_index < self.n_muscles:
        raise ValueError(
            f"muscle_index must be between 0 and {self.n_muscles - 1}, got {muscle_index}"
        )

    # Interpolate from pre-computed moment arms
    return np.interp(
        joint_angle__deg,
        self.joint_angles__deg,
        self._moment_arms__cm[:, muscle_index],
    )

compute_joint_torque

compute_joint_torque(muscle_forces__N: Union[float, ndarray], joint_angle__deg: float) -> float

Compute net joint torque from muscle forces.

Parameters:

Name Type Description Default
muscle_forces__N float or ndarray

Forces from each muscle in Newtons. If float, assumes single muscle. If array, must have length equal to n_muscles.

required
joint_angle__deg float

Current joint angle in degrees.

required

Returns:

Type Description
float

Net joint torque in N⋅cm.

Notes

Positive torque indicates rotation in the positive joint angle direction.

Source code in myogen/simulator/core/force/biomechanics.py
def compute_joint_torque(
    self, muscle_forces__N: Union[float, np.ndarray], joint_angle__deg: float
) -> float:
    """
    Compute net joint torque from muscle forces.

    Parameters
    ----------
    muscle_forces__N : float or np.ndarray
        Forces from each muscle in Newtons. If float, assumes single muscle.
        If array, must have length equal to n_muscles.
    joint_angle__deg : float
        Current joint angle in degrees.

    Returns
    -------
    float
        Net joint torque in N⋅cm.

    Notes
    -----
    Positive torque indicates rotation in the positive joint angle direction.
    """
    forces = np.asarray(muscle_forces__N)

    if forces.ndim == 0:  # Single muscle
        if self.n_muscles != 1:
            raise ValueError(f"Single force provided but {self.n_muscles} muscles defined")
        forces = np.array([forces])
    elif len(forces) != self.n_muscles:
        raise ValueError(
            f"Force array length ({len(forces)}) must match number of muscles ({self.n_muscles})"
        )

    # Get moment arms for current joint angle
    moment_arms = np.array(
        [self.get_moment_arm(joint_angle__deg, i) for i in range(self.n_muscles)]
    )

    # Compute torques (force × moment arm)
    muscle_torques = forces * moment_arms

    # Sum torques (considering muscle action directions would require additional info)
    # For now, assume all muscles act in same direction
    return np.sum(muscle_torques)

get_muscle_length_trajectory

get_muscle_length_trajectory(joint_angle_trajectory__degrees: JOINT_ANGLE__ARRAY, muscle_index: int) -> ndarray

Compute muscle length trajectory from joint angle trajectory.

Parameters:

Name Type Description Default
joint_angle_trajectory__degrees JOINT_ANGLE__ARRAY

Array of joint angles over time.

required
muscle_index int

Index of the muscle.

required

Returns:

Type Description
ndarray

Muscle length trajectory in cm.

Source code in myogen/simulator/core/force/biomechanics.py
def get_muscle_length_trajectory(
    self, joint_angle_trajectory__degrees: JOINT_ANGLE__ARRAY, muscle_index: int
) -> np.ndarray:
    """
    Compute muscle length trajectory from joint angle trajectory.

    Parameters
    ----------
    joint_angle_trajectory__degrees : JOINT_ANGLE__ARRAY
        Array of joint angles over time.
    muscle_index : int
        Index of the muscle.

    Returns
    -------
    np.ndarray
        Muscle length trajectory in cm.
    """
    return np.array(
        [
            self.compute_muscle_length(angle, muscle_index)
            for angle in joint_angle_trajectory__degrees
        ]
    )

get_biomechanical_summary

get_biomechanical_summary() -> dict

Get summary of biomechanical parameters.

Returns:

Type Description
dict

Dictionary containing key biomechanical parameters.

Source code in myogen/simulator/core/force/biomechanics.py
def get_biomechanical_summary(self) -> dict:
    """
    Get summary of biomechanical parameters.

    Returns
    -------
    dict
        Dictionary containing key biomechanical parameters.
    """
    min_angle, max_angle = self._joint_geometry.range_of_motion__degrees

    return {
        "joint_type": self._joint_type,
        "n_muscles": self.n_muscles,
        "joint_radius__cm": self._joint_geometry.radius__cm,
        "joint_range_of_motion__degrees": (min_angle, max_angle),
        "use_simplified_model": self._use_simplified,
        "moment_arm_range__cm": {
            f"muscle_{i}": (
                np.min(self._moment_arms__cm[:, i]),
                np.max(self._moment_arms__cm[:, i]),
            )
            for i in range(self.n_muscles)
        },
        "optimal_muscle_lengths__cm": [
            geom.optimal_length__cm for geom in self._muscle_geometries
        ],
    }

Grid / Neo utilities

create_grid_signal

create_grid_signal(signal: ndarray, grid_shape: tuple[int, int], sampling_rate: Quantity, units: Quantity | str = mV, t_start: Quantity = 0 * s, electrode_positions: list[tuple[float, float]] | None = None, ied: float | None = None, **kwargs) -> AnalogSignal

Create an AnalogSignal with grid metadata annotations for electrode arrays.

This function creates a standard Neo AnalogSignal with grid structure stored in annotations, making it NWB-compatible while preserving spatial information.

Parameters:

Name Type Description Default
signal ndarray

Signal data with shape (time, rows, cols) or (time, n_electrodes). If 3D, will be flattened to (time, n_electrodes) for storage.

required
grid_shape tuple[int, int]

Shape of the electrode grid as (rows, cols).

required
sampling_rate Quantity

Sampling rate of the signal.

required
units Quantity or str

Units of the signal.

pq.mV
t_start Quantity

Start time of the signal.

0*pq.s
electrode_positions list[tuple[float, float]]

Physical (x, y) positions of each electrode in mm. If None, positions are computed from grid_shape and ied.

None
ied float

Inter-electrode distance in mm. Used to compute electrode_positions if not provided directly.

None
**kwargs

Additional arguments passed to AnalogSignal constructor or stored as annotations.

{}

Returns:

Type Description
AnalogSignal

Neo AnalogSignal with grid metadata in annotations: - 'grid_shape': (rows, cols) tuple - 'electrode_positions': list of (x, y) tuples in mm - 'ied': inter-electrode distance in mm (if provided)

Examples:

>>> import numpy as np
>>> import quantities as pq
>>> from myogen.utils.neo import create_grid_signal, signal_to_grid
>>>
>>> # Create grid data (time, rows, cols)
>>> data = np.random.rand(1000, 8, 8)
>>> signal = create_grid_signal(
...     data,
...     grid_shape=(8, 8),
...     sampling_rate=2048 * pq.Hz,
...     ied=8.0,  # 8mm inter-electrode distance
... )
>>>
>>> # Access as grid
>>> grid = signal_to_grid(signal)  # shape: (1000, 8, 8)
>>>
>>> # Access single electrode
>>> row, col = 2, 3
>>> electrode_idx = row * 8 + col
>>> single = signal[:, electrode_idx]
Source code in myogen/utils/neo.py
@beartowertype
def create_grid_signal(
    signal: np.ndarray,
    grid_shape: tuple[int, int],
    sampling_rate: pq.Quantity,
    units: pq.Quantity | str = pq.mV,
    t_start: pq.Quantity = 0 * pq.s,
    electrode_positions: list[tuple[float, float]] | None = None,
    ied: float | None = None,
    **kwargs,
) -> AnalogSignal:
    """
    Create an AnalogSignal with grid metadata annotations for electrode arrays.

    This function creates a standard Neo AnalogSignal with grid structure stored
    in annotations, making it NWB-compatible while preserving spatial information.

    Parameters
    ----------
    signal : np.ndarray
        Signal data with shape (time, rows, cols) or (time, n_electrodes).
        If 3D, will be flattened to (time, n_electrodes) for storage.
    grid_shape : tuple[int, int]
        Shape of the electrode grid as (rows, cols).
    sampling_rate : pq.Quantity
        Sampling rate of the signal.
    units : pq.Quantity or str, default=pq.mV
        Units of the signal.
    t_start : pq.Quantity, default=0*pq.s
        Start time of the signal.
    electrode_positions : list[tuple[float, float]], optional
        Physical (x, y) positions of each electrode in mm. If None,
        positions are computed from grid_shape and ied.
    ied : float, optional
        Inter-electrode distance in mm. Used to compute electrode_positions
        if not provided directly.
    **kwargs
        Additional arguments passed to AnalogSignal constructor or stored
        as annotations.

    Returns
    -------
    AnalogSignal
        Neo AnalogSignal with grid metadata in annotations:
        - 'grid_shape': (rows, cols) tuple
        - 'electrode_positions': list of (x, y) tuples in mm
        - 'ied': inter-electrode distance in mm (if provided)

    Examples
    --------
    >>> import numpy as np
    >>> import quantities as pq
    >>> from myogen.utils.neo import create_grid_signal, signal_to_grid
    >>>
    >>> # Create grid data (time, rows, cols)
    >>> data = np.random.rand(1000, 8, 8)
    >>> signal = create_grid_signal(
    ...     data,
    ...     grid_shape=(8, 8),
    ...     sampling_rate=2048 * pq.Hz,
    ...     ied=8.0,  # 8mm inter-electrode distance
    ... )
    >>>
    >>> # Access as grid
    >>> grid = signal_to_grid(signal)  # shape: (1000, 8, 8)
    >>>
    >>> # Access single electrode
    >>> row, col = 2, 3
    >>> electrode_idx = row * 8 + col
    >>> single = signal[:, electrode_idx]
    """
    signal_array = np.asarray(signal)

    # Handle units from input data
    if hasattr(signal, "units") and units == pq.mV:
        units = signal.units

    # Flatten 3D to 2D if necessary
    if signal_array.ndim == 3:
        time_points, rows, cols = signal_array.shape
        if (rows, cols) != grid_shape:
            raise ValueError(
                f"Signal shape {(rows, cols)} does not match grid_shape {grid_shape}"
            )
        signal_2d = signal_array.reshape(time_points, rows * cols)
    elif signal_array.ndim == 2:
        signal_2d = signal_array
        rows, cols = grid_shape
        expected_channels = rows * cols
        if signal_2d.shape[1] != expected_channels:
            raise ValueError(
                f"Signal has {signal_2d.shape[1]} channels but grid_shape {grid_shape} "
                f"expects {expected_channels} channels"
            )
    else:
        raise ValueError(
            f"Signal must be 2D or 3D array, got {signal_array.ndim}D with shape {signal_array.shape}"
        )

    # Compute electrode positions if not provided
    if electrode_positions is None and ied is not None:
        rows, cols = grid_shape
        electrode_positions = []
        for r in range(rows):
            for c in range(cols):
                x = c * ied
                y = r * ied
                electrode_positions.append((x, y))

    # Separate Neo kwargs from annotation kwargs
    neo_kwargs = {}
    annotation_kwargs = {}
    neo_params = {"name", "description", "file_origin", "array_annotations", "copy"}

    for key, value in kwargs.items():
        if key in neo_params:
            neo_kwargs[key] = value
        else:
            annotation_kwargs[key] = value

    # Create the AnalogSignal
    analog_signal = AnalogSignal(
        signal_2d * units if not hasattr(signal_2d, "units") else signal_2d,
        sampling_rate=sampling_rate,
        t_start=t_start,
        **neo_kwargs,
    )

    # Add grid annotations
    analog_signal.annotate(
        grid_shape=grid_shape,
        electrode_positions=electrode_positions,
        ied=ied,
        **annotation_kwargs,
    )

    return analog_signal

signal_to_grid

signal_to_grid(signal: AnalogSignal, time_slice: slice | None = None) -> ndarray

Convert a grid-annotated AnalogSignal back to 3D grid format.

Parameters:

Name Type Description Default
signal AnalogSignal

AnalogSignal with 'grid_shape' annotation.

required
time_slice slice

Time slice to extract. If None, returns all time points.

None

Returns:

Type Description
ndarray

Data in grid format with shape (time, rows, cols).

Raises:

Type Description
ValueError

If the signal doesn't have grid_shape annotation.

Examples:

>>> grid = signal_to_grid(signal)
>>> grid.shape  # (time, rows, cols)
(1000, 8, 8)
>>>
>>> # Get single time point
>>> frame = signal_to_grid(signal, time_slice=slice(100, 101))
>>> frame.shape
(1, 8, 8)
Source code in myogen/utils/neo.py
@beartowertype
def signal_to_grid(signal: AnalogSignal, time_slice: slice | None = None) -> np.ndarray:
    """
    Convert a grid-annotated AnalogSignal back to 3D grid format.

    Parameters
    ----------
    signal : AnalogSignal
        AnalogSignal with 'grid_shape' annotation.
    time_slice : slice, optional
        Time slice to extract. If None, returns all time points.

    Returns
    -------
    np.ndarray
        Data in grid format with shape (time, rows, cols).

    Raises
    ------
    ValueError
        If the signal doesn't have grid_shape annotation.

    Examples
    --------
    >>> grid = signal_to_grid(signal)
    >>> grid.shape  # (time, rows, cols)
    (1000, 8, 8)
    >>>
    >>> # Get single time point
    >>> frame = signal_to_grid(signal, time_slice=slice(100, 101))
    >>> frame.shape
    (1, 8, 8)
    """
    if "grid_shape" not in signal.annotations:
        raise ValueError(
            "Signal does not have 'grid_shape' annotation. "
            "Use create_grid_signal() to create grid-annotated signals."
        )

    grid_shape = signal.annotations["grid_shape"]
    rows, cols = grid_shape

    data = signal.magnitude
    if time_slice is not None:
        data = data[time_slice]

    if data.ndim == 1:
        # Single time point
        return data.reshape(1, rows, cols)
    else:
        return data.reshape(-1, rows, cols)

get_electrode

get_electrode(signal: AnalogSignal, row: int, col: int) -> AnalogSignal

Extract a single electrode's signal from a grid-annotated AnalogSignal.

Parameters:

Name Type Description Default
signal AnalogSignal

AnalogSignal with 'grid_shape' annotation.

required
row int

Row index of the electrode.

required
col int

Column index of the electrode.

required

Returns:

Type Description
AnalogSignal

Single-channel AnalogSignal for the specified electrode.

Examples:

>>> electrode_signal = get_electrode(signal, row=2, col=3)
>>> electrode_signal.shape
(1000, 1)
Source code in myogen/utils/neo.py
@beartowertype
def get_electrode(
    signal: AnalogSignal,
    row: int,
    col: int,
) -> AnalogSignal:
    """
    Extract a single electrode's signal from a grid-annotated AnalogSignal.

    Parameters
    ----------
    signal : AnalogSignal
        AnalogSignal with 'grid_shape' annotation.
    row : int
        Row index of the electrode.
    col : int
        Column index of the electrode.

    Returns
    -------
    AnalogSignal
        Single-channel AnalogSignal for the specified electrode.

    Examples
    --------
    >>> electrode_signal = get_electrode(signal, row=2, col=3)
    >>> electrode_signal.shape
    (1000, 1)
    """
    if "grid_shape" not in signal.annotations:
        raise ValueError("Signal does not have 'grid_shape' annotation.")

    rows, cols = signal.annotations["grid_shape"]
    if row < 0 or row >= rows:
        raise ValueError(f"Row {row} out of bounds for grid with {rows} rows")
    if col < 0 or col >= cols:
        raise ValueError(f"Column {col} out of bounds for grid with {cols} columns")

    channel_idx = row * cols + col
    return signal[:, channel_idx]

get_row

get_row(signal: AnalogSignal, row: int) -> AnalogSignal

Extract all electrodes from a specific row.

Parameters:

Name Type Description Default
signal AnalogSignal

AnalogSignal with 'grid_shape' annotation.

required
row int

Row index to extract.

required

Returns:

Type Description
AnalogSignal

Multi-channel AnalogSignal for all electrodes in the row.

Examples:

>>> row_signal = get_row(signal, row=2)
>>> row_signal.shape  # (time, n_cols)
(1000, 8)
Source code in myogen/utils/neo.py
@beartowertype
def get_row(
    signal: AnalogSignal,
    row: int,
) -> AnalogSignal:
    """
    Extract all electrodes from a specific row.

    Parameters
    ----------
    signal : AnalogSignal
        AnalogSignal with 'grid_shape' annotation.
    row : int
        Row index to extract.

    Returns
    -------
    AnalogSignal
        Multi-channel AnalogSignal for all electrodes in the row.

    Examples
    --------
    >>> row_signal = get_row(signal, row=2)
    >>> row_signal.shape  # (time, n_cols)
    (1000, 8)
    """
    if "grid_shape" not in signal.annotations:
        raise ValueError("Signal does not have 'grid_shape' annotation.")

    rows, cols = signal.annotations["grid_shape"]
    if row < 0 or row >= rows:
        raise ValueError(f"Row {row} out of bounds for grid with {rows} rows")

    start_idx = row * cols
    end_idx = start_idx + cols
    return signal[:, start_idx:end_idx]

get_column

get_column(signal: AnalogSignal, col: int) -> AnalogSignal

Extract all electrodes from a specific column.

Parameters:

Name Type Description Default
signal AnalogSignal

AnalogSignal with 'grid_shape' annotation.

required
col int

Column index to extract.

required

Returns:

Type Description
AnalogSignal

Multi-channel AnalogSignal for all electrodes in the column.

Examples:

>>> col_signal = get_column(signal, col=3)
>>> col_signal.shape  # (time, n_rows)
(1000, 8)
Source code in myogen/utils/neo.py
@beartowertype
def get_column(
    signal: AnalogSignal,
    col: int,
) -> AnalogSignal:
    """
    Extract all electrodes from a specific column.

    Parameters
    ----------
    signal : AnalogSignal
        AnalogSignal with 'grid_shape' annotation.
    col : int
        Column index to extract.

    Returns
    -------
    AnalogSignal
        Multi-channel AnalogSignal for all electrodes in the column.

    Examples
    --------
    >>> col_signal = get_column(signal, col=3)
    >>> col_signal.shape  # (time, n_rows)
    (1000, 8)
    """
    if "grid_shape" not in signal.annotations:
        raise ValueError("Signal does not have 'grid_shape' annotation.")

    rows, cols = signal.annotations["grid_shape"]
    if col < 0 or col >= cols:
        raise ValueError(f"Column {col} out of bounds for grid with {cols} columns")

    channel_indices = [r * cols + col for r in range(rows)]
    return signal[:, channel_indices]

The deprecated GridAnalogSignal compatibility class is intentionally excluded.