# NeuroMotion: Open-source Simulator with Neuromechanical and Deep Network Models to Generate Surface EMG signals during Voluntary Movement

Shihan Ma, Irene Mendez Guerra, Arnault Hubert Caillet, Jiamin Zhao, Alexander Kenneth Clarke, Kostiantyn Maksymenko, Samuel Deslauriers-Gauthier, Xinjun Sheng, Xiangyang Zhu, Dario Farina

<figure><img src="image_1.png" /></figure>

## To cite this version:

Shihan Ma, Irene Mendez Guerra, Arnault Hubert Caillet, Jiamin Zhao, Alexander Kenneth Clarke, et al.. NeuroMotion: Open-source Simulator with Neuromechanical and Deep Network Models to Generate Surface EMG signals during Voluntary Movement. PLoS Computational Biology, 2024, 10.1371/journal.pcbi.1012257. [https://doi.org/10.1371/journal.pcbi.1012257](https://doi.org/10.1371/journal.pcbi.1012257)

HAL Id: hal-04390169

https://inria.hal.science/hal-04390169v2

Submitted on 20 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
METHODS

# NeuroMotion: Open-source platform with neuromechanical and deep network modules to generate surface EMG signals during voluntary movement

Shihan Ma¹,², Irene Mendez Guerra¹, Arnault Hubert Caillet¹, Jiamin Zhao², Alexander Kenneth Clarke¹, Kostiantyn Maksymenko³, Samuel Deslauriers-Gauthier³,⁴, Xinjun Sheng²,⁵, Xiangyang Zhu²,⁵, Dario Farina¹*

1 Department of Bioengineering, Imperial College London, London, United Kingdom, 2 State Key Laboratory of Mechanical System and Vibration, Shanghai Jiao Tong University, Shanghai, China, 3 Neurodec, Sophia Antipolis, France, 4 Inria Centre at Université Côte d'Azur, Nice, France, 5 Meta Robotics Institute, Shanghai Jiao Tong University, Shanghai, China

* d.farina@imperial.ac.uk

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

OPEN ACCESS

**Citation:** Ma S, Mendez Guerra I, Caillet AH, Zhao J, Clarke AK, Maksymenko K, et al. (2024) NeuroMotion: Open-source platform with neuromechanical and deep network modules to generate surface EMG signals during voluntary movement. PLoS Comput Biol 20(7): e1012257. https://doi.org/10.1371/journal.pcbi.1012257

**Editor:** Adrian M Haith, Johns Hopkins University, UNITED STATES

**Received:** October 9, 2023

**Accepted:** June 15, 2024

**Published:** July 3, 2024

**Peer Review History:** PLOS recognizes the benefits of transparency in the peer review process; therefore, we enable the publication of all of the content of peer review and author responses alongside final, published articles. The editorial history of this article is available here: https://doi.org/10.1371/journal.pcbi.1012257

**Copyright:** © 2024 Ma et al. This is an open access article distributed under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All programming code are publicly available on GitHub, at https://github.com/shihan-ma/NeuroMotion under the

## Abstract

Neuromechanical studies investigate how the nervous system interacts with the musculo-skeletal (MSK) system to generate volitional movements. Such studies have been supported by simulation models that provide insights into variables that cannot be measured experimentally and allow a large number of conditions to be tested before the experimental analysis. However, current simulation models of electromyography (EMG), a core physiological signal in neuromechanical analyses, remain either limited in accuracy and conditions or are computationally heavy to apply. Here, we provide a computational platform to enable future work to overcome these limitations by presenting NeuroMotion, an open-source simulator that can modularly test a variety of approaches to the full-spectrum synthesis of EMG signals during voluntary movements. We demonstrate NeuroMotion using three sample modules. The first module is an upper-limb MSK model with OpenSim API to estimate the muscle fibre lengths and muscle activations during movements. The second module is Bio-Mime, a deep neural network-based EMG generator that receives nonstationary physiological parameter inputs, like the afore-estimated muscle fibre lengths, and efficiently outputs motor unit action potentials (MUAPs). The third module is a motor unit pool model that transforms the muscle activations into discharge timings of motor units. The discharge timings are convolved with the output of BioMime to simulate EMG signals during the movement. We first show how MUAP waveforms change during different levels of physiological parameter variations and different movements. We then show that the synthetic EMG signals during two-degree-of-freedom hand and wrist movements can be used to augment experimental data for regressing joint angles. Ridge regressors trained on the synthetic dataset were directly used to predict joint angles from experimental data. In this way, NeuroMotion was able to generate full-spectrum EMG for the first use-case of human forearm electrophysiology during voluntary hand, wrist, and forearm movements. All intermediate variables are available, which allows the user to study cause-effect relationships in the
GNU General Public License v3.0. All experimental
data is available at https://doi.org/10.6084/m9.
figshare.25981795.v1 under the CC BY 4.0
License.

**Funding:** This work was supported by the European Research Council (ERC) under the Synergy Grant Natural BionicS (810346 to DF), the Engineering & Physical Sciences Research Council (EPSRC) Transformative Healthcare for 2050 project Non-Invasive Single Neuron Electrical Monitoring (NISNEM Technology, grant EP/ T020970/1 to DF), and by the National Natural Science Foundation of China (52227808 to XZ). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** • I have read the journal's policy and the authors of this manuscript have the following competing interests: Kostiantyn Maksymenko and Samuel Deslauriers-Gauthier are the co-founders of the company Neurodec; all other authors have no competing interests.

complex neuromechanical system, fast iterate algorithms before collecting experimental
data, and validate algorithms that estimate non-measurable parameters in experiments. We
expect this modular platform will enable validation of generative EMG models, complement
experimental approaches and empower neuromechanical research.

Author summary

Neuromechanical studies investigate how the nervous system and musculoskeletal system interact to generate movements. Such studies heavily rely on simulation models, which provide non-measurable variables to complement the experimental analyses. However, the simulation models of surface electromyography (EMG), the core physiological signal widely used in neuromechanical analyses, are limited to static conditions. We bridged this gap by proposing NeuroMotion, the first full-spectrum EMG simulator that can be used to generate EMG signals during voluntary movements. NeuroMotion integrates a muscu- loskeletal model, a neural network-based EMG generator, and advanced motoneuron models. With representative applications of this simulator, we show that it can be used to investigate the variabilities of EMG signals during voluntary movement. We also demon- strate that the synthetic signals generated by NeuroMotion can be used to augment experi- mental data for regressing joint angles. We expect the functionality provided by NeuroMotion, which is provided open-source, will stimulate progress in neuromechanics.

Introduction

Human neuromechanics is the discipline that combines neuroscience and biomechanics to reach a fundamental understanding of the interactions between the nervous, muscular, and skeletal systems during human movements [1]. It allows us to uncover the functions and mechanisms of the nervous system under the production of movements [2, 3] by studying the movements from the perspective of their neural control [4]. Neuromechanical investigations are also important for developing decoding methods by providing the link between neural activities and behaviours. The decoding methods can be used to identify human intent from neural signals generated during a specific behaviour (neural interfaces). Neuromechanical studies have therefore allowed us to address problems ranging from motor control [4], rehabil- itation engineering [5], and human-machine interfaces [6].

Due to the complex interplay of the nervous, muscular, and skeletal systems, researchers have used multiple simulation tools to explore certain aspects of human neuromechanics. A series of models that describe the neuron structures and functions have been proposed, includ- ing the classical Hodgkin-Huxley model [7]. At a macro scale, several analytical and numerical electromyogram (EMG) models have been used to study the electrical outputs produced by skeletal muscles upon activations of populations of neurons [8–11]. In parallel, multiple simu- lation platforms, such as OpenSim [12, 13] and MuJuCo [14, 15], have been widely used to study human body movements during dynamic simulations in biomechanical studies. How- ever, the electrical outputs are rarely simulated together with the biomechanical system during voluntary movements.

OpenSim allows the users to include EMG signals as an additional input for better estimat-
ing the musculo-tendon dynamics and parameters but not to forwardly generate EMG signals
during a movement [12]. Fuglevand's model has been widely used to study muscle force,
motoneuron activities, and EMG signals but is limited to isometric contractions [16]. A recur-
rent neural network has been recently used as a black box to convert motions to the down-
sampled and smoothed muscle activities, without knowing the internal parameters of the
system [17]. This inhibits using the model to study cause-effect relations between system
parameters or to validate algorithms that require internal information. For instance, validating
adaptive decomposition algorithms requires knowing the firing status of each motoneuron.
One primary reason for the lack of an integrated and precise EMG simulator feasible for vol-
untary movements is the absence of models that link EMG generation to movement biome-
chanics. Another reason is that current advanced EMG models are not efficient enough to
adapt to the non-stationary physiological parameters during voluntary movements. These two
challenges have been addressed by our recently released EMG model, BioMime [18], which is
a conditional generative model that takes the physiological parameters as inputs and outputs
the dynamic motor unit action potential (MUAP) signals efficiently.

Using BioMime as a key module, here we propose NeuroMotion, an open-source EMG
generative model that replicates the full-spectrum generation of electrophysiological signals
during voluntary movements. EMG signals are produced by the mixed convolutions of the
MUAPs (electrical potential templates generated when one motoneuron is activated) and
spike trains (continuous discharge timings of a motoneuron) of all active motor units during a
movement. Correspondingly, NeuroMotion consists of three modules. The first module is an
upper-limb musculoskeletal (MSK) model with OpenSim API for defining and visualising the
movement and estimating the muscle fibre lengths and muscle activations during the move-
ment. The muscle fibre lengths are then utilised by the second module, BioMime, to simulate
the dynamic MUAPs during the movement. The third module is a motor unit pool model that
receives the neural inputs derived from the muscle activations and outputs stimulations to
each muscle in the format of spike trains. An overview of NeuroMotion is shown in Fig 1.
With NeuroMotion, users can synthesise large EMG datasets under a vast repertoire of hand
and wrist movements with all intermediate variables available. One potential application is to
fast iterate regression and classification algorithms on a synthetic dataset before collecting
experimental data (this data augmentation example will be used later in the Results). Another
application is the validation of information extraction algorithms (e.g., EMG decomposition
algorithms) when experimental data lacks the ground truth. Data, codes, and instructions are
available at https://github.com/shihan-ma/NeuroMotion.

**Methods**

Overview

NeuroMotion is a surface EMG generative model that is designed to simulate EMG signals
during voluntary human upper limb movements. NeuroMotion takes the kinematics of
human hand, wrist, and forearm as inputs and outputs the synthetic surface EMG signals that
explain the myoelectric activities responsible for the movement. Three modules are incorpo-
rated to complete this full-spectrum simulation, including an upper-limb MSK model with
OpenSim API, BioMime, and a Motor Unit Pool model. Intermediate variables (changes in
physiological parameters, MUAPs, and spike trains) are available in this hierarchical and mod-
ular simulation. Details of the three key modules are described in Section Core modules. An
overview of the workflow of NeuroMotion is shown in Fig 1.

Core modules

**MSK model with OpenSim.** OpenSim is an open-source software platform that is used
for modelling and analysing neuromusculoskeletal systems [12, 13]. With this platform,
<figure><img src="image_5.png" /><figcaption>**Fig 1. Overview of NeuroMotion.** A) NeuroMotion provides an integrated generative model to simulate forearm surface EMG signals during voluntary hand and wrist movements. B) The movements are defined by the kinematics of an upper limb MSK model. Muscle fibre lengths during the movement are estimated by **OpenSim** and imported into **BioMime** to simulate dynamic MUAPs, which are the continuously changing electrical potential templates generated by the activations of the motor units. The other physiological parameters, including the current source propagation velocity and motor unit depth, can be estimated from muscle fibre length and imported to BioMime as well. Neural inputs to each muscle, which can be derived from the normalised muscle activations from OpenSim or set by the users, are further transformed into spike trains by the **Motor Unit Pool** model. Surface EMG signals are finally calculated as the summation of the convolution of the dynamic MUAPs and the spike trains among all active motor units.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012257.g001

researchers can build MSK models, simulate dynamic movements, and perform motion analysis for biomechanical studies that advance movement science. Here, we use OpenSim as a bridge between movements and physiological system states for three purposes: (1) to define and visualise movements of MSK models, (2) to track the changes in muscle fibre lengths during the movement, and (3) to predict the muscle activations across muscles.

Since we mainly focus on simulating EMG signals from the human forearm during hand, wrist, and forearm movements, we chose the ARMs Wrist and Hand Model [19] (ARMs in abbreviation) from the large database of opens-source OpenSim models as the MSK model used in NeuroMotion. The original ARMs model includes 23 degrees of freedom (DoFs), including the full DoFs of finger movements and the flexion/extension and radial/ulnar deviation DoFs for the wrist. We added the pronation/supination DoF to the source files of ARMs model such that a full range of hand, wrist, and forearm movements could be simulated.

The MSK model in NeuroMotion can be customised to the subject's anthropometry in three ways. First, a fully personalised MSK model can be created by acquiring medical images of the MSK system and building physics-based models. This process requires considerable skills and manual interventions, but has been facilitated by the advancement of automatic tools [20, 21] and associated machine learning methods [22]. A more practical approach, which is automated in OpenSim, is to scale each segment in the ARMs model by minimising the distances between the virtual markers placed on ARMs model and the real markers placed on the subject. The ARMs model can also be manually scaled by anthropometric measurements.
The movement defined by the user can be visualised and investigated in the OpenSim GUI. Muscle-tendon lengths during the movement are related to the joint angles and are tracked and extracted using the OpenSim built-in functions. Assuming tendons to be rigid at constant slack length, the muscle fibre lengths are deduced and imported into BioMime to generate the dynamic MUAPs. Furthermore, the individual muscle activations responsible for the movement are estimated with the built-in Static Optimization tool in OpenSim. Static Optimisation solves the inverse dynamics by minimising the cost function of muscle activations under the constraints of muscle activation-to-force conditions. The exported muscle activations are further normalised and imported into the Motor Unit Pool model to act as the neural input to the motoneuron pools, according to the linearity properties of the motoneuron pools [23]. OpenSim's built-in tools have been integrated into NeuroMotion by its Python API.

**BioMime.** After OpenSim captures the changes in a biomechanical system from the joint kinematics, BioMime is used to simulate the corresponding MUAPs during the movement. BioMime is a deep conditional generative model that simulates electrical potential fields given a set of physiological parameters [18]. A schematic of BioMime's training and inference pipeline is shown in Fig 2. Learning from the outputs of its teacher numerical model [10], BioMime captures the relations between the physiological parameters and the MUAP templates and essentially replicates the biophysical properties of the volume conductor of the forearm. When the time-varying parameters during a movement are imported into BioMime, the output will change accordingly. Therefore, BioMime can be used to generate a sequence of MUAP templates during any movement as long as the parameter changes are available and plausible. The computational cost to simulate a movement in high temporal resolution is extremely low given BioMime's ultra-fast inference speed (0.287 seconds per muscle per condition [18]).

The public BioMime model was trained on a dataset generated by a numerical model with a specific forearm anatomy. Therefore, the BioMime model can be regarded as a subject-specific surrogate model of the volume conductor with this forearm anatomy. The user can keep this well-trained model or train their individual BioMime. Detailed instructions on how to train BioMime can be found in https://github.com/shihan-ma/BioMime. The newest version of BioMime supports changes in seven physiological parameters, including the conductivity of the fat layer and the fibre number, depth and medial-lateral position, innervation zone, conduction velocity, and fibre length of each motor unit (Details in the supplementary of [18]). Some of these parameters can be measured experimentally. For example, the location and the fibre length of a motor unit could be estimated by using ultrasound [24, 25]. The location of the innervation zone and the muscle fibre conduction velocity can be estimated, for example by using high-density EMG recordings [26, 27]. Some parameters can also be obtained from the

<figure><img src="image_6.png" /><figcaption>**Fig 2.** **Schematic of BioMime.** BioMime is a conditional generative model trained in an adversarial way. The generator takes the specified conditions and latent features of MUAPs as the inputs and outputs the simulated MUAP signals. The latent features can be encoded from existing MUAPs or sampled from a prior distribution. The discriminator distinguishes real samples from fake samples conditioned on the specified parameters. Only the generator is used in NeuroMotion.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012257.g002
synthetic movement of the MSK model. As described in Section MSK model with OpenSim, muscle fibre length can be estimated by OpenSim. The depth of a motor unit territory and the motor unit conduction velocity can be approximated from the muscle fibre length under some assumptions (see Section Estimate parameter changes).

As a conditional generative model trained in an adversarial manner, BioMime has an encoder-decoder structured generator and a discriminator. During the inference, only the generator is required. BioMime can be used to generate new MUAPs by morphing an existing MUAP or by sampling from a prior distribution. The difference is that by morphing an existing MUAP, the output MUAPs keep the properties of the input MUAP that are separated from the seven physiological parameters, while these properties are sampled from a prior distribution by sampling. When these latent properties are obtained, MUAPs that change under dynamic movements can be simulated by passing the latent and the changing physiological parameters into the decoder. Since the seven parameters explain the bulk of variations of the MUAP templates of the specific subject [8, 9, 28], the difference between the MUAPs generated by morphing or by sampling is small. NeuroMotion supports both of these two strategies and provides a dataset of MUAPs paired with their physiological parameters for generation under the morphing pattern.

**Motor unit pool.** The purpose of including the motor unit pool model in NeuroMotion is twofold: (1) to initialise the properties of the motor units within one muscle, and (2) to simulate the motoneuron activities and generate their spike trains. In NeuroMotion, we implemented two types of motor unit pool models, the classical Fuglevand's model proposed in [16] and a cohort of leaky fire-and-integrate (LIF) neuron models adapted from [29].

The properties of the motor units in the two models are initialised following the rules below. In both models, the motor unit size (number of muscle fibres within each unit) is proportional to the amplitude of the motor unit twitch force [16, 30], and the summation of the motor unit sizes meets the expected total number of fibres in one muscle (Eq 1). In the classical model, the peak twitch force is exponentially distributed while in the LIF-based model, the twitches are distributed following the linear-exponential function in Eq 2. Eq 2 was derived from experimental measurements in the literature on human forearm muscles [16, 31, 32] following the method described in [29]. The total number of fibres in one muscle was estimated as the ratio between muscle physiological cross-sectional area and the average fibre cross-sectional area for typical human forearm muscles [16, 33]. In the classical Fuglevand's model, the conduction velocity is normally distributed within the common range and then sorted from small to large to be positively related to the motor unit size. In the LIF-based model, we adapted the mathematical relationships between the axonal conduction velocity and the neuron surface area from [34] to estimate the fibre conduction velocity. In both models, the depth and medial-lateral positions of motor units are uniformly distributed within the muscle territories. The innervation zone and fibre length are normally distributed within the common ranges.

$$MN_{size} = \frac{\overline{f^{tw}}(j)}{\sum_{k=1}^{N} \overline{f^{tw}}(k)} \cdot N_f, j \in [1, N] \quad (1)$$

where $\overline{f^{tw}}(j)$ is the estimated twitch force of the jth motor unit in the population, N the number of motor units, and $N_f$ the total number of fibres in one muscle.

$$\overline{f^{tw}}(j) = 0.81 \cdot (18.51 \cdot \frac{j}{N} + 104.10^{\frac{j^{4.83}}{N}}), j \in [1, N] \quad (2)$$
In both models, the motoneurons innervating a muscle are sequentially recruited in terms of the intensity of the neural input. The smallest motor unit is first recruited with the lowest impulse response amplitudes. Once recruited, the motoneuron starts to fire regularly. The fir-ing rates of the recruited motor units increase with the input to the motoneuron pool until the peak firing rates are reached. The two models both predict the specific discharge behaviour (spike trains) of the motoneurons from the intensity of the neural input. Specifically, the firing state of a motor unit is decided by the recruitment threshold of the motor unit in the classical Fuglevand's model [16]. The recruitment thresholds are exponentially distributed such that few motor units have high thresholds and many motor units have low thresholds. In the LIF model-based approach, each LIF model captures the complex nonlinear MN dynamics by using a parallel combination of a leaky resistor and a capacitor. A spike is generated when the transmembrane voltage reaches a certain threshold, followed by a refractory during the action potential. With physiologically realistic distributions of the electrophysiological properties and a further interpolation between the motoneuron-specific electrophysiological properties sup-ported by experimental data [29] in the LIF model, the predictions of both models are consis-tent with the onion skin theory [35] and Henneman's size principle of sequential motoneuron recruitment [36].

**Toolbox functions**

In this section, we introduce the basic functions provided by NeuroMotion to define the movement of the ARMs model (Section Define movement), track the changes in physio-logical parameters (Section Estimate parameter changes), set common drives to motoneu-ron pools (Section Set neural inputs to motoneuron pools), and define the structure of motor unit pools (Section Configure motor unit pools). Python code examples for each function are displayed.

**Define movement.** Movements are simulated by using the ARMs model in NeuroMotion. A python class *MSKModel* was implemented to handle the related functions. NeuroMotion supports three approaches to defining the movement of the ARMs model. The most basic way is by assigning joint angles to each DoF by using the function *'load*(*Fig 3A*)*. The joint angles can be obtained from motion capture data (like in [37]) or be measured from sensors (as by data glove and angle sensors). There are 24 DoFs in the ARMs model and each can be changed separately within a predefined range. NeuroMotion automatically checks whether each input joint angle is within the correct range. It is possible that a movement defined by a set of joint angles is not feasible even though each joint angle is reasonably assigned. For example, a combination of joint angles may result in a sudden and unlikely change in muscle fibre lengths. We encourage the users to visualise the movement in OpenSim GUI and to check if there are such aberrations. The motion file required by OpenSim can be generated from joint angles by using the function *'write哲' in NeuroMotion (Python code example at Line 10 in Code Block 1).*

The second way to define a movement is by interpolating between predefined poses. Neuro-Motion provides eight default poses, including hand open/grasp, wrist flexion/extension, wrist radial/ulnar deviation, and forearm pronation/supination (Fig 3B). Users can define their cus-tomised poses by creating new poses in OpenSim GUI. A movement can be simulated by set-ting one pose at each stage from the predefined poses and setting the durations of each transition (Python code example at Lines 1-5 in Code Block 1). The eight default poses can be combined before interpolation by concatenating the names of the poses, e.g., 'open+flex' means hand open and wrist flexion at the same time. The joint angles are summed after the combination. Users may want to define the movement directly from the motion capture data,
<figure><img src="image_7.png" /><figcaption>**Fig 3.** Illustration of the function of movement definition in NeuroMotion. The MSK models are driven by the joint angles in OpenSim. NeuroMotion provides three ways to define a movement. A) The most basic way to define a movement is by setting the joint angles directly. In this example, the ARMs model is driven by the metacarpophalangeal (MCP) and proximal interphalangeal (PIP) joints of the five digits. B) Eight default poses are provided. Users can combine these poses and interpolate between them to simulate a smooth movement. C) Motion capture data can be used to drive the MSK model. Markers should be aligned between experimental and simulation settings.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012257.g003

as the third way displayed in Fig 3C. This is possible with OpenSim GUI. The users need to place virtual markers at the same anatomical positions on the MSK model as the experimental positions. Then the MSK model can be scaled and the movement can be driven by the motion capture data.

```python
# build msk model and simulate a movement
msk = MSKModel()
poses = ['default', 'default+flex', 'default']
durations = [2] * 2  # 2 seconds for each transition
fs = 5  # 5 Hz
ms
msk.sim概述 fs, poses, durations
ms_lens = msk mov2len(ms labels=ms labels)
changes = msk.len2params()
msk.概述 .mov ('./res/mov.mot')
```

**Code Block 1.** Example python codes for defining a wrist flexion movement by interpolation and for acquiring changes in physiological parameters.

**Estimate parameter changes.** BioMime can simulate MUAPs during the changes in seven physiological parameters. Users can import the parameters measured during an experiment. If these physiological parameters are not available from the experiment, NeuroMotion provides an alternative method to approximate some of the parameters. Muscle fibre lengths during a movement of the MSK model can be estimated by using the Muscle Analysis Tool in
OpenSim. Two more parameters can be approximated from the muscle fibre length changes. Under the assumption of constant muscle volumes during a movement [11], the cross-sectional area of a muscle changes inversely with the muscle fibre length. Then the conduction velocity can be tracked since it is positively correlated with the cross-sectional area of muscle volumes. If the location of a motor unit consistently changes with the radius of the cross-section, the depth of the motor unit can be estimated as well. Therefore, changes in three physiological parameters can be estimated during a movement. These changes can be obtained by using the function of 'len2params' in NeuroMotion (Python code example at Line 9 in Code Block 1).

**Set neural inputs to motoneuron pools.** Activations that trigger muscle contractions are often represented by spike trains. The spike trains are generated by motoneurons by transforming the input they receive and are convolved with the MUAP templates to generate the interference EMG signals. It still remains a challenge to decode the real neural input to motoneurons from experimental EMG data, and thus it is difficult to set these inputs to their true values during a movement.

NeuroMotion provides three ways to set the input to each motoneuron pool. First, the user can set predefined and commonly used activation profiles, which include constant, trapezoidal, triangular, and sinusoidal activations. The amplitude and duration of each type of activation can be easily changed. Second, as described in Section MSK model with OpenSim, muscle activations during a movement can be estimated by OpenSim with the Static Optimisation Tool. The muscle activations are then normalised to the activations during the maximum voluntary contractions and used as the neural inputs to the muscles. Third, NeuroMotion provides the interfaces for the user to set the muscle activations by the outputs of their algorithms, for example, using the normalised EMG signals recorded during experiments or the signals derived from muscle synergies [38].

**Configure motor unit pools.** NeuroMotion organises the motor unit pool structures by providing a python class *MotoneuronPool*. This class parameterises the electrical and mechanical properties of a motor unit pool, of which each property can be easily customised. Fundamental properties include the motor unit size, conduction velocity, motor unit positions, innervation zones, and fibre lengths. In the implementation of classical Fuglevand's model, the motor unit recruitment threshold, the minimum and maximum firing rate of each motor unit, variability of the inter-pulse intervals, and motor unit twitch force and contraction time can be manually assigned. In the LIF-based model, the motoneuron-specific parameters, such as the resistance, capacity, and time constant, are set given the mathematical descriptions in [29]. Example Python codes that define a classical motor unit pool model, define the input to motoneurons and calculate the spike trains are shown in Code Block 2.

Users may prefer using the classical Fuglevand's model as it is more conceptually simple with fewer parameters. The LIF-based model might be preferred as it provides a more physiological way to predict the discharge activities based on the physiological dynamics of motoneuron membrane charging and discharging and allows a finer control of motoneuron properties by providing more controllable parameters. Users could also add their self-developed models, for example the Izhikevich neuron model [39], into NeuroMotion.

## Results

NeuroMotion provides the input, output, and all intermediate variables to the users, including the joint kinematics and the muscle fibre lengths during a movement, the neural input to each
```python
# test motor unit pool model
mn_pool = MotoneuronPool(num_mu, rr, rm, rp, pfr1, pfrd, mfr1, mfrd,
    ge, c_ipi, frs1, frsd)          # parameters that define the motor unit
    recruitment and firing pattern
# properties
config = edict({
    'numFB': 25000,
    'depth': [20, 30],
    'angle': [20, 30],
    'iz': [0.5, 0.1],
    'len': [0.5, 0.1],
    'cv': [4, 0.5]
})
properties = mn_pool Assign_properties(config, True)
# define neural input
fs = 2048          # Hz
duration = 6       # s
ext = np.linspace(0, 1.0, fs * duration)
# Force and twitches
mn_pool.init_twitches(fs)
mn_pool.init_quisistatic_ef_model()
# spike trains
_, spikes, fr, ipis = mn_pool(generate_spike_trains(ext))
```

<caption>**Code Block 2.** Example python codes for defining a motor unit pool model and for setting the neural input to the pool.</caption>

motoneuron pool, the neural commands to the muscles in the format of spike trains, the dynamically changing MUAPs, and the interference surface EMG signals. Given the above information, NeuroMotion can be used to analyse the impact of a specific parameter on the EMG signals and to provide synthetic datasets for data augmentation and validations of EMG-related algorithms. Here we displayed the variations of MUAPs during dynamic contractions and compared the similarities between MUAPs within and across muscles. A case study was provided in which the synthetic dataset generated by NeuroMotion was used to augment an experimental dataset for improving the accuracy in estimating joint angles from EMG signals.

## Synthetic MUAPs during dynamic contractions

**MUAPs when one, two, and three parameters change.** NeuroMotion allows researchers to estimate the changes in muscle fibre length, conduction velocity, and motor unit depth during the movement of the ARMs model (Sections MSK model with OpenSim and Estimate parameter changes). Here, we gave examples of the simulated MUAPs during a wrist flexion and extension movement (Fig 4A) with four levels of parameter changes, including changing only muscle fibre length, changing muscle fibre length and conduction velocity, changing muscle fibre length and motor unit depth, and changing all three parameters. One representative MUAP from the ulnar head of Flexor carpi ulnaris (FCU(u)) and its changes over time are displayed in Fig 4B. By changing the muscle fibre length, the duration of the MUAPs slightly varied. Changing the conduction velocity also influenced the duration of the MUAP while changing the motor unit depth had a more substantial impact on the waveforms. Changes due to the three physiological parameters during the movement are visualised in Fig 4C, 4D and 4E, respectively.

**MUAPs across poses.** We then proceeded to study how MUAPs from different muscles change across movements. Three basic movements were simulated by using the tools described in Section Define movement, including sequences of hand-open to hand-grasp, hand-open to wrist flexion and extension, and hand-open to radial and ulnar deviations. We
<figure><img src="image_8.png" /><figcaption>**Fig 4.** Changes in the MUAPs of one representative motor unit in FCU(u) during a wrist flexion/extension movement. A) Sequence of the movement. B) MUAPs generated by NeuroMotion. The MUAPs in the first row were morphed using the muscle fibre length profiles output by the ARMs model using OpenSim, the second row using muscle fibre length and conduction velocity, the third row using muscle fibre length and motor unit depth, and the fourth row using all three parameters. C) Normalised muscle fibre length profiles from the ARMs model. D) Normalised conduction velocity. E) Normalised motor unit depth.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012257.g004

visualised the changes in MUAPs from six superficial muscles. As shown in Fig 5, only the lengths of digitorum muscles (extensor digitorum, ED and flexor digitorum superficialis, FDS) changed during hand grasp and open movement. As a result, MUAPs from ED and FDS displayed large variations while MUAPs from the other muscles maintained their waveforms constant. All six muscles were activated during the flexion and extension movement. Therefore, the duration, amplitude, and waveform of the MUAPs from the six muscles all changed consistently with the movement. The digitorum muscles (ED and FDS) and palmaris longus (PL) muscle contributed less to the radial and ulnar deviations, and thus MUAPs from these three muscles showed fewer variations during wrist deviations.

**In-muscle and cross-muscle MUAP similarities.** Muscles within the forearm have small volumes. Thus, MUAPs within a single forearm muscle might have similar waveforms due to the overlapping motor unit territories. We studied the similarities between MUAPs within a single muscle, across pairs of muscles, and across joint angles during a wrist flexion and extension movement. The similarity was evaluated by the normalised mean square error (NMSE), which is defined by the mean square error between two MUAPs divided by their averaged
<figure><img src="image_9.png" /><figcaption>**Fig 5.** Changes in MUAPs during three movements. One representative MUAP was chosen from each of the six muscles and was consistently morphed during three movements. The first row: normalised muscle fibre length profiles from the ARMs model during grasp-open, flexion-extension, and radial-ulnar deviation. The second to the seventh rows: MUAPs transformed during the three movements in the six muscles. The colours from dark blue to light green indicate the time frames from the beginning to the end of the movement. For each MUAP, a subgrid of three channels were chosen for display, which was consistent across the three movements.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012257.g005

power. The MUAPs were first cropped to a sub-grid where all channels had the maximum amplitude above 75% of the average across channels.

The similarities between MUAPs of three muscles, an extensor (Extensor carpi radialis brevis, ECRB) and two heads of one flexor (Flexor carpi ulnaris ulnar and humeral heads, FCU(u) and FCU(h)), are shown by the confusion matrices in Fig 6. MUAPs within the same muscle generally showed higher similarity than MUAPs across different muscles. MUAPs within the two heads of the FCU muscle were more similar than MUAPs within the FCU and a more distant muscle ECRB. Different MUAPs within the same muscle (off-diagonal) could have similar waveforms since multiple combinations of physiological parameters can contribute to similar MUAPs and different MUAPs may overlap in their territories due to the small muscle volume in the forearm.

We were also interested in how movements influence the similarities of the MUAPs and whether the changes in the similarities during the movement are consistent across MUAPs within each muscle. Fig 7 shows the similarities between MUAPs morphed during a wrist flexion/extension movement and their baseline shapes (MUAPs at zero joint angles). MUAPs were gradually morphed away from their baseline shapes when the joint flexion or extension angle increased. The similarities of MUAPs across joint angles were consistent within the
<figure><img src="image_10.png" /><figcaption>**Fig 6.** Similarity and dissimilarity between MUAPs within/across muscles in confusion matrices. Deeper colours in the first and the third columns indicate a lower normalised mean square error (NMSE) and higher similarity. Black denotes NMSE > 0.2 (less similar) and white NMSE < 0.2 (more similar) in the second and the fourth columns. The left two columns show similarity between MUAPs within single muscles while the right two columns show similarity between MUAPs across two muscles.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012257.g006

MUAPs of each muscle with small variances. MUAPs in extensors showed higher levels of dis-similarity during wrist extension than during wrist flexion, and the opposite for MUAPs in flexor muscles.

## Synthetic EMG signals for data augmentation

**Simulation process.** The synthetic EMG signals generated by NeuroMotion can be used to augment experimental datasets for myoelectric control. In this preliminary case study, we simulated hand and wrist movement following the experiments in [40]. In the experiments, six subjects performed two motor tasks with self-selected speeds, including wrist flexion/extension and MCP flexion/extension movements. Five trials were performed for each task. The duration of each trial was 15 seconds. Six EMG sensors (Delsys Trigno Wireless, Delsys Inc.) were placed on the six forearm muscles according to SENIAM guidelines [41]. The six muscles include flexor carpi radialis (FCR), flexor digitorum superficialis (FDS), flexor carpi ulnaris (FCU), extensor carpi ulnaris (ECU), extensor digitorum (ED), and extensor carpi radialis longus (ECRL). Proper positioning of each EMG sensor was chosen by physically palpating the muscle during sustained isometric contraction and visually confirming the EMG signal. The surface EMG signals were recorded at 2000 Hz and the wrist and MCP joint angles at 40 Hz.

The EMG signals were normalised by the maximum EMG signals during maximum voluntary contractions and then used as the neural inputs to the corresponding muscles to produce spike trains. The joint angles were used to drive the ARMs model, where the output muscle fibre lengths were imported to BioMime to simulate dynamic MUAPs. Note that the latent representation was initialised for each motor unit by randomly sampling from the Normal
<figure><img src="image_11.png" /><figcaption>**Fig 7.** Similarity between MUAPs across joint angles during a wrist flexion/extension movement. The joint angles were normalised between -1.0 (flexion) and 1.0 (extension). MUAPs at zero flexion/extension were used as the baseline shapes, which were compared with the continuously morphed MUAPs during the movement to calculate the NMSE (dissimilarity metric, higher values indicate higher differences). NMSEs were averaged within each muscle at each time step, with the variances shown by the green shaded area.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012257.g007

distribution and then fixed during the simulation. This ensures that the MUAPs are only changed by the physiological parameters during the movement. Finally, the synthetic EMG signals were simulated by convolving the spike trains and the MUAPs. We simulated a total of 320 channels of surface EMG for all trials of the six subjects. The simulated channels are uniformly placed around the forearm and cover the majority of the forearm muscles. Six channels were selected for each subject by matching the coordinates between the experimental electrodes and the simulated electrodes. We also checked the amplitudes of the selected channels during specific movements. For instance, the signals collected by channels placed on the flexors should have a higher amplitude during wrist flexion compared to the amplitude during wrist extension. Root mean square (RMS) values were computed from 200-ms intervals with 50-ms overlapping. Each RMS segment was labelled with the two joint angles.

**Case studies of data augmentation.** We investigated whether the simulated signals have practical uses, for example, to augment the experimental dataset in this case study. The assumption here is that even though the modules used in the simulation were not customised to each subject, given the common patterns of muscle activation during movements, the simulated EMG signals could capture the general characteristics of each movement that would be useful for regressing the joint angles. We show that the synthetic EMG signals are of practical use in two ways.

First, we trained ridge regressors on the synthetic dataset from each trial of each subject. The trained regressors were directly applied to regress the joint angles from the experimental data of the same subject and the same trial. We observed that for each of the subjects, there was at least one trial, in which the joint angles predicted by the regressors trained on the synthetic data were highly consistent with the experimental joint angles. The Pearson correlation coefficients (PCC) are 0.52 ± 0.14 (mean and standard deviation across subjects) and 0.61 ± 0.14 for regressing the wrist and MCP angles, respectively. Taking one subject as an example, we found consistent regression results when the simulation was repeated five times. As might be expected, we found that the data produced by NeuroMotion closely matched the real data when the regressors trained on the synthetic data performed well on the real data. An example of the regression results is shown in Fig 8. To account for the fact that the joint angles mainly consist of low-frequency components, we estimated a lower bound of such predictions by comparing the predictions with the randomly shuffled ground truth labels. An extreme case
<figure><img src="image_12.png" /><figcaption>**Fig 8.** Representative results of using NeuroMotion's output to train regression algorithms. A) The regressors trained on the synthetic dataset could predict wrist and MCP joint angles from the experimental data. The ground truth joint angles are displayed in black and the predictions in blue. B) The simulated six-channel EMG signals (blue) and the experimental signals (black) of the same trials in A. C) the RMS extracted from the simulated signals (blue) and the experimental signals (black) of the same trials in A.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012257.g008

is that if the joint angles are constant, shuffling the labels of joint angles will not affect the pre-
diction accuracy. The random shuffling was repeated ten times for each trial. The PCC
between the prediction and the randomly shuffled angles has a large variance and contains
both negative and positive values, resulting in an averaged PCC around zero (−0.14 ± 3.68).
The PCC between the predictions and the ground truth joint angles is significantly higher than
the PCC between the prediction and the randomly shuffled angles following the Student's t-
test (p < 0.005).

Second, we show that the synthetic data can be potentially used to augment the experimen-
tal data to improve the regression accuracy. Specifically, we randomly selected two trials of
experimental data and two trials of synthetic data to form a new training dataset. The ridge
regressors trained on this augmented dataset were tested on one trial of the experimental data,
which was different from the trials in the training dataset. For all six subjects, we were able to
augment the experimental data with data generated by NeuroMotion such that the resultant
regressor performed better than a regressor trained on the experimental data alone. The per-
formance of the regressor was significantly improved when trained on these augmented data-
sets (Student t-test with *p* < 0.05). The PCC of training on the augmented dataset versus
training on the full experimental dataset is 0.77 ± 0.08 versus 0.75 ± 0.08 for wrist and
0.72 ± 0.06 and 0.66 ± 0.07 for MCP.

Discussion

We have presented NeuroMotion, which provides the first open-source simulator for neuro-
mechanical investigations by modelling the electrical outputs during voluntary human move-
ments. With the three key modules, NeuroMotion provides a wealth of functions that give full
freedom to users to simulate signals by using the default settings or by customising the config-
urations. All input, intermediate, and output variables are available in this hierarchical and
full-spectrum simulation, indicating multiple potential usages of NeuroMotion.

Rich variances of MUAPs during dynamic contractions

NeuroMotion provides three ways to define the movement of the ARMs model. The easiest
way is to interpolate between (the combinations of) the eight default poses. Examples in Figs 4
and 5 showed simulations of MUAPs by using the interpolation tools. Physiological parame-
ters were smoothly changed during the movement, which resulted in a continuous variation of
MUAP waveforms. The variations of synthetic MUAPs are similar to the variations of MUAPs
observed during experiments in a qualitative way [42–44]. For example, the shortening of
muscle fibres or increase of conduction velocity during natural muscle contractions reduced
the duration of a MUAP; MUAPs of muscles that contribute less to a movement showed fewer
variations (PL in wrist deviation).

Compared with the limited muscles studied under a few discretised joint angles during experiments, NeuroMotion can simulate MUAPs from eight forearm superficial muscles under any voluntary movements in a high temporal resolution. Such functionality allows users to test the assumptions made in the simulations (e.g., changes in parameters) or study the rich variances of MUAPs within or across muscles as displayed in Fig 6. The ability of NeuroMotion to generate dynamic MUAPs consistent with muscle functionalities demonstrates its superiority to previous models that either adapted the classical cylindrical model to discretised stages [45] or empirically transformed the MUAP waveforms [46]. We expect that NeuroMotion will meet the demand for providing synthetic EMG signals during dynamic muscle contractions for validating adaptive decomposition algorithms [42]. For example, NeuroMotion was used in [47] to generate data during wrist flexion and extension at constant forces. The simulated data were used to test the online learning EMG decomposition algorithm.

We also showed that the generated MUAPs were muscle-specific with low similarity across different muscles (Fig 6). The variations of MUAPs during a movement were consistent within muscle but different across muscles (Fig 7). When the wrist flexion/extension angle increased, MUAPs were gradually transformed and deviated away from their initial version, leading to an increasing difference (higher NMSE). The similarity at each joint angle showed small variances across MUAPs within each muscle, which indicates that deformations of fibres within one muscle are similar. MUAPs from flexors experienced more variations during flexions while MUAPs from extensors showed more variations during extensions. This observation fits well with the fact that flexors/extensors manifest a more drastic contract during flexions/extensions. These results evidence the potential of using NeuroMotion to investigate the variations of MUAPs within/across muscles during voluntary movements.

Potential as a test-bed for cause-effect studies and data augmentation for regression analysis

With access to all inputs, outputs, and internal parameters in NeuroMotion, users can study cause-effect relationships between the joint kinematics, neural inputs, physiological parameters, muscle forces, MUAPs, spike trains, and surface EMG signals. As a simple demonstration of such studies, we showed how MUAP changes differently with one to three physiological parameters changed during a wrist flexion/extension movement (Fig 4). Coupled with experimental MUAPs observed at the same muscle under discretised movements, such simulations will help test the assumptions made in Section Estimate parameter changes and elucidate the contributions of the physiological parameters to the MUAPs. Another example displays the transformations of MUAPs under three wrist-and-hand movements (Fig 5). This allows the investigation of how MUAP changes when the commonly used forearm DoFs are activated, which will further provide insights into designing and improving current adaptive decomposition algorithms.

NeuroMotion can also be used to study the “inverse” relations between the outputs and the
inputs of the system, for example, to pre-train algorithms that regress joint angles from EMG
signals and to augment the experimental dataset. Here we provide a preliminary study to show
that ridge regressors trained on the synthetic dataset can be potentially applied to regress joint angles from experimental signals. Furthermore, we observed that the performance of some regressors could be improved when trained on the augmented dataset. These findings are consistent with the perspective of using synthetic datasets to augment the limited actual dataset and facilitate algorithm training [10, 48]. One important observation is that the effectiveness of data augmentation is correlated with subjects. This is expected, as the BioMime module in NeuroMotion only captures the geometry of one specific anatomy from the database, which may deviate from the anatomy of the specific subjects. We discuss this limitation in the next section in more detail.

We also observed that the variance of the augmentation performance across trials for one subject could be large. One possible reason is that we used the normalised EMG signals recorded on each muscle as the neural input, which largely affects the properties of the synthetic EMG signal. Electrode shifts across trials, misalignment between the electrodes and muscles, and crosstalk due to the volume conduction all make the accurate estimation of neural inputs challenging. Moreover, the perspective that motor neurons are controlled in the group of muscles is being challenged by new experimental findings that motor neurons are controlled in clusters that might distribute across muscles [49]. In this case, it is even more difficult to accurately group motor neurons into clusters and estimate the neural inputs to the clusters. Even though, we found that synthetic signals of some channels (muscles) constantly show a higher correlation with their experimental counterparts. This can be explained by the muscle's location, i.e., some muscles are easier to find by palpation. For example, for channels that collect signals from ED and ECRL, the PCC is above 0.7 for wrist regression and above 0.90 for MCP regression. The PCC is lower for slender muscles that are difficult to separate from other muscles (0.40 for PL muscle for wrist angle regression). All these uncertainties of experimental signals relate to the general reasoning behind building forward models like NeuroMotion. The flexibility provided by NeuroMotion allows the investigation of the actual relationship between the neural input and the joint kinematics, for instance, by adjusting the neural inputs to make the synthetic EMG signals more similar to the experimental data under the same protocols. This will in turn give insights into the real relations between joint kinematics and neural inputs during voluntary movements.

## Limitations and future directions

As the first integrated open-source generative model to simulate EMG signals during voluntary movements, NeuroMotion still has some limitations that should be noted. With its modular structure, the intrinsic limitations of NeuroMotion come from the three modules, of which the error is trackable but aggregates in the computation of the predicted EMG signals.

First, the current BioMime model in NeuroMotion is used as a replica of a specific volume conductor but not a universal one. Therefore, there are errors introduced when NeuroMotion is used to simulate signals for specific subjects. A related limitation is that the locations of the electrodes are fixed (10 × 32 channels around the forearm) and cannot be specified. This can be addressed by training a subject-specific BioMime with different electrode configurations. Building a universal BioMime model that can be adapted to specific subjects and allows random samplings of electrical potentials on the surface is also in our future work.

Second, the accuracy of the estimated muscle fibre length will affect MUAP generation, as the muscle fibre lengths are the inputs to the BioMime module to indicate physiological parameter changes. Forearm muscle fibre length is very difficult to estimate experimentally, although some algorithms can predict lower limb muscle fascicle length by tracking the endpoints in the recorded ultrasound images [50]. In NeuroMotion, the muscle fibre lengths are
derived from the geometrical muscle-tendon lengths computed with the ARMs model in OpenSim. Under the assumption that the tendons are rigid and at constant slack length [51], the method does not account for the dynamic changes in muscle fibre lengths to achieve the physiological force equilibrium between the muscle and the extensible tendon. The simplification of neglecting the compliance of tendons, taken for computational convenience and modelling simplicity, may be a source of error when estimating muscle fibre lengths, as quantified in [51], and discussed for forearm muscle-tendon systems in [52].

Third, the accuracy of the spike trains estimated from the Motor Unit Pool module is affected by both the neural inputs to the module and the module itself. The neural inputs estimated from the muscle activations computed by the OpenSim's Static Optimisation tool are sensitive to the aforementioned inaccuracies in muscle fibre length [13] and the intrinsic limitations of Hill-type muscle modelling [53]. This includes multiscale modelling simplifications [54], physiological and numerical instabilities [55], and limited phenomenological modelling of the muscle's activations dynamics [53, 56, 57]. The Static Optimisation method also computes muscle activations, and hence, neural inputs without considering the concurrent dynamics of the proprioceptive, spinal and propriospinal circuits and reflexes during voluntary movement. These concurrent dynamics can be modelled and were demonstrated to have effects in the simulation of neuromuscular and musculoskeletal function [58-63]. The deafferented approach taken in the present study does not consider the sensory feedback related to passive changes of fibre length, which would result in inaccuracies in the estimated muscle fibre lengths and muscle activations and thus simulated EMG as well. The EMG-informed estimation of neural inputs (used in Section Synthetic EMG signals for data augmentation) is sensitive to the normalisation and the experimental settings, for example, the electrode placement and crosstalk [64]. Inaccuracy in the estimated spike trains also comes from the limitations of motor unit pool models. The two models used in NeuroMotion are designed for isometric contractions and do not consider changes in motor neuron discharge patterns associated with voluntary movement [16, 29, 44]. These models are also not complex enough to physiologically describe the intrinsic discharge properties and mechanisms of motor neurons (e.g., the nonlinear transformation of common synaptic inputs into total dendritic membrane current discussed in [65]).

Even with all the limitations of each module, the modular structure of NeuroMotion allows it to be improved with the emergence of new measuring and modelling approaches in the EMG field. For example, *in vivo* measurements of muscle activations, muscle fibre length, and other physiological parameters would contribute to a more accurate estimation of neural inputs and MUAPs. More complex modellings of motor neuron pools could account for nonlinearities in spiking dynamics and afferent feedback activities. We would expect NeuroMotion to be frequently updated with the advancements in each module and be validated in future by comparing the synthetic signals with the experimental data.

## Conclusion

We proposed NeuroMotion, an open-source EMG simulator that generates surface EMG signals during human forearm movements. For the first time, NeuroMotion provides a core resource for the neuromechanics community to address the problem of simulating physiological electrical outputs during biomechanical movements and allows a full-spectrum simulation from movements to neural commands and EMG signals. All intermediate variables (kinematics, dynamic MUAPs, spike trains) are available during the simulation, which makes the model flexible to use for the users' purposes. We demonstrated that one potential application is to use synthetic data to augment the experimental data for training regression algorithms.
We would expect that with the functionality and extensibility provided by NeuroMotion, users
can customise the way they utilise NeuroMotion that will ultimately yield progress in the neu-
romechanics fields.

Author Contributions

**Conceptualization:** Shihan Ma, Dario Farina.

**Data curation:** Shihan Ma, Jiamin Zhao, Kostiantyn Maksymenko, Samuel Deslauriers-Gauthier.

**Formal analysis:** Shihan Ma, Dario Farina.

**Funding acquisition:** Xinjun Sheng, Xiangyang Zhu, Dario Farina.

**Investigation:** Shihan Ma, Irene Mendez Guerra, Arnault Hubert Caillet.

**Methodology:** Shihan Ma, Arnault Hubert Caillet.

**Project administration:** Xinjun Sheng, Xiangyang Zhu, Dario Farina.

**Resources:** Shihan Ma, Dario Farina.

**Software:** Shihan Ma.

**Supervision:** Dario Farina.

**Validation:** Shihan Ma.

**Visualization:** Shihan Ma.

**Writing – original draft:** Shihan Ma.

**Writing – review & editing:** Shihan Ma, Irene Mendez Guerra, Arnault Hubert Caillet, Alexander Kenneth Clarke, Dario Farina.

References

1. Valero-Cuevas, F. Fundamentals of Neuromechanics. Springer; 2016.

2. Ting LH, McKay JL. Neuromechanics of muscle synergies for posture and movement. Curr Opin Neuro-biol. 2007; 17(6):622–628. https://doi.org/10.1016/j.conb.2008.01.002 PMID: 18304801

3. Schwartz AB. Movement: how the brain communicates with the world. Cell. 2016; 164(6):1122–1135. https://doi.org/10.1016/j.cell.2016.02.038 PMID: 26967280

4. Nishikawa K, Biewener AA, Aerts P, Ahn AN, Chiel HJ, Daley MA, et al. Neuromechanics: an integrative approach for understanding motor control. Integr Comp Biol. 2007; 47(1):16–54. https://doi.org/10.1093/icb/icm024 PMID: 21672819

5. Ting LH, Chiel HJ, Trumbower RD, Allen JL, McKay JL, Hackney ME, et al. Neuromechanical principles underlying movement modularity and their implications for rehabilitation Neuron. 2015; 86(1):38–54. https://doi.org/10.1016/j.neuron.2015.02.042 PMID: 25856485

6. Farina D, Jiang N, Rehbaum H, Holobar A, Graimann B, Dietl H, et al. The extraction of neural informa-tion from the surface EMG for the control of upper-limb prostheses: emerging avenues and challenges. IEEE Trans Neural Syst Rehabil Eng. 2014; 22(4):797–809. https://doi.org/10.1109/TNSRE.2014.2305111 PMID: 24760934

7. Hodgkin AL, Huxley AF. A quantitative description of membrane current and its application to conduc-tion and excitation in nerve. J Physiol. 1952; 117(4):500. https://doi.org/10.1113/jphysiol.1952.sp004764 PMID: 12991237

8. Farina D, Rainoldi A. Compensation of the effect of sub-cutaneous tissue layers on surface EMG: a sim-ulation study. Med Eng Phys. 1999; 21(6-7):487–497. https://doi.org/10.1016/S1350-4533(99)00075-2 PMID: 10624744

9. Farina D, Mesin L, Martina S, Merletti R. A surface EMG generation model with multilayer cylindrical description of the volume conductor. IEEE Trans Biomed Eng. 2004; 51(3):415–426. https://doi.org/10.1109/TBME.2003.820998 PMID: 15000373
10. Maksymenko K, Clarke AK, Mendez Guerra I, Deslauriers-Gauthier S, Farina D. A myoelectric digital twin for fast and realistic modelling in deep learning. Nat Commun. 2023; 14(1):1600. https://doi.org/10.1038/s41467-023-37238-w PMID: 36959193

11. Mesin L, Joubert M, Hanekom T, Merletti R, Farina D. A finite element model for describing the effect of muscle shortening on surface EMG. IEEE Trans Biomed Eng. 2006; 53(4):593–600. https://doi.org/10.1109/TBME.2006.870256 PMID: 16602565

12. Delp SL, Anderson FC, Arnold AS, Loan P, Habib A, John CT, et al. OpenSim: open-source software to create and analyze dynamic simulations of movement. IEEE Trans Biomed Eng. 2007; 54(11):1940–1950. https://doi.org/10.1109/TBME.2007.901024 PMID: 18018689

13. Seth A, Hicks JL, Uchida TK, Habib A, Dembia CL, Dunn JJ, et al. OpenSim: Simulating musculoskeletal dynamics and neuromuscular control to study human and animal movement. PLoS Comput Biol. 2018; 14(7):e1006223. https://doi.org/10.1371/journal.pcbi.1006223 PMID: 30048444

14. Wang H, Caggiano V, Durandau G, Sartori M, Kumar V. MyoSim: Fast and physiologically realistic MuJoCo models for musculoskeletal and exoskeletal studies. In: 2022 International Conference on Robotics and Automation (ICRA). IEEE; 2022. p. 8104–8111.

15. Todorov E, Erez T, Tassa Y. Mujoco: A physics engine for model-based control. In: 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems. IEEE; 2012. p. 5026–5033.

16. Fuglevand AJ, Winter DA, Patla AE. Models of recruitment and rate coding organization in motor-unit pools. J Neurophysiol. 1993; 70(6):2470–2488. https://doi.org/10.1152/jn.1993.70.6.2470 PMID: 8120594

17. Schmidt MD, Glasmachers T, Iossifidis I. The concepts of muscle activity generation driven by upper limb kinematics. Biomed Eng Online. 2023; 22(1):63. https://doi.org/10.1186/s12938-023-01116-9 PMID: 37355651

18. Ma S, Clarke AK, Maksymenko K, Deslauriers-Gauthier S, Sheng X, Zhu X, et al. Conditional generative models for simulation of EMG during naturalistic movements. arXiv:2211.01856 [preprint], 2022. Available from https://arxiv.org/abs/2211.01856.

19. McFarland DC, Binder-Markey BI, Nichols JA, Wohlman SJ, De Bruin M, Murray WM. A musculoskeletal model of the hand and wrist capable of simulating functional tasks. IEEE Trans Biomed Eng. 2022; 70(5):1424–1435. https://doi.org/10.1109/TBME.2022.3217722

20. Modenese L, Renault J B. Automatic generation of personalised skeletal models of the lower limb from three-dimensional bone geometries. J Biomech. 2021; 116:110186. https://doi.org/10.1016/j.jbiomech.2020.110186 PMID: 33515872

21. Valente G, Crimi G, Vanella N, Schileo E, Taddei F. nmsBuilder: Freeware to create subject-specific musculoskeletal models for OpenSim. Comput Methods Programs Biomed. 2017; 152:85–92. https://doi.org/10.1016/j.cmpb.2017.09.012 PMID: 29054263

22. Saxby DJ, Killen BA, Pizzolato C, Carty C, Diamond L, Modenese L, et al. Machine learning methods to support personalized neuromusculoskeletal modelling. Biomech Model Mechanobiol. 2020; 19:1169–1185. https://doi.org/10.1007/s10237-020-01367-8 PMID: 32676934

23. Farina D, Negro F. Common synaptic input to motor neurons, motor unit synchronization, and force control. Exerc Sport Sci Rev. 2015; 43(1):23–33. https://doi.org/10.1249/JES.0000000000000032 PMID: 25390298

24. Maitland S, Hall J, McNeill A, Stenberg B, Schofield I, Whittaker R. Ultrasound-guided motor unit scanning electromyography. MMuscle Nerve. 2022; 66(6):730–735. https://doi.org/10.1002/mus.27720 PMID: 36106775

25. Hauraix H, De Monsabert BG, Herbaut A, Berton E, Vigouroux L. Force-length relationship modeling of wrist and finger flexor muscles. Med Sci Sports Exerc. 2018; 50(11):2311–2321. https://doi.org/10.1249/MSS.0000000000001690 PMID: 29933345

26. Zhang C, Dias N, He J, Zhou P, Li S, Zhang Y. Global innervation zone identification with high-density surface electromyography. IEEE Trans Biomed Eng. 2019; 67(3):718–725. https://doi.org/10.1109/TBME.2019.2919906 PMID: 31150334

27. Martinez-Valdes E, Farina D, Negro F, Del Vecchio A, Falla D. Early motor unit conduction velocity changes to high-intensity interval training versus continuous training. Med Sci Sports Exerc. 2018; 50(11):2339–2350. https://doi.org/10.1249/MSS.0000000000001705 PMID: 29975302

28. Merletti R, Conte LL, Avignone E, Guglielminotti P. Modeling of surface myoelectric signals. I. Model implementation. IEEE Trans Biomed Eng. 1999; 46(7):810–820. https://doi.org/10.1109/10.771191 PMID: 10396899

29. Caillet AH, Phillips AT, Farina D, Modenese L. Estimation of the firing behaviour of a complete motoneuron pool by combining electromyography signal decomposition and realistic motoneuron modelling. PLoS Comput Biol. 2022; 18(9):e1010556. https://doi.org/10.1371/journal.pcbi.1010556 PMID: 36174126
30. Caillet AH, Phillips AT, Farina D, Modenese L. Motoneuron-driven computational muscle modelling with motor unit resolution and subject-specific musculoskeletal anatomy. PLoS Comput Biol. 2023; 19(12): e1011606. https://doi.org/10.1371/journal.pcbi.1011606 PMID: 38060619

31. Romaiguère P, Vedel JP, Pagni S, Zenatti A. Physiological properties of the motor units of the wrist extensor muscles in man. Exp Brain Res. 1989; 78:51–61. PMID: 2591518

32. Riek S, Bawa P. Recruitment of motor units in human forearm extensors. J Neurophysiol. 1992; 68 (1):100–108. https://doi.org/10.1152/jn.1992.68.1.100 PMID: 1517816

33. Feinstein B, Lindegård B, Nyman E, Wohlfart G. Morphologic studies of motor units in normal human muscles. Cells Tissues Organs. 1955; 23(2):127–142. https://doi.org/10.1159/000140989 PMID: 14349537

34. Caillet AH, Phillips AT, Farina D, Modenese L. Mathematical relationships between spinal motoneuron properties. Elife. 2022; 11:e76489. https://doi.org/10.7554/eLife.76489 PMID: 35848819

35. De Luca CJ, Hostage EC. Relationship between firing rate and recruitment threshold of motoneurons in voluntary isometric contractions. J Neurophysiol. 2010; 104(2):1034–1046. https://doi.org/10.1152/jn. 01018.2009 PMID: 20554838

36. Henneman E. Recruitment of motoneurones: the size principle. Prog Clin Neurophysiol. 1981; 9:26–60.

37. Sgambato BG, Hasbani MH, Barsakcioglu DY, Ibáñez J, Jakob A, Fournelle M, et al. High performance wearable ultrasound as a human-machine interface for wrist and hand kinematic tracking. IEEE Trans Biomed Eng. 2023.

38. Zhao J, Yu Y, Wang X, Ma S, Sheng X, Zhu X. A musculoskeletal model driven by muscle synergy-derived excitations for hand and wrist movements. J Neural Eng. 2022; 19(1):016027. https://doi.org/ 10.1088/1741-2552/ac4851 PMID: 34986472

39. Izhikevich EM. Simple model of spiking neurons. IEEE Trans Neural Netw. 2003; 14(6):1569–1572. https://doi.org/10.1109/TNN.2003.820440 PMID: 18244602

40. Zhao J, Yu Y, Sheng X, Zhu X. Consistent control information driven musculoskeletal model for multiday myoelectric control. J Neural Eng. 2023. https://doi.org/10.1088/1741-2552/acef93 PMID: 37567218

41. Hermens HJ, Freriks B, Disselhorst-Klug C, Rau G. Development of recommendations for sEMG sensors and sensor placement procedures. J Electromyogr Kinesiol. 2000; 10(5):361–374. https://doi.org/ 10.1016/S1050-6411(00)00027-4 PMID: 11018445

42. Glaser V, Holobar A. Motor unit identification from high-density surface electromyograms in repeated dynamic muscle contractions. IEEE Trans Rehabil Eng. 2018; 27(1):66–75. PMID: 30571641

43. Kramberger M, Holobar A. On the prediction of motor unit filter changes in blind source separation of high-density surface electromyograms during dynamic muscle contractions. IEEE Access. 2021; 9:103533–103540. https://doi.org/10.1109/ACCESS.2021.3099015

44. Oliveira AS, Negro F. Neural control of matched motor units during muscle shortening and lengthening at increasing velocities. J Appl Physiol. 2021; 130(6):1798–1813. https://doi.org/10.1152/japplphysiol. 00043.2021 PMID: 33955258

45. Glaser V, Farina D, Holobar A. Simulations of high-density surface electromyograms in dynamic muscle contractions. In: 2017 39th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). IEEE; 2017. p. 3453–3456.

46. Chen C, Ma S, Sheng X, Farina D, Zhu X. Adaptive real-time identification of motor unit discharges from non-stationary high-density surface electromyographic signals. IEEE Trans Biomed Eng. 2020; 67 (12):3501–3509. https://doi.org/10.1109/TBME.2020.2989311 PMID: 32324538

47. Mendez Guerra I, Barsakcioglu DY, Farina D. Wearable neural interfaces: real-time identification of motor neuron discharges in dynamic motor tasks. bioRxiv. 2024; p. 2024–02. Available from https:// www.biorxiv.org/content/10.1101/2024.02.05.578874v1.

48. Wen S, Yin A, Furlanello T, Perich M, Miller L, Itti L. Rapid adaptation of brain–computer interfaces to new neuronal ensembles or participants via generative modelling. Nat Biomed Eng. 2021; p. 1–13. https://doi.org/10.1038/s41551-021-00811-z PMID: 34795394

49. Hug F, Avrillon S, Ibáñez J, Farina D. Common synaptic input, synergies and size principle: Control of spinal motor neurons for movement generation. J Physiol. 2023; 601(1):11–20. https://doi.org/10.1113/ JP283698 PMID: 36353890

50. Farris DJ, Lichtwark GA. UltraTrack: Software for semi-automated tracking of muscle fascicles in sequences of B-mode ultrasound images. Comput Methods Programs Biomed. 2016; 128:111–118. https://doi.org/10.1016/j.cmpb.2016.02.016 PMID: 27040836

51. Millard M, Uchida T, Seth A, Delp SL. Flexing computational muscle: modeling and simulation of muscu-lotendon dynamics. J Biomech Eng. 2013; 135(2):021005. https://doi.org/10.1115/1.4023390 PMID: 23445050
52. Zajac FE. How musculotendon architecture and joint geometry affect the capacity of muscles to move and exert force on objects: a review with application to arm and forearm tendon transfer design. J Hand Surg. 1992; 17:5:799-804. https://doi.org/10.1016/0363-5023(92)90445-U PMID: 1401783

53. Caillet AH, Phillips AT, Carty C, Farina D, Modenese L. Hill-type computational models of muscle-tendon actuators: a systematic review. bioRxiv. 2022; p. 2022-10. Available from https://www.biorxiv.org/content/10.1101/2022.10.14.512218v2.

54. Huijing PA, Slawnych MP. Modeling of homogeneous muscle: Is it realistic to consider skeletal muscle as a lumped sarcomere or fiber? In: Biomechanics and Neural Control of Posture and Movement. Springer; 2000. p. 92-99.

55. Yeo SH, Verheul J, Herzog W, Sueda S. Numerical instability of Hill-type muscle models. Journal of the Royal Society Interface. 2023; 20(199):20220430. https://doi.org/10.1098/rsif.2022.0430 PMID: 36722069

56. Perreault EJ, Heckman CJ, Sandercock TG. Hill muscle model errors during movement are greatest within the physiologically relevant range of motor unit firing rates. J Biomech. 2003; 36(2):211-218. https://doi.org/10.1016/S0021-9290(02)00332-9 PMID: 12547358

57. Hussein M, Shebl S, Elnemr R, Elkaranshawy H. A new muscle activation dynamics model, that simulates the calcium kinetics and incorporates the role of store-operated calcium entry channels, to enhance the electromyography-driven hill-type models. J Biomech Eng. 2022; 144(1):011002. https://doi.org/10.1115/1.4051718 PMID: 34251438

58. Song S, Kidziński L, Peng XB, Ong C, Hicks J, Levine S, et al. Deep reinforcement learning for modeling human locomotion control in neuromechanical simulation. J Neuroeng Rehabil. 2021; 18:1-17. https://doi.org/10.1186/s12984-021-00919-y PMID: 34399772

59. Winters JM. An improved muscle-reflex actuator for use in large-scale neuromusculoskeletal models. Ann Biomed Eng. 1995; 23:359-374. https://doi.org/10.1007/BF02584437 PMID: 7486344

60. Elias LA, Watanabe RN, Kohn AF. Spinal mechanisms may provide a combination of intermittent and continuous control of human posture: predictions from a biologically based neuromusculoskeletal model. PLoS Comput Biol. 2014; 10(11):e1003944. https://doi.org/10.1371/journal.pcbi.1003944 PMID: 25393548

61. Proctor J, Holmes P. Reflexes and preflexes: on the role of sensory feedback on rhythmic patterns in insect locomotion. Biol Cybern. 2010; 102:513-531. https://doi.org/10.1007/s00422-010-0383-9 PMID: 20358220

62. Geyer H, Herr H. A muscle-reflex model that encodes principles of legged mechanics produces human walking dynamics and muscle activities. IEEE Trans Rehabil Eng. 2010; 18(3):263-273. PMID: 20378480

63. Wang J, Qin W, Sun L. Terrain adaptive walking of biped neuromuscular virtual human using deep reinforcement learning. IEEE Access. 2019; 7:92465-92475. https://doi.org/10.1109/ACCESS.2019.2927606

64. Gallina A, Disselhorst-Klug C, Farina D, Merletti R, Besomi M, Holobar A, et al. Consensus for experimental design in electromyography (CEDE) project: High-density surface electromyography matrix. J Electromyogr Kinesiol. 2022; 64:102656. https://doi.org/10.1016/j.jelekin.2022.102656 PMID: 35344841

65. Caillet AH, Phillips AT, Modenese L, Farina D. NeuroMechanics: electrophysiological and computational methods to accurately estimate the neural drive to muscles in humans in vivo. J Electromyogr Kinesiol. 2024; p. 102873. https://doi.org/10.1016/j.jelekin.2024.102873 PMID: 38518426