Shihan Ma

Irene Mendez Guerra

Arnault Hubert Caillet

Jiamin Zhao

Alexander Kenneth Clarke

Kostiantyn Maksymenko

Samuel Deslauriers-Gauthier

Xinjun Sheng

Xiangyang Zhu

Dario Farina
email: d.farina@imperial.ac.uk

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

Neuromechanical studies investigate how the nervous system interacts with the musculoskeletal (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, are mainly limited to static contractions and cannot fully represent the dynamic modulation of EMG signals during volitional movements. Here, we overcome these limitations by presenting NeuroMotion, an open-source simulator that provides a full-spectrum synthesis of EMG signals during voluntary movements. NeuroMotion is comprised of three 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 BioMime, a deep neural network-based EMG generator that receives nonstationary physiological parameter inputs, such as 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. Here we also provide representative applications of NeuroMotion. We first show how simulated 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 regression. Ridge regressors trained on the synthetic dataset were directly used to predict joint angles from experimental data. NeuroMotion is the first full-spectrum EMG generative model to simulate 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 complex neuromechanical system, fast iterate algorithms before collecting experimental data, and validate algorithms that estimate non-measurable parameters in experiments. We expect this full-spectrum model will complement experimental approaches and facilitate neuromechanical research.

October 6, 2023 1/20

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 musculoskeletal model, a neural network-based EMG generator, and an advanced motoneuron model. 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 demonstrate that the synthetic signals generated by NeuroMotion can be used to augment experimental 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 [START_REF] Ting | Neuromechanics of muscle synergies for posture and movement[END_REF][START_REF] Schwartz | Movement: how the brain communicates with the world[END_REF], by studying the movements from the perspective of their neural control [START_REF] Nishikawa | Neuromechanics: an integrative approach for understanding motor control[END_REF]. 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 [START_REF] Nishikawa | Neuromechanics: an integrative approach for understanding motor control[END_REF], rehabilitation engineering [START_REF] Ting | Neuromechanical principles underlying movement modularity and their implications for rehabilitation[END_REF], and human-machine interfaces [START_REF] Farina | The extraction of neural information from the surface EMG for the control of upper-limb prostheses: emerging avenues and challenges[END_REF].

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, including the classical Hodgkin-Huxley model [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]. At a macro scale, several analytical and numerical electromyography (EMG) models have been used to study the electrical outputs produced by skeletal muscles upon activations of populations of neurons [START_REF] Farina | Compensation of the effect of sub-cutaneous tissue layers on surface EMG: a simulation study[END_REF][START_REF] Farina | A surface EMG generation model with multilayer cylindrical description of the volume conductor[END_REF][START_REF] Maksymenko | A myoelectric digital twin for fast and realistic modelling in deep learning[END_REF][START_REF] Mesin | A finite element model for describing the effect of muscle shortening on surface EMG[END_REF]. In parallel, multiple simulation platforms, such as OpenSim [START_REF] Delp | OpenSim: open-source software to create and analyze dynamic simulations of movement[END_REF][START_REF] Seth | OpenSim: Simulating musculoskeletal dynamics and neuromuscular control to study human and animal movement[END_REF], have been widely used to study human body movements during dynamic simulations in biomechanical studies. However, 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 estimating the musculo-tendon dynamics and parameters but does not forwardly generate EMG signals during a movement [START_REF] Delp | OpenSim: open-source software to create and analyze dynamic simulations of movement[END_REF]. Fuglevand's model has been widely used to study muscle force, motoneuron activities, and EMG signals mainly under isometric contractions [START_REF] Fuglevand | Models of recruitment and rate coding organization in motor-unit pools[END_REF]. A recurrent neural network has been recently used as a black box to convert motions to the downsampled and smoothed muscle activities without knowing the internal parameters of the system [START_REF] Schmidt | The concepts of muscle activity generation driven by upper limb kinematics[END_REF]. One primary reason for the lack of an integrated and precise EMG simulator feasible for voluntary movements is the absence of models that link EMG generation to movement biomechanics. Another reason is that current advanced EMG models are not efficient enough to adapt to the non-stationary physiological parameters during a voluntary movement. These two challenges have been addressed by our recently released EMG model, BioMime [START_REF] Ma | Conditional Generative Models for Simulation of EMG During Naturalistic Movements[END_REF], which is a conditional

Core Modules

MSK Model with OpenSim

OpenSim is an open-source software platform that is used for modelling and analysing neuromusculoskeletal systems [START_REF] Delp | OpenSim: open-source software to create and analyze dynamic simulations of movement[END_REF][START_REF] Seth | OpenSim: Simulating musculoskeletal dynamics and neuromuscular control to study human and animal movement[END_REF]. With this platform, 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 [START_REF] Mcfarland | A Musculoskeletal Model of the Hand and Wrist Capable of Simulating Functional Tasks[END_REF] (ARMs in abbreviation) from the large database of opens-source OpenSim models as the October 6, 2023 3/20 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 [START_REF] Modenese | Automatic generation of personalised skeletal models of the lower limb from three-dimensional bone geometries[END_REF][START_REF] Valente | nmsBuilder: Freeware to create subject-specific musculoskeletal models for OpenSim[END_REF] and associated machine learning methods [START_REF] Saxby | Machine learning methods to support personalized neuromusculoskeletal modelling[END_REF].

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 fibre 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 October 6, 2023 4/20 motoneuron pools, according to the linearity properties of the motoneuron pools [START_REF] Farina | Common synaptic input to motor neurons, motor unit synchronization, and force control[END_REF].

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 [START_REF] Ma | Conditional Generative Models for Simulation of EMG During Naturalistic Movements[END_REF]. A schematic of BioMime's training and inference pipeline is shown in Figure 2. Learning from the outputs of its teacher numerical model [START_REF] Maksymenko | A myoelectric digital twin for fast and realistic modelling in deep learning[END_REF], 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 [START_REF] Ma | Conditional Generative Models for Simulation of EMG During Naturalistic Movements[END_REF]).

Input MUAP

Output MUAP Discriminator

Prior distribution 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.

Specified conditions

True / False Generator

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 [START_REF] Ma | Conditional Generative Models for Simulation of EMG During Naturalistic Movements[END_REF]). 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 [START_REF] Maitland | Ultrasound-guided motor unit scanning electromyography[END_REF][START_REF] Hauraix | Force-length relationship modeling of wrist and finger flexor muscles[END_REF].

The location of the innervation zone and the muscle fibre conduction velocity can be estimated, for example by using high-density EMG recordings [START_REF] Zhang | Global innervation zone identification with high-density surface electromyography[END_REF][START_REF] Martinez-Valdes | Early motor unit conduction velocity changes to high-intensity interval training versus continuous training[END_REF]. Some parameters can also be obtained from the 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 October 6, 2023 5/20

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. These properties are sampled from a prior distribution by sampling. Therefore, generating MUAPs by sampling only uses the decoder. Since the seven parameters explain the bulk of variations of the MUAP templates of the specific subject, 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 [START_REF] Fuglevand | Models of recruitment and rate coding organization in motor-unit pools[END_REF] and a cohort of leaky fire-and-integrate (LIF) neuron models adapted from [START_REF] Caillet | Estimation of the firing behaviour of a complete motoneuron pool by combining electromyography signal decomposition and realistic motoneuron modelling[END_REF].

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 [START_REF] Fuglevand | Models of recruitment and rate coding organization in motor-unit pools[END_REF][START_REF] Caillet | Motoneuron-driven computational muscle modelling with motor unit resolution and subject-specific musculoskeletal anatomy[END_REF], and the summation of the motor unit sizes meets the expected total number of fibres in one muscle (Equation 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 Equation 2. Equation 2 was derived from experimental measurements in the literature on human forearm muscles [START_REF] Fuglevand | Models of recruitment and rate coding organization in motor-unit pools[END_REF][START_REF] Romaiguère | Physiological properties of the motor units of the wrist extensor muscles in man[END_REF][START_REF] Riek | Recruitment of motor units in human forearm extensors[END_REF] following the method described in [START_REF] Caillet | Estimation of the firing behaviour of a complete motoneuron pool by combining electromyography signal decomposition and realistic motoneuron modelling[END_REF]. 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 [START_REF] Fuglevand | Models of recruitment and rate coding organization in motor-unit pools[END_REF][START_REF] Feinstein | Morphologic studies of motor units in normal human muscles[END_REF]. 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 [START_REF] Caillet | Mathematical relationships between spinal motoneuron properties[END_REF] 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.

M N size = f tw (j) N k=1 f tw (k) • N f , j ∈ [1, N ] (1)
where 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.

f tw (j) = 0.81 • (18.51 • j N + 104.10 j N 4.83 ), j ∈ [1, N ] (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 firing 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 [START_REF] Fuglevand | Models of recruitment and rate coding organization in motor-unit pools[END_REF]. The recruitment thresholds are exponentially distributed such that few motor units have October 6, 2023 6/20 high threshold and many motor units have low threshold. 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 supported by experimental data [START_REF] Caillet | Estimation of the firing behaviour of a complete motoneuron pool by combining electromyography signal decomposition and realistic motoneuron modelling[END_REF] in the LIF model, the predictions of both models are consistent with the onion skin theory [START_REF] Luca | Relationship between firing rate and recruitment threshold of motoneurons in voluntary isometric contractions[END_REF] and Henneman's size principle of sequential motoneuron recruitment [START_REF] Henneman | Recruitment of motoneurones: the size principle[END_REF].

Toolbox Functions

In

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 mov' (Figure 3a). The joint angles can be obtained from motion capture data (like in [START_REF] Sgambato | High Performance Wearable Ultrasound as a Human Machine Interface for wrist and hand kinematic tracking[END_REF]) or be measured from sensors (as by data glove 1 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 mov' 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. NeuroMotion provides eight default poses, including hand open/grasp, wrist flexion/extension, wrist radial/ulnar deviation, and forearm pronation/supination (Figure 3b). Users can define their customised poses by creating new poses in OpenSim GUI. A movement can be simulated by setting 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, as the third way displayed in Figure 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.

1 # build msk model and simulate a movement

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 [START_REF] Mesin | A finite element model for describing the effect of muscle shortening on surface EMG[END_REF], 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). 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 [START_REF] Zhao | A musculoskeletal model driven by muscle synergy-derived excitations for hand and wrist movements[END_REF].

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 [START_REF] Caillet | Estimation of the firing behaviour of a complete motoneuron pool by combining electromyography signal decomposition and realistic motoneuron modelling[END_REF]. 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.

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 motoneuron pool, the neural commands to the muscles in the

MUAPs across Poses

We

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, which is defined by the mean square error between two MUAPs divided by their averaged 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 Figure 6. MUAPs within the same muscle 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 the ECRB muscles.

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

the joint flexion or extension angle increased. The similarities of MUAPs across joint angles were consistent within the MUAPs of each muscle with small variances. MUAPs in extensors showed higher levels of dissimilarity during wrist extension than during wrist flexion, and the opposite for MUAPs in flexor muscles.

Synthetic EMG signals for data augmentation

The synthetic EMG signals generated by NeuroMotion can be used to augment experimental datasets for myoelectric control. In this case study, we simulated hand and wrist movement following the experiments in [START_REF] Zhao | Consistent control information driven musculoskeletal model for multiday myoelectric control[END_REF]

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 Figure 4 and Figure 5 showed simulations of MUAPs by using the interpolation tools. Physiological parameters 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 [START_REF] Glaser | Motor unit identification from high-density surface electromyograms in repeated dynamic muscle contractions[END_REF][START_REF] Kramberger | On the prediction of motor unit filter changes in blind source separation of high-density surface electromyograms during dynamic muscle contractions[END_REF][START_REF] Oliveira | Neural control of matched motor units during muscle shortening and lengthening at increasing velocities[END_REF]. For example, the shortening of muscle fibres or increasing 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 Figure 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 [START_REF] Glaser | Simulations of high-density surface electromyograms in dynamic muscle contractions[END_REF] or empirically transformed the MUAP waveforms [START_REF] Chen | Adaptive real-time identification of motor unit discharges from non-stationary high-density surface electromyographic signals[END_REF]. We expect that NeuroMotion will meet the demand for providing synthetic EMG signals during dynamic muscle contractions for validating adaptive decomposition algorithms [START_REF] Glaser | Motor unit identification from high-density surface electromyograms in repeated dynamic muscle contractions[END_REF].

We also showed that the generated MUAPs were muscle-specific with low similarity across different muscles (Figure 6). The variations of MUAPs during a movement were consistent within muscle but different across muscles (Figure 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. Test-bed for Cause-Effect Studies and Data Augmentation for

Regression Analysis

With access to all inputs, outputs, and intermediate variables in NeuroMotion, users can study cause-effect relationships between the joint kinematics, neural drives, 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 (Figure 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 (Figure 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 gave an example to show that ridge regressors trained on the synthetic dataset can be directly applied to regress joint angles from experimental signals. Furthermore, we found that the performance of the 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 [START_REF] Maksymenko | A myoelectric digital twin for fast and realistic modelling in deep learning[END_REF][START_REF] Wen | Rapid adaptation of brain-computer interfaces to new neuronal ensembles or participants via generative modelling[END_REF]. Another important observation from the regression example is that the properties of the between the neural input and the joint kinematics, for instance, by adjusting the input drives 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.

First, the current BioMime model in NeuroMotion is used as a replica of a specific volume conductor system but not a universal one. 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. Details of data preparation and model training can be found in [START_REF] Ma | Conditional Generative Models for Simulation of EMG During Naturalistic Movements[END_REF]. Building a universal BioMime model that can be adapted to specific subjects and allows a random sampling of electrical potentials on the surface is also in our future work.

Second, the muscle-related results obtained from MSK simulations are highly sensitive to some MSK parameters [START_REF] Caillet | Hill-type computational models of muscle-tendon actuators: a systematic review[END_REF][START_REF] Navacchia | Prediction of in vivo knee joint loads using a global probabilistic analysis[END_REF][START_REF] Myers | A probabilistic approach to quantify the impact of uncertainty propagation in musculoskeletal simulations[END_REF][START_REF] Valente | Are subject-specific musculoskeletal models robust to the uncertainties in parameter identification?[END_REF], including the muscle-specific maximum isometric force, tendon slack length, optimal fibre length, and the muscle moment arm. These parameters in the ARMs model may deviate from the specific subject, which will affect the estimation of muscle activation. Scaling the ARMs model or building a personalised MSK model can address some of the above issues. The personalised models could be further used to simulate patient's data, such as fatigue and limited mobility in people with spinal cord injuries.

Third, changes in physiological parameters during voluntary movement are estimated by first extracting the changes in muscle fibre length using OpenSim and then approximating the changes in conduction velocity and motor unit depth under physiological assumptions. Currently, the fibre lengths are derived from the geometrical muscle-tendon lengths computed with the ARMs model in OpenSim by assuming the tendons to be rigid and at constant slack length [START_REF] Caillet | Motoneuron-driven computational muscle modelling with motor unit resolution and subject-specific musculoskeletal anatomy[END_REF][START_REF] Millard | Flexing computational muscle: modeling and simulation of musculotendon dynamics[END_REF]. This assumption of stiff tendon [START_REF] Zajac | Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control[END_REF] is acceptable given that the investigated muscles have ratios of tendon slack length to optimal fibre length below 5.0 [START_REF] Zajac | Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control[END_REF]. However, the muscle-tendon interacting dynamics during contraction are not considered under this assumption, which might introduce errors in estimating the changes in muscle fibre length and muscle activation. We expect that the accuracy of NeuroMotion will be improved with the advancement of in vivo measurements of physiological parameters during movements.

Conclusion

We

  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 dynamic contractions. 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 movement. 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 Figure 1

Fig 1 .

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

Fig 2 .

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

1Fig 3 . 2 msk

 32 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 MCP (metacarpophalangeal) and PIP (proximal interphalangeal) 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.

1 # 13 # 17 #

 11317 test motor unit pool model 2 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 . as si g n _ p r o p e r t i e s ( config , True ) Force and twitches 18 mn_pool . init_twitches ( fs )mn_pool . i n i t _ q u i s i s t a t i c _ e f _ m o d e l ()

  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 ContractionsMUAPs 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 (Figure4a) 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 Figure4b. 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 Figure4c, d, and e, respectively.

Fig 4 .

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

Figure 7 Fig 5 .

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

Fig 6 .

 6 Fig 6. Similarity 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.

Fig 7 .

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

RFig 8 .

 8 Fig 8. Representative results of using NeuroMotion's output to train regression algorithms and to augment the experimental dataset. a), The regressors trained on the synthetic dataset could predict wrist and MCP joint angles from the experimental data, with the Pearson correlation coefficient (R) 0.62 and 0.72, respectively. b), When trained on the dataset augmented by the synthetic signals, the ridge regressors achieved better performance compared with the regressors trained only on the experimental data.

  synthetic EMG signals largely depend on neural inputs/spike trains. It still remains a major challenge to decode the real neural inputs during a voluntary movement, and thus impossible to set an exactly accurate neural input in NeuroMotion. Here, we used the normalised EMG signals recorded on each muscle as the neural inputs. The flexibility provided by NeuroMotion allows the investigation of the actual relationship October 6, 2023 15/20

  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' own 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 neuromechanics fields.

  . 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 correspond to the myoelectric activities responsible for the movement. Three modules are incorporated 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 modular simulation. Details of the three key modules are described in Section Core Modules. An overview of the workflow of NeuroMotion is shown in Figure 1.

  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 physiological parameters (Section Estimate Parameter Changes), set common drives to motoneuron 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.

  Set Neural Inputs to Motoneuron PoolsActivations 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.

	October 6, 2023	8/20

  Code Block 2. Example Python codes for defining a motor unit pool model and for setting the input drive to the pool.

	_ , spikes , fr , ipis = mn_pool . g e n e r a t e _ s p i k e _ t r a i n s ( ext )

Author summary