Automated estimation of ankle muscle EMG envelopes and resulting plantar-dorsi flexion torque from 64 garment-embedded electrodes uniformly distributed around the human leg

The design of personalized movement training and rehabilitation pipelines relies on the ability of assessing the activation of individual muscles concurrently with the resulting joint torques exerted during functional movements. Despite advances in motion capturing, force sensing and bio-electrical recording technologies, the esti- mation of muscle activation and resulting force still relies on lengthy experimental and computational procedures that are not clinically viable. This work proposes a wearable technology for the rapid, yet quanti- tative, assessment of musculoskeletal function. It comprises of (1) a soft leg garment sensorized with 64 uniformly distributed electromyography (EMG) electrodes, (2) an algorithm that automatically groups electrodes into seven muscle-specific clusters, and (3) a EMG-driven musculoskeletal model that estimates the resulting force and torque produced about the ankle joint sagittal plane. Our results show the ability of the proposed technology to automatically select a sub-set of muscle-specific electrodes that enabled accurate estimation of muscle excitations and resulting joint torques across a large range of biomechanically diverse movements, un- derlying different excitation patterns, in a group of eight healthy individuals. This may substantially decrease time needed for localization of muscle sites and electrode placement procedures, thereby facilitating applicability of EMG-driven modelling pipelines in standard clinical protocols.


Introduction
Human movement emerges from mechanical forces generated by muscle-tendon units (MTUs) spanning multiple skeletal degrees of freedom (Winter, 2009). The design of effective training and rehabilitation pipelines relies on the ability to assess in vivo musculoskeletal function i.e., the mechanical force generated by individual MTUs and the resulting torque exerted on biological joints during functional movements (Maggioni et al., 2016). In this context, the ability to perform quantitative, yet rapid, muscle-level functional assessment would have a large impact on domains ranging from clinical-decision making to rehabilitation, from sport training, to injury prevention at the home-and the work-place (Lloyd, 2021).
However, despite advances in movement measuring technologies such as motion capture (e.g., camera, inertial measurement unit-based), force sensing (e.g., in-ground force plates, instrumented treadmills, sensorized pressure insoles) and bio-electrical recording (e.g., highdensity electromyography, HD-EMG), the progress towards musculoskeletal functional assessment has been hampered by the lack of techniques that can measure individual MTU activation and resulting force in a rapid and quantitative way (Modenese et al., 2018).
In clinical settings, where rapidity is crucial, musculoskeletal function is often assessed via functional assessment categories (FAC) (Mehrholz et al., 2007). FAC-based assessment is however subjective and approximative. Here repeatability in subjects' assessments across evaluators is not assured. Moreover, FAC-based assessment does not provide information on how individual skeletal muscles contribute to the observed movement.
Conversely, a more quantitative picture of an individual's musculoskeletal function could be achieved using advanced motion capture technologies, and computer-based EMG-driven neuromusculoskeletal (NMS) modeling techniques (Sartori and Sawicki, 2021;Durandau et al., 2018). This could provide indicators of an individual's MTU activation, force and resulting joint dynamics. However, the EMG-dependence in these muscle-driven modeling pipelines implies the manual localisation of multiple muscles and subsequent electrodes' placement. In this, human error as well as differences in individual experience can affect placement accuracy and repeatability (Stegeman and Hermens, 2007). Moreover, the time needed for muscle localization and placement procedures, hampers the applicability of EMG-driven NMS modelling pipelines in standard clinical protocols, where large numbers of patients are assessed within stringent time windows.
This work proposes a technology for the fast and quantitative assessment of musculoskeletal function, at the MTU-level, in the intact human in vivo. This comprises a soft and sensorized leg garment that incorporates 64 uniformly distributed EMG electrodes, an algorithm that automatically groups garment-embedded electrodes into musclespecific clusters, and a EMG-driven NMS model that uses automatically derived muscle-specific EMG linear envelopes to estimate the mechanical force produced by individual MTUs as well as the resulting plantar-dorsi flexion torque (Lloyd and Besier, 2003;Winter, 2009) (see Fig. 1). In this context, a crucial challenge is that of assuring high-fidelity in the automatically derived muscle-specific EMG linear envelopes, something which is further challenged by the textile and stretchable nature of our proposed garment and by the use of dry-electrodes.
The novelty of our work is the combination of garment-embedded electrodes with non-negative matrix factorization (NNMF)-based automated electrode clustering and EMG-driven NMS modelling pipelines.
Several studies used textile-embedded EMG sensors in different domains, such as prosthetic control, rehabilitation and sport applications (Guo et al., 2019). Textile-embedded low-dimensional sets of bipolar electrodes were previously presented in commercial products (e.g. Myotech and Athos), which allowed measuring lower limb activations from muscles of interest during dynamic conditions including running or cycling (Finni et al., 2007;Coleman, 2021). Textile-embedded bipolar electrodes were also used to classify hand and arm gestures (Assad et al., 2013) as well as for the control of a myoelectric prosthetic hand (Alizadeh-Meghrazi et al., 2022). Textile-embedded dry HD-EMG electrodes for the upper limb were used to record the activation from 100 flat knitted dry electrodes during nine isometric contraction tasks of hand and wrist. The recorded EMGs were finally used to myoelectrically control the hand and wrist functions of an active prosthesis (Farina et al., 2010). Textile-embedded dry HD-EMG electrodes were used to monitor the activity of a target muscle during isometric and dynamic tasks for upper and lower limb, respectively (Cerone et al., 2021).
Signal-driven techniques including K-means vector quantization (Staudenmann et al., 2009), watershed image segmentation (Vieira et al., yyyy) as well as principal component analysis (PCA) and NNMF (Gallina et al., 2018), were previously applied on HD EMGs to identify sub-compartments of activation in single skeletal muscles. These approaches were susceptible to over-clusterization and were employed in single-muscle isometric contractions. Although, NNMF was recently utilised to study muscle compartmentalization during dynamic tasks, this was done on single muscles (Huang et al., 2016). Furthermore, all these methodologies were disjointed with subsequent EMG-driven NMS modelling pipelines. Both PCA and NNMF were employed to investigate multi-muscle synergies during dynamic movement. However, only bipolar EMGs were adopted and manually placed on targeted muscles (Ivanenko et al., 2004;Bizzi and Cheung, 2013).
In this work, we propose an NNMF-based EMG clustering methodology to extract multi-muscles sites from 64 equally distributed electrodes around the leg. The proposed approach relies on the a priori knowledge of ankle muscle synergies and modularity during walking at one fixed speed Cappellini et al., 2006). First, we hypothesize that the automatically extracted muscle sites (derived from walking muscle synergies at one speed) provide accurate estimates of muscle-specific EMG envelopes (by direct comparison EMG envelopes from manually placed electrodes) during movements biomechanically different and underlying different muscle synergies, which were not used for the initial clustering. Second, we hypothesize that the automatically extracted muscle-specific EMG envelopes provide the fidelity necessary to translate into dynamically consistent plantar-dorsi flexion torque profiles across the same subject population and movement repertoire.

Experimental procedures
The study was approved by the Natural Sciences and Engineering Sciences Ethics Committee of the University of Twente (reference number 2020.37). Data were recorded from eight healthy subjects (age = 27.6 ± 1.2 years, height = 170.4 ± 6.9 cm, weight = 68.8 ± 11.7 kg) that volunteered after signing an informed consent. Subjects performed nine tasks as detailed in Table 1. All the tasks were performed on an instrumented split-belt treadmill (M-Gait, MotekForce Link, Netherlands) that recorded ground reaction forces (GRFs) at 1024 Hz. EMGs were detected from the right leg using the sensorized 0-series garment ( Fig. 2) described in the Section 1 of the supplementary material. Before applying the garment, the leg was shaved and the skin was cleaned with alcohol. A ground electrode, not embedded in the garment, was attached to a bony area on the right lateral epicondyle or patella. The EMGs were sampled at 2048 Hz from 64 dry floating electrodes using a multi-channel amplifier (REFA, TMSi, The Netherlands). All EMG channels were amplified against the average of all connected floating inputs, i.e., average reference mode, with exception of the ground electrode that was not an active input. It was required to keep patient and amplifier potentials at around the same level. Movement data were recorded (128 Hz) using retroreflective markers placed on the subject's bony landmarks and anatomical segments (Sartori et al., 2014) using eight optical motion capture cameras (Qualisys Oqus, Sweden).

Data processing
Kinematic and kinetic data were low-pass filtered (6 Hz) with a zerolag 2 nd order Butterworth filter. EMGs signals were amplified with a gain of 20 (REFA, TMSi, The Netherlands). The raw monopolar EMG signals were automatically inspected to remove channels with very large voltage fluctuations due to movement artifacts. Channels with voltage > 0.35 mV or with standard deviations > 20 scaled median absolute deviations from the median among channels were zeroed (Gwin et al., 2010;Schlink et al., 2020). The number of zeroed channels across all subjects and tasks ranged between 1 and 9 with an average of 3.1 ± 2.2. A re-referencing processing was applied to the resulting raw EMGs data. The mean across the non-fluctuating channels was subtracted from the original signals. EMGs were then high-pass filtered (20 Hz) using a zerolag 2 nd order Butterworth filter and full-wave rectified. To maximize EMG-to-force relationship, linear envelopes were obtained via low-pass filtering at 6 Hz with a zero-lag 2 nd order Butterworth filter. EMG envelopes were normalized using the maximum values extracted from all performed tasks via a moving median window (30 ms). Zeroed channels were reconstructed via bilinear interpolation (Bovik, 2009) of the four neighboring (top, bottom, right and left) electrode envelopes. If one of the neighboring channels was silent, i.e. showed a flat line signal, the reconstruction was not applied and the channel was kept to zero.

Multi-channel EMG clustering 2.3.1. Non-negative matrix factorization
Standard NNMF (Lee and Seung, 1999) was extended to enforce partbased extraction by incorporating gradient smoothing and independence constraints into the factorized bases. This prioritizes factorization into components that reflect contiguous regions of activation, thereby identifiable by the activity of a whole muscle. Multiplicative rules were adjusted to enforce weightings locality, by including a gradient-based smoothing constraint (Joshi et al., 2010): where A represented the weightings, B the primitives, and G = − ▽ 2 A and β, γ were empirically determined constants to ensure positive values as previously presented (Joshi et al., 2010).

Associate garment channels to functional groups via NNMF
This is based on two sequential stages. In the first-stage NNMF for each subject, NNMF was applied to all 64 EMG envelopes recorded from three cycles of slow walking. Two sets of primitives/weights were extracted to distinguish between dorsi-and plantar-flexion functional channel groups. The primitives were averaged over the three gait cycles. The average primitive that had bigger area-under-the-curve in the interval between the 40% and 60% of the gait cycle was labelled as reference plantar-flexor primitive. The other primitive was labelled as reference dorsi-flexor primitive. Labelled reference averaged primitives were input to the second-stage NNMF.
Second-stage NNMF located seven sets of primitives/weights to be subsequently associated to muscles including tibialis anterior (TA), extensor hallucis longus (EHL) as dorsi-flexors and peroneus longus (PL), brevis (PB), gastrocnemius medialis (GM), lateralis (GL), and  soleus (SOL) as plantar-flexors (see Section 2.3.3) (Gonzalez-Vargas et al., 2015;Ivanenko et al., 2004). However, the second-stage NNMF is employed differently with respect to the first-stage NNMF. In the first stage, the NNMF was used to extract conventional walking muscle synergies which differentiated the 64 electrodes in two functional electrode clusters. Electrodes that activated together were grouped because the underlying muscles were involved in the same functional movement, i.e.. dorsi-or plantar-flexion. In the second-stage, NNMF is applied to decompose the activation of the entire leg, represented by 64 channels, into its single components, i.e.., the muscles. Here, NNMF was used to extract distinct anatomical clusters, i.e.. groups of close electrodes that activate together because they are placed on top of individual muscles. Garment electrodes were distributed into two groups of 32 channels: garment's four upper and lower electrode rows. In light of SENIAM guidelines (Stegeman and Hermens, 2007), we assumed that the upper rows were covering the areas ranging from the distal half of the GM to the tibial tuberosity. The lower rows covered from the end of distal half of the GM to the lateral malleolus ( Fig. 4.C). NNMF was applied respectively to EMG envelopes measured from the 32-upper and 32lower electrodes separately. Across all subjects, four sets of primitives/weights were always extracted from the 32-upper channels, while three sets of primitives/weights were always extracted from the 32lower channels. The extracted primitives were averaged over the three gait cycles and compared, for shape similarity, with reference averaged dorsi-plantar flexion profiles from first-stage NNMF (see Fig. 3.B).
The Pearson correlation coefficient (R) was computed between second-stage and first-stage dorsi-flexor primitives. Across all subjects, second-stage primitives with R > 0.8 were assigned to the dorsi-flexor group. The remaining primitives were initially considered as plantarflexors. R was computed between second-stage and first-stage plantarflexor primitives. All plantar-flexor primitives with R < 0.35 were discarded, i.e., EMG envelopes contributed by discarded second-stage plantar-flexor primitives were reconstructed (by multiplying primitives with related weightings) and subtracted from the original EMG envelopes as input to second-stage NNMF. Subsequently, NNMF was iteratively applied on the updated set of EMG envelopes until all the R values for plantar-flexor similarity were >=0.35 (Fig. 3.B).
As for the dorsi-flexor threshold (R > 0.8), the 0.35-threshold used for plantar-flexor discrimination was chosen to be unique across all subjects. The lowest R-value observed between first-versus second-stage plantar-flexor primitives, was selected and employed for the entire experimentation. For each subject, the lowest threshold was always associated with the PL and ranged between 0.37 (subject 2) and 0.74 (subject 8) with an average of 0.57 ± 0.12. Excluding subjects 2, 5 and 7 (lowest thresholds of 0.37, 0.48, 0.43), all remaining subjects had lowest threshold always > 0.56.

Clusters-to-muscle mapping
Each of the seven second-stage primitives were associated to a set of 64 weightings. Each weighting reflected how much a given electrode was represented by the associated primitive. For each primitive, all channels with a weighting > 70% of the channel set's maximal weighting were kept to form a functional cluster whose EMG envelopes are most represented by the associated primitive (see Fig. 3.C).
The location of each of the seven functional clusters within the garment was computed using the center of gravity (CoG) (Rojas- Martínez et al., 2012). A heuristic rule was defined to map the position of each cluster's CoG to a specific leg muscle's anatomical location. In this context, the most distal plantar-flexion cluster's CoG was always associated with PB. The most distal dorsi-flexion cluster's CoG was always associated with EHL. The most proximal dorsi-flexion cluster's CoG was always associated with TA. The most proximal plantar-flexion cluster's CoG was always associated with PL. For the remaining clusters' CoG, the GL was always associated to the cluster's CoG located in the upper part of the garment and being the closest to the PL cluster's CoG. The GM's CoG was always assumed to be adjacent to the GL's CoG in the medial direction. The SOL was always assumed to be in the lower part of the garment. Finally, for each task, the normalized (Section 2.2) linear envelopes of the selected electrodes from the same cluster were averaged to form muscle-specific EMG envelopes.

Manually selected channels
To validate the automated clustering procedure (Section 2.3.3), two adjacent electrodes belonging to the same column of the grid were manually selected as the closest to each muscle belly following the SENIAM recommendations (Stegeman and Hermens, 2007). This was done for TA, GL, GM, PL and SOL muscles. Afterward, the manual muscle-specific EMG envelopes were obtained by subtracting the rereferenced (see Section 2.2) raw EMG signals of the muscle-specific electrode pair, and by processing the resulting signal as explained in Section 2.2.

Human movement modeling
The open-source software OpenSim (Delp et al., 2007) was used to create subject-specific musculoskeletal geometry models. A generic musculoskeletal geometry model (Delp et al., 1990) was scaled to match each subject's anthropometry using 3-D marker trajectories derived from static standing trials. Marker trajectories and GRFs recorded during dynamic trials were used to solve inverse kinematics (IK) and inverse dynamics (ID), and to derive leg joint angles and reference ankle torques.

EMG-driven NMS modeling
Joint torques were also computed via the open-source CEINMS software toolbox we previously developed (Pizzolato et al., 2015;Sartori et al., 2012). MTU lengths and moment arms were obtained from the subject-specific geometry model as previously presented (Sartori et al., 2012). Muscle-specific EMG envelopes and MTU kinematics were used to drive forward Hill-type muscle models (see Table 2). Resulting MTU force estimates were multiplied with MTU moment arms to estimate net ankle plantar-flexion torque (Sartori et al., 2012).
For each subject, the EMG-driven NMS model was calibrated to identify the value of muscle-tendon force-generating parameters that vary non-linearly across subjects' anthropometries. This was framed as an optimization that minimized the mean squared error between reference and EMG-driven NMS model estimated torques over a range of calibration trials including static standing pose (3s), two cycles of  backward and forward walking (see Table 1). For subject 1 and 8 the static task was not used due to marker occlusion. Two calibrated EMGdriven NMS models were obtained for each subject using muscle-specific EMG envelopes derived from manually and automatically selected electrodes.

Validation procedures
The first test assessed whether first-stage NNMF (Section 2.3.2) could be applied to a set of 64 EMG channels uniformly distributed throughout the leg and extract muscles synergies with a structure resembling that of conventional synergies, i.e., those derived using bipolar electrodes manually located on targeted muscles. Reference dorsi-flexor and plantar-flexor primitives (Section 2.3.2) were averaged across all subjects and gait cycles during slow walking and compared with state-of-art synergies from studies during walking tasks at similar speeds (Ivanenko et al., 2004;Clark et al., 2010;Gonzalez-Vargas et al., 2015). Primitives were normalized with respect to their maximum value. Finally, the coefficient of determination (R 2 ) were computed between the subjects' dorsi-flexor and plantar-flexor average profiles and the ones from the mentioned studies.
The second test assessed the ability of the clustering algorithm (Section 2.3) to select muscle-specific electrodes that were valid for "unseen" and biomechanically different movements, i.e., movements that were not used for the automated clustering, which underlay different muscle synergies and force profile. Across all subjects, muscles and tasks, R 2 and root mean squared error (RMSE) were calculated, along with median and percentiles (25% and 75%), between EMG envelopes derived from automatically and manually selected channels.
The third test assessed the ability to estimate ankle joint torques via subject-specific NMS models when driven by EMG envelopes derived from manually and automatically selected electrodes respectively. The average across all cycles of manually and automatically selected EMGdriven torque estimations were compared with the averaged reference torque using R 2 and normalized RMSE (NRME)
Second test: Fig. 5 displays subject 1's mean profiles for automatically and manually derived EMG envelopes across muscles and motor tasks. Fig. 6 displays the distribution R 2 and RMSE between automatically and manually derived EMG envelopes across all muscles, subjects and tasks. Error metrics are grouped in three categories: metrics from (1) FW1, i.e., the task used for the clusterization, (2) from other walking tasks FW3, FW5, BW1, BW3 (Table 1) and (3) from running and side stepping. Across all categories RMSEs were always < 0.2. R 2 were always > 0.7 for at least 56% of all trials, with the exception of PL, which had R 2 > 0.7 for 28% of trials in category (1) and 26% of all trials in category (2). In category (3) TA had R 2 > 0.7 for 16% of all trials. In category (2) GL had R 2 > 0.7 for 38% of all trials. Supplementary Table 1 reports mean R 2 and RMSEs in details.
Third test: Fig. 7 displays EMG-driven model-based estimates of torques against reference profiles across all subjects, trials and tasks that were not used for model calibration. Fig. 8 shows the distribution of R 2 and NRMSEs across all trials not used for calibration. Motor tasks were divided into two categories: (1) tasks used for the model calibration such as FW and BW, (2) tasks not used for the model calibration including RN and SS. In category (1), the 77% and 80% of NRMSE values for,

Fig. 5.
Manually (E m -blue lines) and automatically (E a -red lines) retrieved normalized EMG envelopes (E). Mean (solid lines) and standard deviation (dashed lines) across all movements cycles are plotted for subject 1, across all movements (rows) and muscles (columns). The red shaded area highlights the EMG envelopes of the task used in the channel selection method. Data are reported over the movement cycle event with 0% being heel-strike. Fig. 6. Distribution of coefficient of determination (R 2 ) and root mean square error (RMSE) between automatically and manually derived EMG envelopes across all subjects. Three categories are reported: (1) task used for channel selection (forward walking at 1 km/h, FW1 -black bars), (2) additional walking tasks, FW3, FW5, back walking at 1 km/h and 3 km/h (FW-BW -dark grey) and (3) tasks including running and side stepping (RN-SS -light grey). Blue lines refer to normalized R 2 and RMSE distributions considering all trials, tasks and subjects. Vertical blue lines show median (solid line) and percentiles at 25% (dashed line) and 75% (dash and dot line). respectively automatically and manually derived EMG-driven torque estimations, were always < 0.1. In category (2), the 69% and 78% of NRMSE values, for automatically and manually derived EMG-driven torque estimations, were always < 0.1.
In terms of R 2 , For category (1), automatically and manually derived EMG-driven torque estimations presented, respectively, 91% and 92% of the trials with R 2 value > 0.7. In category (2) the percentage of trials with R 2 > 0.7 was 66% and 64% for automatically and manually derived EMG-driven torque estimations respectively. Supplementary Table 2 reports detailed similarity metrics.

Discussions
We developed a technology that linked a soft leg garment equipped with 64 uniformly distributed EMG dry electrodes with a signal-driven pipeline for the automatic selection of muscle-specific electrodes and for the model-based estimation of ankle torque during a set of biomechanically different motor tasks. These movements were chosen to cover a broad range of movement conditions (e.g., speeds) and excitation patterns of lower-limb muscles. However, even if the muscle-specific clusters are extracted from a single task, i.e., a slow walking movement, we wanted to show how the obtained clusters can generalize to different movements underlying different muscle synergies.
The proposed technology prevents from having to manually locate leg muscles and manually apply bi-polar electrodes on target locations. Instead, a single garment is directly worn, enabling rapid placement of 64 EMG channels. Moreover, the proposed clustering algorithm required a minimal EMG-dataset (i.e., from three gait cycles) and was based on two sequential EMG envelopes factorization stages and on a single set of R-thresholds that could be used across all subjects and muscles (see second-stage NNMF in Section 2.3.2). This simplifies applicability to different subjects, crucial for translation towards out-of-the-lab settings.
In this context, a low R-threshold (0.35) was favoured for secondstage NNMF plantar-flexor primitive discrimination because of the observed R variability. Across all subjects and second-stage NNMF plantar-flexor primitives, R varied between 0.37 < 0.57 ± 0.12 < 0.99. Variability was linked to the fact that plantar-flexor synergies contributed to the activation of larger sets of muscles, i.e., five muscles in our study, than dorsi-flexor synergies which, though synergistic, presented differences in peak activation timing and shape as previously reported (Whittle, 2014).
The first test (Section 3, Supplementary fig. 1) showed that first-stage NNMF could factorise the garment-recorded 64 EMG envelopes into synergies that were in line with literature (Clark et al., 2010;Gonzalez-Vargas et al., 2015;Ivanenko et al., 2005). This provides evidence that reference plantar-dorsi flexion synergies could be derived from the garment's large sets of uniformly distributed electrodes that did not target a piori chosen muscles. This enabled using a priori knowledge on ankle muscle modularity during slow walking to initially distribute 64 electrodes into two macro groups of ankle antagonist muscles. Moreover, this enabled a second-stage NNMF to further decompose input EMGs into seven groups of synergistic electrodes (Fig. 3).
The second test showed that our proposed clustering used a priori information on muscle modularity during walking and, when applied to motor tasks capturing different muscle synergies, it led to electrophysiologically and dynamically consistent muscle-specific EMG envelopes (Fig. 5). The least favourable similarity metrics were always obtained for the peroneus longus (Fig. 6). This may be linked to the muscle proximity to tibialis anterior and gastrocnemius lateralis. Future work will assess whether a garment with increased electrode density in the correspondence of the leg lateral sites may help improve localisation of small muscles such as the peroneus longus.
Remarkably, the automatically selected channels always provided EMG envelopes with the level of fidelity required for driving forward subject-specific EMG-driven NMS models and estimate EMG-dependent ankle plantar-flexion torques across all conditions (Figs. 7, 8). This was also true for movements not used for the model calibration (Fig. 7). Torque estimates from manually and automatically derived EMG envelopes displayed large similarity (Fig. 8, R 2 > 0.72, NRMSE < 0.21 across all subjects and trials). However, results showed underestimation of reference plantar-flexor peak torque during walking with discrepancies increasing with walking speed. This may be caused by inability of the calibration procedure to find a unique set of model parameters that can simultaneously fit five motor tasks, as well as by the possible delay  (2) tasks not used for model calibration (lower row). Category (1) plots comprised all trials except the one used for calibration. Vertical lines indicate the median (solid line) and the percentiles at 25% (dashed line) and 75% (dash and dot lines).
introduced during the synchronization of motion tracking and EMG systems during the data recording. Future works should systematically investigate the effect of the chosen tasks in the calibration outcome as well as quantify and correct the delay between EMG and motion tracking systems' signals.
During running and forward walking at 5 km/h, torque oscillations were visible in the swing phase. However, these were present when using both manually and automatically derived EMGs. This may be due to movement artefacts both in EMGs and reflective markers located on fatty areas, which may be translated into unwanted oscillations in MTU kinematics and torque. Future work should develop wireless garments and use inertial measurement units located on bony landmarks to minimise artefacts.
Study limitations included the garment being produced in one size only, thereby not assuring similar levels of fit across all subjects. Future work should investigate the ability of the proposed clustering methodology to generalize across different body parts as well as across individuals with neuromuscular injuries (e.g., post-stroke individuals), who underlie abnormal activation patterns (Ting et al., 2015). Additionally, future studies will investigate how the NNMF-based methodology can be employed online, thereby removing the deed for defining the a priori dimensionality assumptions on the number of active muscles during walking (i.e., seven) and the number of conventional walking muscle synergies (i.e., two). Moreover, our re-referencing technique might lead to distortion in the monopolar EMG signals. However, our first test results (Section 3, Supplementary fig. 1) were in line with the literature. Future work will systematically assess the effect of the rereferencing method on monopolar EMG signals against data referenced using a common reference electrode. Furthermore, the proposed technology is not fully portable as it relies on the use of lab-based camera systems and force plates. Future work will employ fully wearable sensors for kinematic and kinetic measurements.

Conclusions
We developed a smart wearable garment equipped with 64 EMG channels uniformly covering the human leg. We developed a pipeline that enabled the automated identification of muscle-specific electrodes and the subsequent estimation of realistic muscle excitations and resulting forces exerted by the leg musculoskeletal system. The approach was validated across eight healthy subjects and across seven biomechanically diverse movements. The proposed technology has the potential to prevent time-consuming muscle manual localization procedures and facilitate the use of advanced EMG-driven modelling pipelines in the clinics as well as in recreational and occupational domains.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.