In the lower limb dataset, we recorded data from nineteen young volunteers (14 male and 5 females, with age 26.0 ± 3.0 (MEAN ± SD) (years)) and fourteen aged volunteers (12 male and 2 female, with age70.0 ± 6.0 (years)). The physical activity of the participants was assessed according to the international physical activity questionnaire (IPAQ) (Booth, 2000), and they were then ranked with minimally active (low) or active (moderate or high) scores (Sember et al., 2020). Based on volunteer reports, none of them did have any balance impairment, neurological disorders, muscular injuries, or the intake of medications for body balance.

In the upper limb dataset, we recorded data from five healthy male volunteers with an average age of 21.3 ± 2.8  (years), weight of 71.0 ± 3.4 (kg), and height of 174.3 ± 2.6 (cm).

In both datasets, all subjects gave informed consent to the experimental procedure. The procedures were confirmed with the Declaration of Helsinki and were approved by Politecnico di Torino Research Ethics Committee.

In the lower limb dataset, sEMG signals were recorded in signal differential mode. All signals were amplified by a between-individuals variable factor, ranging from 5,000 to 10,000, to maximize the analog gain without saturation. The force signals were recorded by a piezoelectric force plate (9286AA Kistler, Zurich, Switzerland). Both sEMG and force signals were sampled synchronously at 2048 Hz with a 12-bit analog-to-digital converter (EMG-USB, OTBioelettronica, and LISiN, Politecnico di Torino, Turin, Italy). We used linear array electrodes to record the calf muscles’ activity: sEMG from tibialis anterior and medial and lateral gastrocnemius muscles were detected with three arrays of 16 electrodes, each with 10 mm inter-electrode distance (IED), whereas two arrays of four electrodes each with 10 mm IED were used to sample EMGs medially and laterally from soleus (cf Figure 1) in (dos Anjos et al., 2017). For the gastrocnemius muscle, the most proximal electrode was located 2 cm distal to the popliteal fossa, and the arrays were aligned parallel to the longitudinal axis of each gastrocnemius head. For the tibialis anterior muscle, the array was aligned 1 cm laterally and parallel to the tibial crest, with the most proximal electrode located 2 cm distal to the fibula’s head. The soleus muscle’s lower border of the medial and lateral arrays was positioned 3 cm distal to the medial gastrocnemius myotendinous junction (dos Anjos et al., 2017).

In the upper limb dataset, sEMG signals were recorded from the Biceps Brachii (BB), Brachioradialis (BR), Triceps Brachii lateral (TBL), and medial head (TBM) during isometric voluntary flexions-extensions with the angle of the elbow at 90 degrees. For recording from BB, we used an HD-sEMG array of 64 circular electrodes disposed into five columns and 13 rows with 8 mm IED (one missing corner electrode). For BR, TBL, and TBM recording, we used three linear arrays of 8 electrodes with an IED of 5 mm. The main IZ was located for each muscle before placing the electrode array. Depending on the subject’s anatomical features, we placed the adhesive arrays either proximally or distally from the main IZ location. We used paste to clean the skin before the recording (Meditec-every, Parma, Italy). We recorded data in monopolar mode with the same amplifier used in the lower limb experiments, allowing for the synchronous recording of HD-sEMG and force signals after amplifying the latter by a factor of 100 (Force Amplifier MISO-II, LISiN, Politecnico di Torino, Italy). The force signal was shown to the subjects as real-time feedback (Botter et al., 2011).

In the lower limb experiments, participants were asked to stand upright on a force plate with eyes open and arms alongside the body. Ground reaction forces were measured in two conditions while subjects: i) kept the position of their center of pressure at 65% of the distance between the tip of the calcaneus bone and the tip of the third metatarsal head and; ii) stood at ease. Both trials lasted 60 s, and center-of-pressure visual feedback was provided in the first trial. This task was selected to ensure a somewhat high degree of calf muscle active loading while not threatening stability, particularly for aged individuals. Five-minute intervals were applied between trials, and their order was randomized.

In the upper limb dataset, each of the volunteers was asked to do three maximum voluntary isometric flexions and extension contractions for 5 s before starting the experiment, and the highest was chosen as the Maximum Voluntary Contraction (MVC). After a few pieces of training for each subject, they were asked to perform a series of flexion-extension force ramps for 25 s. Each series consisted of four isometric ramps from different percentages (30%, 50%, and 70%) extension MVC to flexion MVC and back. The single differential signals were used to reduce muscle crosstalk (Botter et al., 2011; Jafari et al., 2014).

The lower limb dataset was recorded in single differential mode, and the signals were digitally bandpass filtered with a fourth-order Butterworth filter in the 15–350 frequency band. Given that we were interested in estimating the temporal force profile, regardless of its absolute value, the center of pressure position in the sagittal plane was taken as the ankle force signal (Morasso and Schieppati, 1999). In the upper limb dataset, the monopolar HD-sEMG signals were digitally bandpass filtered (20–450 Hz with a fourth-order Butterworth filter), and the force signal was low pass filtered (cut off at 1 Hz with a fourth-order non-casual Butterworth filter).

This section briefly introduces the least square (LS) problem and its solution based on calculus.

Suppose that we have linear equations A x = b , where matrix A has n rows (i.e., equations) and m columns (i.e., unknowns), x is the unknown vector with m rows, and b is an n − v e c t o r . These equations have a solution if b is a linear combination of columns of A .

For most cases, we need to find an x for r = A x − b , that minimizes the residuals. We choose x that can minimize the norm of the residual, A x − b . Minimizing the norm of the residual is similar to minimizing its square. The problem of finding m − v e c t o r x ^ that minimize A x − b 2 is called a least squares problem. x ^ = argmin x A x − b 2 (1)

This residual is an affine function. Affine functions are considered linear functions, and they are convex. To solve this problem, we must calculate the gradient. ∂ f ∂ x i x ^ = 0 , i = 1 , … , m (2)

Or we can consider it in vector notation ∇ f x ^ = 0 . (3)

The gradient in the matrix form is described below ∇ f x = 2 A T A x − b = 0 (4)

The solution can be presented as the following: A T A x ^ = A T b (5)

This equation is called a normal equation, and with the assumption that columns of A are linearly independent, the solution is calculated as below: x ^ = A T A − 1 A T b (6)

If the matrix A is not full rank, the normal equation must be solved with Moore-Penrose pseudo-inverse with the help of the singular value decomposition (SVD) (Strang et al., 1993). When the SVD is calculated, the reciprocal of the non-zero items of the middle matrix is calculated while keeping the zero items, and multiplying SVD matrices from right to left is then used as the pseudo-inverse matrix.

The proposed method had the minimum adj. R ² of 80.13% in the entire test sets of the upper and lower-limb datasets (Tables 1, 2), showing a proper goodness-of-fit (Figures 4–10). Since the number of samples was much more than the number of parameters of the models, adj. R ² and R ² indices were very similar (Eq. 22; Tables 1, 2). Based on the Bland Altman plots of Figure 4), the residual signal of the proposed method has better homogenous variance in the upper-limb dataset compared with the lower-limb dataset. The regression line of the residual signal of the upper-limb dataset was y = 0.006848–0.06733×x, while it was y = −0.0003940 + 0.00007351×x for the lower-limb dataset, where parameters x and y are the measured and residual force signals.

The regression line of the scatter plots of Figures 6, 7 were y = 0.008371 + 0.9175 ×x (R ² = 85.6%), and y = 0.007631 + 0.9612×x (R ² = 96.1%), respectively, showing proper fitness of the predicted force signal (x) vs. the measured signal (y) in the upper and lower-limb datasets. In the entire datasets, the maximum running time of the proposed algorithm was 40.8 us and 17.5 us on the training and testing datasets, respectively. Thus, it is suitable for real-time applications.

The cut-off frequency of 2.0 Hz was used in our study for envelope estimation. The sEMG amplitude during isometric and quasi-static contractions was shown to have a frequency below 5–10 Hz, and the optimal cut-off frequencies between 2–3 Hz were provided in the literature (Staudenmann et al., 2010; Ranaldi et al., 2022). Such a cut-off frequency compromises sEMG dynamics, muscle force difference, and the time lag between those signals (Staudenmann et al., 2010). Different cut-off frequencies could affect the smoothness resulting in correlations among estimated sEMG amplitude and any target signal, especially in non-slowly varying signals (Ranaldi et al., 2022). Moreover, non-causal digital filters were used in our study. In a real-time application, however, causal filters must be used. Since both filters of the sEMG and muscle force were non-causal, such a modification should not be problematic.

Squared or absolute sEMG signals could be used for envelope detection. While the former is optimal for normal distribution, the latter is used when the distribution of the sEMG signals is more centrally peaked (a.k.a., Laplacian) (Hogan and Mann, 1980; Clancy and Hogan, 1999; Clancy et al., 2004). The distribution of our signal epochs was mostly Laplacian, as identified by Goodness-of-fit Test Statistics for the Laplace Distribution (Puig and Stephens, 2000). Moreover, the variation of the absolute signal is lower than that of the squared signal, resulting in a superior signal-to-noise ratio (Clancy and Hogan, 1999). Thus, the absolute signal (a.k.a., full-wave rectifier) was used in our study for envelope detection.

The proposed algorithm is a convex representation of models 2 and 4, in which the natural logarithm of the sEMG envelopes is used. The logarithmic transformation is primarily used in statistics to reduce the skewness of the data (West, 2022), resulting in a more normally distributed dataset, which is preferred in clinical practice (Feng et al., 2014). Such a transformation is also helpful to use LS, a convex optimization method, rather than stochastic optimization methods (e.g., PSO), which are time-consuming and get stuck in local minima (Beck, 2014) (Tables 1, 2). Moreover, the logarithmic transformation, a member of the Box-Cox family of transformations, could improve the performance of the regression models in system identification (Keene, 1995; Ljung, 1999).

The proposed model could classify active and minimally active subjects of the lower limb dataset, mainly based on the importance of the tibialis anterior muscle during sway protocol. The algorithm’s sensitivity is 100% to identify active subjects, while its specificity is 78% to detect minimally-active subjects. Tibialis anterior muscle was shown as an essential muscle during quiet standing sway in healthy and also subjects with neurodegenerative diseases (such as Parkinson’s disease) (Vette et al., 2017) (Warnica et al., 2014). The soleus is one of the principal sources of proprioceptive information during standing (Di Giulio et al., 2009), whose exercises might help to reduce the risk of falling in the elderly (Lee and Yoo, 2017). Figure 10 shows that the tibialis anterior (RL) or medial soleus (LL) are mainly active during sway protocol. They are the antagonist’s muscles and have important roles in postural control (Nagai et al., 2012).

In the proposed algorithm, the matrix A in Eq. 6 is estimated using the envelope of the sEMG signals (Eq. 8). We used the stable solution of Eq. 6 (i.e., A\b) instead of calculating the inverse matrix (A^TA)^−1, which is highly affected by the condition number of the matrix (Pintelon and Schoukens, 2012; Beck, 2014). The condition number of the matrix A equals the square root of the condition number of the matrix (A^TA). The condition number of the positive semidefinite matrix (A^TA) equals the maximum eigenvalue of (A^TA) divided by the minimum eigenvalue of (A^TA) (Beck, 2014).

We further calculated the correlation between the adj. R ² and the condition number of the matrix A (Eqs 6, 8). The average condition number of matrix A was 763 ± 464. No significant correlation was found between the goodness-of-fit and the condition number of the matrix in the entire dataset (r = −0.263; p-value = 110). Our method not only had proper goodness-of-fit but also did not have a significant bias in the entire dataset (p-value = 0.610).

In our model (Eq. 7), we also estimated the intercept point (w₀). If it is subtracted from the measured force signal (f(t)), the mean value of the measured signal is removed, and the number of estimated parameters is decreased. It is aligned with the recommendation of the system identification guidelines to remove the measured signals detrend before system identification.

The same LS model was used in our method for the upper and lower-limb datasets, showing the generalization capability of the proposed algorithm. However, our method’s goodness-of-fit and residual signal analysis were better on the upper-limb dataset than on the lower-limbs dataset (Tables 1, 2; Figure 4). It could be that the number of involved muscles in the upper-limb dataset was less than that of the lower-limb dataset. Moreover, unlike black-box models such as ANN or SVM (Mokri et al., 2022), our method is a grey-box model, in which the model interpretation is possible (Figures 9, 10). The activity of each muscle could be provided for each active muscle during recording to provide insights into load-sharing problems (Rojas-Martínez et al., 2019). Such information is also helpful in prosthesis control to identify major active muscles when a limited number of electrodes are used in practice.

Traditional EMG amplitude processing methods were used for comparison with our proposed model. Since analyzed methods had statistically significantly different results using the rm-ANOVA (F (4,567,728) = 969.856; p-value < 0.001), the pairwise comparison was performed, and to reduce the Type I error, the Bonferroni correction was used in our study and adj. p-value was reported. Overall, the proposed method significantly outperformed the other methods (adj. p-value < 0.05), in terms of the goodness-of-fit measures. The comparison between ANN, OLS, and SVM methods and the proposed method on a recording from the upper-limb dataset was further provided in Figure 5. ANN does not provide proper fitness among the analyzed models, as shown in the Bland-Altman plot (Figure 1). It could be because the ANN method is a black-box method, and they are not usually acceptable in clinical applications.

The OLS did not have a good Bland Altman plot either (Figure 2), showing that linear weights were unsuitable for the EMG-force problem. However, OLS (a.k.a. Multiple Linear Regression) has been widely used in the literature. It was shown in the literature that the EMG-force relationship (i.e., the sEMG envelope, which is smoothed rectified EMG by itself, and muscle force) needs not be linear (Fuglevand et al., 1993; Hof, 1997), especially in broad force range. The muscle force has different biomechanical components, including the sigmoid shape between excitation and muscle force (Fuglevand et al., 1993), which could be approximated using the exponential forms.

One of the limitations of our study is that we did not use a random force trajectory. It was shown that such a force profile could provide rich information about muscle excitation (Wang et al., 2021), which is the focus of our future work. Using proper signal processing methods, it is possible to identify not only the frequency spectrum of the trajectory force signal to stimulate different frequencies equally, but we can also provide training and test sets without a medium-to-high degree of similarity to reduce such a bias. However, the average correlation between the epochs of estimation and test sets of the measured force signals was −0.0046 ± 0.6812 and −0.0048 ± 0.5194 in the upper and lower limbs datasets, respectively, showing low-to-medium correlation of the measured force signal in estimation and test sets in our study.

Also, the analyzed methods did not consider the electromechanical delay (Lacourpaille et al., 2013). Although properly incorporating the delay parameter could improve the goodness-of-fit of the models (Johns et al., 2016), it increases the complexity of the methods (Corcos et al., 1992; Isermann and Münchhof, 2011) that might not be suitable for real-time applications. Having implemented classical signal conditioning methods on the sEMG signals, a delay is introduced as a solution to this problem, although not a systematic method. An accurate and efficient estimation of this parameter will also be a focus of our future activities. Moreover, in our study, a representative sEMG channel was selected for each muscle. Alternatively, principal component analysis (PCA) could improve the force estimation from HDsEMG signals (Staudenmann et al., 2006).

Moreover, we did not compare our method with sEMG superimposition, which uses motor control information rather than traditional EMG amplitude processing, which has been recently used in muscle force estimation (Savc et al., 2018). In principle, it is possible to use parts of the sEMG decomposition algorithms (such as, (Mohebian et al., 2019) to estimate muscle activity index as to estimate muscle force. However, the running time and the required number of electrodes of sEMG decomposition-based methods are problematic in real-time force prediction compared with fast and reliable EMG amplitude-based methods, such as our proposed method.

It was shown in the literature that amplitude cancellation distorts the spectrum of the rectified sEMG signal (Dideriksen and Farina, 2019). Moreover, many factors influence the relation between EMG amplitude and force, including the amount of crosstalk from nearby muscles, the number and the discharge rate of recruited motor units, Skin-electrode contact, Interelectrode distance, Electrode size and shape, the orientation of the recording array compared with the muscle fibers, Length of the fibers, the shape of the volume conductor, and thickness and inhomogeneities of the subcutaneous tissue layers. Moreover, the same control strategy may generate signals with different amplitude trends based on the motor unit locations (Farina et al., 2004). Moreover, optimal smoothing and envelope detection depends on the load-sharing scenario (Staudenmann et al., 2010). Thus, a specific sEMG amplitude-force relation cannot have general validity and could have a subject-by-subject and muscle-by-muscle relationship. It was also shown in the literature that the relationship between force and absolute sEMG is linear for small muscles with narrow motor unit recruitment force ranges, while it is non-linear for larger muscles with wide motor unit recruitment force ranges (such as proximal leg or arm muscles) (Zhou and Rymer, 2004).

However, our proposed method had an average goodness-of-fit (R ²) of 96.77 ± 1.67 (%) and 91.08 ± 6.84 (%) for the upper and lower-limb datasets, respectively. The subjects-by-subject variations are acceptable, as SD values are less than the uncertainty on the mean values (i.e., MEAN divided by the number of subjects).

The primary elbow flexors are the biceps brachii, brachialis, and brachioradialis, while the triceps brachii is the primary extensor. Moreover, anconeus could contribute to elbow extension (Day, 2009). The brachialis muscle is a fusiform muscle located deep to the biceps brachii. Since the sEMG of such deep muscles cannot be collected with surface electrodes (e.g., brachialis muscle) remains a substantial limitation and a significant cause of the error (Botter et al., 2011). However, it could be assumed that the approximation error term includes the activity of the deep muscle if such an activity is not correlated with other agonist or antagonist muscles since the residuals must be uncorrelated with the predictors in regression analysis. Triceps brachii is the antagonist, and brachialis is a synergist with biceps brachii. Thus, the condition applies to elbow flexion and extension.

For the lower limb muscles, we did not miss any major muscles concurring with the production of plantar flexion torque. Moreover, during quiet standing, dorsiflexors are well-acknowledged to be silent (Di Giulio et al., 2009). Therefore, the sEMG signals we collected from gastrocnemii and soleus represent the net excitation commanding plantar flexion. However, our sEMG signals were likely not sensitive to changes in excitation of the proximal soleus region, which is covered by the gastrocnemius heads. Of more critical concern, though, is collecting sEMG sensitive to mediolateral differences in soleus excitation, which have been shown to change during quiet and perturbed standing (Cohen et al., 2020; Cohen et al., 2021). With our experimental protocol, sEMG signals were most sensitive and specific to the target muscles, minimizing both Types I and II errors (Vieira and Botter, 2021).

The required sample size of the experiment was calculated based on the expected goodness-of-fit (e.g., R2), the number of predictors, and the statistical power in our regression analysis (Alasuutari et al., 2008; Stanley, 2022). For the upper limb dataset, with the expected R ² of 0.964 (Jafari et al., 2014), five variables, and a statistical power of 80%, the minimum sample size was 8. For the lower limb dataset, however, with the expected R ² of 0.960 (Mokri et al., 2022), 11 variables, and a statistical power of 80%, the minimum sample size was 14. In our study, five subjects participated in the upper limb study, while 33 subjects were enrolled in the lower limb study. Thus, the reliability of generalization of the lower-limb experiment is expected, while this is not fully guaranteed in the upper-limb dataset.