Source code for pymuscle.potvin_motor_neuron_pool
import numpy as np
from numpy import ndarray
from .model import Model
[docs]class PotvinMotorNeuronPool(Model):
"""
Encapsulates the motor neuron portion of the motor unit model.
The name of each parameter as it appears in Potvin, 2017 is in parentheses.
If a parameter does not appear in the paper but does appear in the
accompanying Matlab code, the variable name from the Matlab code is used in
the parentheses.
:param motor_unit_count: Number of motor units in the muscle (n)
:param max_recruitment_threshold:
Max excitation required by a motor unit within the pool before
firing (RR)
:param firing_gain:
The slope of firing rate by excitation above threshold (g)
:param min_firing_rate:
The minimum firing rate for a motor neuron above threshold (minR)
:param max_firing_rate_first_unit:
Max firing rate for the first motor unit (maxR(1))
:param max_firing_rate_last_unit:
Max firing rate for the last motor unit (maxR(last))
:param pre_calc_firing_rates:
Whether to build a dict mapping excitation levels to firing rates for
each motor neuron. This can speed up simulation at the cost of
additional memory.
:param pre_calc_resolution:
Step size for excitation levels to pre-calculate (res)
:param pre_calc_max: Highest excitation value to pre-calculate
Usage::
from pymuscle import PotvinMotorNeuronPool
motor_unit_count = 60
pool = PotvinMotorNeuronPool(motor_unit_count)
excitation = np.full(motor_unit_count, 10.0)
firing_rates = pool.step(excitation)
"""
def __init__(
self,
motor_unit_count: int = 120,
max_recruitment_threshold: int = 50,
firing_gain: float = 1.0,
min_firing_rate: int = 8,
max_firing_rate_first_unit: int = 35,
max_firing_rate_last_unit: int = 25,
pre_calc_firing_rates: bool = True,
pre_calc_resolution: float = 0.1,
pre_calc_max: float = 70.0
):
self._recruitment_thresholds = self._calc_recruitment_thresholds(
motor_unit_count,
max_recruitment_threshold
)
self._peak_firing_rates = self._calc_peak_firing_rates(
max_firing_rate_first_unit,
max_firing_rate_last_unit,
max_recruitment_threshold,
self._recruitment_thresholds
)
# Assign additional non-public attributes
self._max_recruitment_threshold = max_recruitment_threshold
self._firing_gain = firing_gain
self._min_firing_rate = min_firing_rate
# Assign public attributes
self.motor_unit_count = motor_unit_count
# Pre-calculate firing rates for all motor neurons across a range of
# possible excitation levels.
self._firing_rates_by_excitation = None
if pre_calc_firing_rates:
self._build_firing_rates_cache(pre_calc_max, pre_calc_resolution)
def _build_firing_rates_cache(
self,
pre_calc_max: float,
pre_calc_resolution: float
) -> None:
"""
Pre-calculate and store firing rates for all motor neurons across a
range of possible excitation levels.
"""
self._firing_rates_by_excitation = {}
# TODO: This is a hack. Maybe memoize vs pre-calculate?
# Is there a good LRU cache decorator I can drop in?
# Maybe https://docs.python.org/3/library/functools.html#functools.lru_cache
resolution_places = len(str(pre_calc_resolution).split(".")[1])
excitations = np.zeros(self.motor_unit_count)
excitation_values = np.arange(
0.0,
pre_calc_max + pre_calc_resolution,
pre_calc_resolution
)
for i in excitation_values:
i = round(i, resolution_places)
excitations += pre_calc_resolution
self._firing_rates_by_excitation[i] = \
self._inner_calc_firing_rates(
excitations,
self._recruitment_thresholds,
self._firing_gain,
self._min_firing_rate,
self._peak_firing_rates
)
def _calc_firing_rates(self, excitations: ndarray) -> ndarray:
"""
Calculates firing rates on a per motor neuron basis for the given
array of excitations.
"""
assert (len(excitations) == len(self._recruitment_thresholds))
if self._firing_rates_by_excitation:
excitation = excitations[0] # TODO - Support variations
firing_rates = self._firing_rates_by_excitation[excitation]
else:
firing_rates = self._inner_calc_firing_rates(
excitations,
self._recruitment_thresholds,
self._firing_gain,
self._min_firing_rate,
self._peak_firing_rates
)
return firing_rates
@staticmethod
def _calc_peak_firing_rates(
max_firing_rate_first_unit: int,
max_firing_rate_last_unit: int,
max_recruitment_threshold: int,
recruitment_thresholds: ndarray,
) -> ndarray:
"""
Calculate peak firing rates for each motor neuron
"""
firing_rate_range = (
max_firing_rate_first_unit - max_firing_rate_last_unit
)
first_recruitment_thresh = recruitment_thresholds[0]
recruitment_thresh_range = (
max_recruitment_threshold - first_recruitment_thresh
)
temp_thresholds = (
recruitment_thresholds - first_recruitment_thresh
)
temp_thresholds /= recruitment_thresh_range
return (
max_firing_rate_first_unit - (firing_rate_range * temp_thresholds)
)
@staticmethod
def _calc_recruitment_thresholds(
motor_unit_count: int,
max_recruitment_threshold: int
) -> ndarray:
"""
Calculate recruitment thresholds for each motor neuron
"""
motor_unit_indices = np.arange(1, motor_unit_count + 1)
r_log = np.log(max_recruitment_threshold)
r_exponent = (r_log * (motor_unit_indices)) / (motor_unit_count)
return np.exp(r_exponent)
@staticmethod
def _inner_calc_firing_rates(
excitations: ndarray,
thresholds: ndarray,
gain: float,
min_firing_rate: int,
peak_firing_rates: ndarray
) -> ndarray:
firing_rates = excitations - thresholds
firing_rates += min_firing_rate
below_thresh_indices = firing_rates < min_firing_rate
firing_rates[below_thresh_indices] = 0
firing_rates *= gain
# Check for max values
above_peak_indices = firing_rates > peak_firing_rates
firing_rates[above_peak_indices] = peak_firing_rates[above_peak_indices]
return firing_rates
[docs] def step(self, motor_pool_input: ndarray) -> ndarray:
"""
Advance the motor neuron pool simulation one step.
Returns firing rates on a per motor neuron basis for the given
array of excitations.
"""
return self._calc_firing_rates(motor_pool_input)