Elbow Joint Stiffness Functional Scales Based on Hill’s Muscle Model and Genetic Optimization

Abstract

The ultimate goal of rehabilitation engineering is to provide objective assessment tools for the level of injury and/or the degree of neurorehabilitation recovery based on a combination of different sensing technologies that enable the monitoring of relevant measurable variables, as well as the assessment of non-measurable variables (such as muscle effort/force and joint mechanical stiffness). This paper aims to present a feasibility study for a general assessment methodology for subject-specific non-measurable elbow model parameter prediction and elbow joint stiffness estimation. Ten participants without sensorimotor disorders performed a modified “Reach and retrieve” task of the Wolf Motor Function Test while electromyography (EMG) data of an antagonistic muscle pair (the triceps brachii long head and biceps brachii long head muscle) and elbow angle were simultaneously acquired. A complete list of the Hill’s muscle model and passive joint structure model parameters was generated using a genetic algorithm (GA) on the acquired training dataset with a maximum deviation of 6.1% of the full elbow angle range values during the modified task 8 of the Wolf Motor Function Test, and it was also verified using two experimental test scenarios (a task tempo variation scenario and a load variation scenario with a maximum deviation of 8.1%). The recursive least square (RLS) algorithm was used to estimate elbow joint stiffness (Stiffness) based on the estimated joint torque and the estimated elbow angle. Finally, novel Stiffness scales (general patterns) for upper limb functional assessment in the two performed test scenarios were proposed. The stiffness scales showed an exponentially increasing trend with increasing movement tempo, as well as with increasing weights. The obtained general Stiffness patterns from the group of participants without sensorimotor disorders could significantly contribute to the further monitoring of motor recovery in patients with sensorimotor disorders.

1. Introduction

Functional assessment is a basic and indispensable component of neurorehabilitation that identifies sensorimotor deficits in patients, determines priorities, and plans therapeutic protocols, i.e., it monitors the effects of therapy. Clinical scales such as the Fugl-Meyer Assessment (FMA) and the Wolf Motor Function Test (WMFT) are frequently used to assess motor functionality in patients after stroke [1,2]. However, clinical scales are not sensitive enough to capture the quality of sensory and motor performance. To improve our understanding of the mechanisms that drive motor recovery, we need to distinguish between ‘true neurological repair’ (i.e., restitution), in which neurological impairments are restored towards normal, and behavioral compensation strategies [3]. Technological progress and development have enabled the production of various portable and built-in sensors that accurately record kinetic and kinematic parameters, as well as additional parameters capable of objective assessment of the sensorimotor function of the upper extremities [4]. New devices and technologies that contain different sensor systems for measuring trajectories, joint angles, forces, and electromyography (EMG) signals, as well as systems that provide estimations of non-measurable variables (such as muscle effort/force and joint mechanical stiffness or impedance) are of great use for precise kinetic and kinematic analysis of motor behavior. Kinematic data can be obtained during performance of a specific functional task or with specially designed non-functional assays, and depending on the analysis, it is possible to determine whether a given movement is compensatory or becoming more similar to a normal movement. Recovery trials need to consider serially applied kinematic/kinetic measurements alongside clinical assessments to distinguish between restitution and compensation [5].

Among non-measurable variables, the human ability to adequately generate an interaction force is recognized as one of the basic parameters for the analysis of human movement dynamics and its application in neurorehabilitation. Pathological conditions, such as post-stroke, lead to a reduction in joint stiffness and force due to the improper activation of antagonist muscle pairs. While several studies have investigated how to detect the rotational stiffness of lower-limb joints during locomotion, exploiting the cyclicity of gait or implementing musculoskeletal models driven by muscle activity, less attention has been given to the real-time estimation of upper limb stiffness modulation during non-cyclic movements [6].

Basic approaches have been proposed to estimate joint stiffness and force during human arm movement [7]. Stiffness estimation for short-range human arm movement with small variations in the interaction force can be performed using linear models. On the other hand, stiffness estimation for long-distance hand movements with large variations in the forces of interaction with the environment involves nonlinear models using musculoskeletal modeling [8,9].

A common method for estimating joint stiffness during small movements of the human arm using a linear model was proposed by Perreault [10,11]. This approach is based on applying stochastic perturbation to the human hand using a robotic manipulator and measuring the restoring forces imposed by the hand movement. Ajoudani et al. extended the research of Perreault et al. [10,11] to include muscle activation using EMG signals as additional parameters to assess endpoint stiffness. They concluded that endpoint stiffness of the human arm is related to the muscle co-contraction index and arm configuration. In [12], the Jacobian matrix of the human arm and the co-contraction index of the dominant antagonist muscle pair define the endpoint stiffness. Changing the activation of the dominant muscles directly affects the volume of the stiffness ellipsoid, while changing the position of the arm affects the direction of its axis. The identified endpoint stiffness model of the human arm in real time is used in the tele-impedance of the robot. In a continuation of this research, Ajoudani et al. [13] presented the application of new models of the musculoskeletal system of the human arm based on the muscular Jacobian and the synergistic contribution of muscle activities (the long heads of the biceps brachii and triceps brachii) to estimate the endpoint stiffness. The disadvantage of this approach is that it does not include the joint moment within the muscular Jacobian matrix. The evaluation of the complete joint stiffness matrix is incomplete because the endpoint stiffness from the end effector space needs to be transferred to the larger joint space. To identify the actual total joint stiffness of the human hand, they developed a new technique based on nullspace stiffness [14]. It is a complex two-stage identification procedure where, in the first stage, shoulder–hand rotational perturbations are applied with wrist fixation, and in the second stage, translational and rotational perturbations are applied without a wrist brace. On the other hand, Stanev et al. [15] analyzed the problem of redundancy at the muscle level. They proposed a method to estimate the endpoint and joint stiffness of the redundant musculoskeletal model for any arm movement. Since the CNS uses different strategies to coordinate muscles during arm movement, they proposed a nullspace solution for muscle redundancy and calculation of the joint stiffness according to defined tasks and muscle constraints. The presented results are based on simulation of the arm and leg. In the first part of the analysis, a simplified arm model with three degrees of freedom and nine muscles without joint nonlinearity was used, and in the second part of the analysis, the existing OpenSim model was used without any improvements. In [16], a new model was constructed based on the main axes of the stiffness ellipsoid and their lengths, i.e., the Eigen decomposition of the stiffness matrix, based on the arm configuration. Eigen decomposition was much easier compared with Jacobians, especially the muscle Jacobian. They defined four parameters that are closely related to muscle strength and skeletal dimensions. In this way, the stiffness calculation was simplified and allowed for personalized matching of an individual’s physical interaction characteristics. An average model error of about 20% existed due to the neglect of the influence of external forces and their assumption on muscle synergy. This approach is not suitable for other applications that are more sensitive to variation in the stiffness value. Most studies are based on a perturbation experiment, where joint stiffness is estimated from a linear model using small displacement and small variations in muscle activity. This approach is not suitable for estimating joint stiffness in real-world scenarios. Since the human arm behaves nonlinearly on the large scale of hand movement and interaction force represented in real life, nonlinear modeling of the human hand is required [17].

Nonlinear models of the human arm require modeling of the musculoskeletal system, which is predominantly based on Hill’s muscle model and the passive joint structure model [18,19]. This model uses the muscle activation of the antagonistic muscle pair as the input parameter, whereby the muscle activation is estimated using EMG signals of antagonistic muscles during targeted movement. The values of the Hill’s muscle model and passive joint structure model parameters of upper extremities have been shown through several studies [20,21,22,23,24,25,26,27,28], but none of these publications provide a full list of range values of the parameters.

The research performed in this paper takes into account the nonlinearity of human muscles and considers a stiffness estimation method based on Hill’s muscle model and the passive joint structure model, resulting in a complete list of verified model parameters during the elbow flexion/extension task, as well as elbow stiffness curves (so-called functional scales) for different task tempos and loads. It is worth noting that our approach provides continuous estimation of stiffness in time, which is highly beneficial for human–robot interaction tasks [8,29] and for planning the stiffness trajectory of an artificial variable-stiffness elbow joint that is used as a prosthesis [30]. This paper provides a feasibility study for a general assessment methodology that overcomes one-size-fits-all challenges through the application of a genetic algorithm (GA) for subject-specific elbow model parameter prediction. We used an approach based on a commonly accepted joint/muscle physical model [8,18,19,20,21,22,23,24,25,26,27,28], where we derived a subject-specific GA-based model, and then, estimated joint stiffness based on experiments using wearable sensors. This approach gives us the ability to predict and measure stiffness and other non-measurable quantities based on the model, even when measurements are not directly available. The GA-based methodology was applied to a population without sensorimotor disorders during a “Reach and retrieve” task of the Wolf Motor Function Test—WMFT (modified to include elbow extension in addition to elbow flexion) to examine the feasibility of the proposed objective functional assessment, considering that the same methodology could be applied to patients with sensorimotor disorders and upper extremity weakness in the future. The main contributions of the paper are as follows: (1) identification of the range of the Hill’s muscle model and passive joint structure model parameters using an elbow flexion/extension task in a group of ten healthy volunteers without sensorimotor disorders based on GA, providing of a full list of the range values of the parameters; (2) verification of the identified model parameters; (3) continuous-time elbow joint stiffness estimation based on verified subject-specific model parameters; (4) the definition of novel scales for assessment based on elbow joint stiffness estimation in scenarios of different task tempos and elbow loads.

2. Materials and Methods

2.1. Participants

Ten healthy volunteers without sensorimotor disorders participated in the study: 5 males (age 33.6 ± 4.27) and 5 females (age 28 ± 9.35). The participants verbally confirmed that they had no previous injuries to the elbow joint, nor did they currently feel any discomfort. Before the experiment, the selected body segment parameters (BSP) of each participant were measured using the Hanavan model of the human body [31] (Table 1).

Based on BSP parameters, we estimated, for each participant, the mass (m) of the forearm–hand segment, the distance from the elbow joint axis to the center of the forearm mass (r_cm), and the moment of inertia (M) of the forearm-hand segment. The estimated parameters were used as input parameters for Hill’s muscle model and the passive joint structure model (see Section 2.4.1 and Section 2.4.2).

2.2. Experiment Description

At the beginning of the experiment, the participant was sitting on a chair at a table with their arm extended in front of them. The arm was supported by the elbow on the table so that only the forearm was allowed to move in the horizontal plane. The wrist joint was immobilized to avoid its independent movement. This position was taken as the initial position of the arm (zero angles in the elbow, see Figure 1). The monitor for the feedback display of the movement tempo was positioned in front of the participant.

The experimental protocol of the performed study was based on the “Reach and retrieve” task of WMFT (task no. 8, where the participant attempts to pull the weight across the table by flexing the elbow) [32]. In the performed study, the “Reach and retrieve” task was modified, with additional extension of the elbow with the same load as in the flexion.

All participants were subjected to the same experimental conditions. Each participant was instructed to move their hand while holding the weight from the starting position to a position that corresponded to the maximum flexion of the user’s elbow (the range of elbow motion was approximately from 0° to 150°), i.e., touching the chest with the hand. The far-right position of the “moving dot” on the feedback monitor corresponded to an outstretched human arm (the initial position), while the far-left position corresponded to maximum flexion at the elbow. The participants were instructed to avoid pushing the weight against the table and lifting the weight. Friction between the weights and the table was ignored due to the rigid, smooth finish and different materials. The weight was made of stainless steel and the table had a veneer top.

The overall performed task, the modified task 8 of WMFT (WMFT8_m), included: (1) attempting to push the weight across the table by flexing the elbow and (2) attempting to push the weight across the table by extending the elbow.

During the task activity, the Trigno^TM Avanti (Delsys, Natick, MA, USA) system was used for the data acquisition of: (1) surface EMG activity of two antagonistic muscle pairs of the upper arm (triceps brachii long head and biceps brachii long head muscle) and (2) the angle in the elbow. Trigno Avanti Duo EMG units (Delsys, Natick, MA, USA) were used for EMG recordings and a goniometer (Biometrics, Ladysmith, VA, USA) with Trigno Avanti Goniometer Adapter Sensor (Delsys, Natick, MA, USA) was used for angle measurement. The positions of all sensors are shown in Figure 1. EMGworks Acquisition software (Delsys, Natick, MA, USA) was used for simultaneous data acquisition from all sensors. The sampling rate for data acquisition was approximately 2000 Hz.

The overall experimental scenario included the following phases:

• Maximum voluntary contraction (MVC) data acquisition for the triceps brachii long head and biceps brachii long head muscle against static resistance [33,34,35,36]. The recorded values of the MVC were used in all experiments for normalization of the recorded EMG signals (see Section 2.3).

• Training phase. In the beginning, the participant’s arm was in the initial position (Figure 1A). The participant repeated the WMFT8_m rotation (while pulling the 0.5 kg weight) of the elbow joint following the tempo of the “moving dot” shown on the feedback display. The tempo of the “moving dot” was set to 15 beats per minute (bpm), 30 bpm, 45 bpm, and 60 bpm and it was changed for each 10 WMFT8_m movement, respectively. This phase was used for estimating the parameters of Hill’s model (see Section 2.4.1) by GA, as described in Section 2.4.3.

• Task frequency variation experiment (E1). In the beginning, the participant’s arm was in the initial position (Figure 1A). The participant repeated the WMFT8_m rotation (while pulling the 0.5 kg weight) of the elbow joint following the tempo of the “moving dot” shown on the feedback display. The tempo of the “moving dot” was set to 15 bpm, 30 bmp, 45 bmp, and 60 bmp. The participant repeated the WMFT8_m movement for each tempo 10 times in three series. The resting pause between series was 1 min. The resting pause between different tempos was 5 min.

• Load variation experiment (E2). In the beginning, the participant’s arm was in the initial position (Figure 1A). The participant repeated the WMFT8_m rotation of the elbow joint following the tempo of the “moving dot” shown on the feedback display. The tempo of the “moving dot” was set to 30 bmp. The participant repeated the same WMFT8_m movement with different “pulling” weights: 0.25, 0.5, 0.75, and 1 kg. The participant repeated the WMFT8_m movement for each weight 10 times in three series. The resting pause between series was 1 min. The resting pause between different weights was 5 min.

The experimental design was approved by the ethical review board of the University of Belgrade—School of Electrical Engineering. Participants were well-informed about the noninvasive protocol, and they signed informed consent forms.

2.3. EMG Processing and Muscle Activation Dynamics

The acquired EMG and MVC signals from antagonistic pairs were processed to achieve muscle activation. First, the EMG and MVC signals were filtered using the Butterworth bandpass filter (10–500 Hz, order = 5) and the Butterworth Notch filter (48–52 Hz, order = 5). The filtered EMG and MVC signals were rectified, giving the Processed EMG signal (Figure 2, left) and Processed MVC signal. The EMG and MVC envelopes were extracted from the Processed EMG and Processed MVC signals. The EMG envelope was normalized by the maximum value of the MVC envelope, giving the neural activation $\mathrm{u}\left(\mathrm{t}\right)$ (Figure 2). Muscle activation $\mathrm{a}\left(\mathrm{t}\right)$ was calculated using the model with the lowest data point optimization time [23], as suggested by Manal et Buchanan [37], given in Equation (1):

(1)$\mathrm{a}\left(\mathrm{t}\right)=\frac{{\mathrm{e}}^{\mathrm{Au}\left(\mathrm{t}\right)}-1}{{\mathrm{e}}^{\mathrm{A}}-1}$

where A = −20 was taken to achieve the range of muscle activation function [0,1].

An example of the specific measurement of the Processed EMG signal and its envelope data for E1 and E2 belonging to participant 9 is presented in Appendix A.

2.4. Algorithm for Elbow Joint Stiffness Estimation

Figure 3 presents a block diagram for the overall process of elbow joint stiffness estimation based on Hill’s muscle model, the passive joint structure model, GA, and the Recursive Least-Squares (RLS) algorithm. The muscle activation function of the biceps brachii long head muscle (${\mathrm{a}}_{\mathrm{AG}}$(t)) and triceps brachii long head (${\mathrm{a}}_{\mathrm{ANT}}$(t)) are the inputs to Hill’s muscle model. The passive joint structure estimates passive torque (τ_p), damping torque (τ_d), and gravity torque (τ_g) based on the estimated joint angle (q) and velocity ($\dot{\mathrm{q}}$). Total joint torque (τ_J) is defined as the sum of torques coming from the antagonistic muscle pair (τ_tAG, τ_tANT) and passive joint structures (τ_p, τ_d, τ_g). GA adjusts the parameters of Hill’s muscle model for the muscle pair and passive joint structures by taking into account the difference between the measured joint angle (q_goni, using a goniometer) and the estimated joint angle (q). The RLS algorithm was used to estimate the joint stiffness (Stiffness) based on the estimated joint torque τ_J and the estimated joint angle q. The following sections will provide more details on all assessment modules.

2.4.1. Hill’s Muscle Model

A typical and widely used representation of the muscle is Hill’s muscle model [38] (Figure 4). The model contains: (1) a contractile element (CE) (a force generator imposed by neural activation u(t)), (2) a parallel element (PE) that has viscoelastic constants that oppose the movement, and (3) a series element (SE) that contributes to tendon elasticity. The total muscle mass is represented as the moment of inertia of the forearm–hand segment M. Neural activation u(t) in the muscle contractile element will shorten the fiber length of the muscle (l_m), while the coupling force F will increase the length of the muscle.

The force generated by CE (F_CE) is given in Equations (2)–(4):

(2)${\mathrm{F}}_{\mathrm{CE}}={\mathrm{F}}_{\max }\mathrm{a}\left(\mathrm{t}\right)\mathrm{f}\left({\mathrm{l}}_{\mathrm{m}}\right)\mathrm{g}\left({\mathrm{v}}_{\mathrm{m}}\right)$

(3)$\mathrm{g}\left({\mathrm{v}}_{\mathrm{m}}\right)=\frac{1-{\mathrm{v}}_{\mathrm{m}}/{\mathrm{v}}_{\max }}{1+\frac{1}{{\mathrm{a}}_{\mathrm{f}}}{\mathrm{v}}_{\mathrm{m}}/{\mathrm{v}}_{\max }}$

(4)$\mathrm{f}\left({\mathrm{l}}_{\mathrm{m}}\right)=\exp \left(-{\left(\frac{\left(\frac{{\mathrm{l}}_{\mathrm{m}}}{{\mathrm{l}}_{\mathrm{mopt}}}\right)-1}{\mathrm{w}}\right)}^2\right)$

where ${\mathrm{F}}_{\max }$ is the maximum isometric force (the maximum force in the static state at MVC), $\mathrm{a}\left(\mathrm{t}\right)$ is the muscle activation function, ${\mathrm{v}}_{\mathrm{m}}$ is the actual muscle velocity, ${\mathrm{v}}_{\max }$ is the maximum unloaded shortening speed, ${\mathrm{a}}_{\mathrm{f}}$ is the Hill’s model parameter, ${\mathrm{l}}_{\mathrm{m}}$ is the actual muscle fiber length, ${\mathrm{l}}_{\mathrm{mopt}}$ is the optimal muscle fiber length, and $\mathrm{w}$ is the width of the force–length relationship.

The total force ${\mathrm{F}}_{\mathrm{PE}}$ of the PE element is the sum of the forces delivered by the stiffness (F_p) and damping (F_d) force (Equation (5)).

(5)${\mathrm{F}}_{\mathrm{PE}}={\mathrm{F}}_{\mathrm{P}}+{\mathrm{F}}_{\mathrm{d}}$

The muscle can act as a spring and has a nonlinear characteristic. The spring nature of the muscle is modeled as stiffness with its force F_p, given in Equations (6) and (7):

(6)${\mathrm{F}}_{\mathrm{P}}\left({\mathrm{l}}_{\mathrm{m}}\right)=\{\begin{array}{cc}0,& {\mathrm{l}}_{\mathrm{m}}<{\mathrm{l}}_{\mathrm{ms}}\\ \frac{{\mathrm{k}}_{\mathrm{ml}}}{{\mathrm{k}}_{\mathrm{me}}}\left(\exp \left({\mathrm{k}}_{\mathrm{me}}\left({\mathrm{l}}_{\mathrm{m}}-{\mathrm{l}}_{\mathrm{ms}}\right)\right)-1\right),& {\mathrm{l}}_{\mathrm{ms}}\le {\mathrm{l}}_{\mathrm{m}}\le {\mathrm{l}}_{\mathrm{mc},}\\ {\mathrm{F}}_{\mathrm{mc}}+{\mathrm{k}}_{\mathrm{m}}\left({\mathrm{l}}_{\mathrm{m}}-{\mathrm{l}}_{\mathrm{mc}}\right),& {\mathrm{l}}_{\mathrm{m}}>{\mathrm{l}}_{\mathrm{mc}}\end{array}$

(7)${\mathrm{F}}_{\mathrm{mc}}=\frac{{\mathrm{k}}_{\mathrm{ml}}}{{\mathrm{k}}_{\mathrm{me}}}\left(\exp \left({\mathrm{k}}_{\mathrm{me}}\left({\mathrm{l}}_{\mathrm{mc}}-{\mathrm{l}}_{\mathrm{ms}}\right)\right)-1\right)$

where ${\mathrm{l}}_{\mathrm{m}}$ is the actual muscle fiber length, ${\mathrm{l}}_{\mathrm{ms}}$ is the passive muscle slack length, ${\mathrm{l}}_{\mathrm{mc}}$ is the length at which the spring becomes linear and is equal to the optimal muscle length ${\mathrm{l}}_{\mathrm{mopt}}$, ${\mathrm{k}}_{\mathrm{m}}$ and ${\mathrm{k}}_{\mathrm{ml}}$ are parallel spring constants, and ${\mathrm{k}}_{\mathrm{me}}$ is an exponential shape parameter. The PE also contains a damping element for which the damping force F_d is given in Equation (8):

(8)${\mathrm{F}}_{\mathrm{d}}\left({\mathrm{v}}_{\mathrm{m}}\right)={\mathrm{B}}_{\mathrm{m}}{\mathrm{v}}_{\mathrm{m}}$

where ${\mathrm{B}}_{\mathrm{m}}$ is the parallel damping constant and ${\mathrm{v}}_{\mathrm{m}}$ is the change in muscle speed over time.

The tendon connects muscle and bone and is modeled as SE. The force in the muscle tendon F_SE is given in Equations (9) and (10):

(9)${\mathrm{F}}_{\mathrm{SE}}\left({\mathrm{l}}_{\mathrm{t}}\right)=\{\begin{array}{cc}0,& {\mathrm{l}}_{\mathrm{t}}<{\mathrm{l}}_{\mathrm{ts}},\\ \frac{{\mathrm{k}}_{\mathrm{tl}}}{{\mathrm{k}}_{\mathrm{te}}}\left(\exp \left({\mathrm{k}}_{\mathrm{te}}\left({\mathrm{l}}_{\mathrm{t}}-{\mathrm{l}}_{\mathrm{ts}}\right)\right)-1\right),& {\mathrm{l}}_{\mathrm{ts}}\le {\mathrm{l}}_{\mathrm{t}}\le {\mathrm{l}}_{\mathrm{tc}}\\ {\mathrm{F}}_{\mathrm{tc}}+{\mathrm{k}}_{\mathrm{t}}\left({\mathrm{l}}_{\mathrm{t}}-{\mathrm{l}}_{\mathrm{tc}}\right),& {\mathrm{l}}_{\mathrm{t}}<{\mathrm{l}}_{\mathrm{tc}}\end{array}$

(10)${\mathrm{F}}_{\mathrm{tc}}=\frac{{\mathrm{k}}_{\mathrm{tl}}}{{\mathrm{k}}_{\mathrm{te}}}\left(\exp \left({\mathrm{k}}_{\mathrm{te}}\left({\mathrm{l}}_{\mathrm{tc}}-{\mathrm{l}}_{\mathrm{ts}}\right)\right)-1\right)$

where ${\mathrm{k}}_{\mathrm{t}}$ and ${\mathrm{k}}_{\mathrm{tl}}$ represent tendon spring constants, ${\mathrm{k}}_{\mathrm{te}}$ is the exponential shape parameter,${\mathrm{l}}_{\mathrm{t}}$ is the tendon length,${\mathrm{l}}_{\mathrm{ts}}$ is the tendon slack length, and ${\mathrm{l}}_{\mathrm{tc}}$ is the length at which the tendon becomes linear.

Having in mind muscle representation given in Equations (2)–(10), the actual muscle fiber length ${\mathrm{l}}_{\mathrm{m}}$, muscle speed ${\mathrm{v}}_{\mathrm{m}}$, and tendon length ${\mathrm{l}}_{\mathrm{t}}$ are given in Equations (11)–(13):

(11)${\mathrm{v}}_{\mathrm{m}}=\frac{1}{{\mathrm{M}}_{\mathrm{m}}}\left({\mathrm{F}}_{\mathrm{SE}}-{\mathrm{F}}_{\mathrm{CE}}-{\mathrm{F}}_{\mathrm{PE}}\right)$

(12)${\mathrm{l}}_{\mathrm{m}}={{\displaystyle \int }}_0^{\mathrm{t}}{\mathrm{v}}_{\mathrm{m}}\mathrm{dt}$

(13)${\mathrm{l}}_{\mathrm{t}}={\mathrm{l}}_{\mathrm{ms}}+{\mathrm{l}}_{\mathrm{ts}}+{\mathrm{r}}_{\mathrm{p}}{\mathrm{q}}_0-{\mathrm{l}}_{\mathrm{m}}$

where ${\mathrm{M}}_{\mathrm{m}}$ is the muscle mass, ${\mathrm{r}}_{\mathrm{p}}$ is the radius of the elbow joint calculated as a half value of the BSP parameter P10 from Table 1, and ${\mathrm{q}}_0$ is the elbow joint angle during the full extension (${\mathrm{q}}_0=0\mathrm{deg}$).

The joint torque ${\mathsf{\tau }}_{\mathrm{t}}$ generated from the muscle and transmitted through the tendon ${\mathrm{F}}_{\mathrm{SE}}$ via ${\mathrm{r}}_{\mathrm{p}}$ is given in Equation (14).

(14)${\mathsf{\tau }}_{\mathrm{t}}={\mathrm{F}}_{\mathrm{SE}}{\mathrm{r}}_{\mathrm{p}}$

2.4.2. Passive Joint Structures

Passive joint structure implies stability in the musculoskeletal model [39,40]. The passive joint torque ${\mathsf{\tau }}_{\mathrm{p}}$, the damping torque ${\mathsf{\tau }}_{\mathrm{d}}$, and the torque originating from gravity ${\mathsf{\tau }}_{\mathrm{G}}$ are given in t Equations (15)–(17):

(15)${\mathsf{\tau }}_{\mathrm{p}}\left(\mathrm{q}\right)={\mathrm{k}}_1\exp \left(-{\mathrm{k}}_2\left(\mathrm{q}-{\mathrm{q}}_2\right)\right)-{\mathrm{k}}_4\exp \left(-{\mathrm{k}}_5\left({\mathrm{q}}_1-\mathrm{q}\right)\right)$

(16)${\mathsf{\tau }}_{\mathrm{d}}\left(\dot{\mathsf{q}}\right)=-{\mathrm{B}}_{\mathrm{p}}\dot{\mathsf{q}}$

(17)${\mathsf{\tau }}_{\mathrm{G}}={\mathrm{mgr}}_{\mathrm{cm}}\sin \left(\frac{\mathsf{\pi }}{2}-\mathrm{q}\right)$

where $\mathrm{m}$ is the mass of the forearm–hand segment, ${\mathrm{B}}_{\mathrm{p}}$ is the joint damping constant, $\mathrm{g}$ is the gravitational constant, ${\mathrm{r}}_{\mathrm{cm}}$ is the distance between the joint axis and the center of mass of the forearm–hand, ${\mathrm{q}}_1$ and ${\mathrm{q}}_2$ represent joint angles where the torque begins to increase,${\mathrm{k}}_1$ and ${\mathrm{k}}_4$ are joint torque constants, and ${\mathrm{k}}_2$ are ${\mathrm{k}}_5$ are dimensionless parameters. The parameters ${\mathrm{k}}_1,{\mathrm{k}}_2$, ${\mathrm{k}}_4$, and ${\mathrm{k}}_5$ define the shape of the torque joint curve given in Equation (15), which fits the torque joint curve experimentally defined and represented in [40].

The total joint torque ${\mathsf{\tau }}_{\mathrm{J}}$ is the sum of the torques originating from the passive joint structure and the torque induced by the antagonist’s muscle ${\mathsf{\tau }}_{\mathrm{tAN}\mathrm{T}}$ subtracted by the torque produced by the agonist’s muscle ${\mathsf{\tau }}_{\mathrm{tAG}}$ (Equation (18)).

(18)${\mathsf{\tau }}_{\mathrm{J}}={\mathsf{\tau }}_{\mathrm{p}}+{\mathsf{\tau }}_{\mathrm{d}}+{\mathsf{\tau }}_{\mathrm{G}}+{\mathsf{\tau }}_{\mathrm{tAN}\mathrm{T}}-{\mathsf{\tau }}_{\mathrm{tAG}}$

According to Newton’s Second Law for rotation, Equation (16), the joint angular velocity $\dot{\mathsf{q}}$ and the joint angle $\mathrm{q}$ can be obtained through integration and double integration of the angular acceleration $\ddot{\mathrm{q}}$, satisfying the initial static condition of the joints ${\dot{\mathrm{q}}}_0=0$ and ${\mathrm{q}}_0=0$ (Equations (19)–(21)):

(19)${\mathsf{\tau }}_{\mathrm{J}}=\mathrm{M}\ddot{\mathrm{q}}$

(20)$\dot{\mathsf{q}}={{\displaystyle \int }}_0^{\mathrm{t}}\ddot{\mathrm{q}}\mathrm{dt}+{\dot{\mathsf{q}}}_0$

(21)$\mathrm{q}={{\displaystyle \int }}_0^{\mathrm{t}}\dot{\mathsf{q}}\mathrm{dt}+{\mathrm{q}}_0$

where M is the moment of inertia of the forearm–hand segment.

2.4.3. GA for Adjustment of Parameters for Hill’s Muscle Model and Passive Joint Structure Model

GA is defined as a multi-objective optimization approach that aims to calculate the optimal combination of Hill’s and the passive joint structure model parameters (defined in Section 2.4.1 and Section 2.4.2) while considering the difference between the measured joint angle (q_goni) and the estimated joint angle (q). Thus, one combination of coefficients in GA notation is one individual. Each individual consists of chromosomes and represents a unique solution within the solution space. For GA, in this paper, each chromosome represents one parameter. Accordingly, the initial values of these parameters and their ranges are defined and encoded into GA using the actual value of the chromosome representation. The population of the designed GA is represented by a set of 20 individuals. By testing GA with different numbers of individuals in one population, it was concluded that 20 individuals are a good compromise between the complexity of GA, its execution time, and optimization quality. To cover the complete solution space, a uniform distribution of individuals was used to generate the initial population. The fitness function defined by Equation (22) was evaluated for each individual per generation:

(22)${\mathrm{F}}_{\mathrm{fitness}}={{\displaystyle \int }}_0^{{\mathrm{t}}_{\mathrm{end}}}\left({\mathrm{q}}_{\mathrm{goni}}-\mathrm{q}\right)\mathrm{dt}\approx \frac{{\mathrm{t}}_{\mathrm{end}}}{2\mathrm{N}}{{\displaystyle \sum }}_{\mathrm{n}=1}^{\mathrm{N}}\left({\mathrm{q}}_{\mathrm{goni}}\left({\mathrm{t}}_{\mathrm{n}}\right)-\mathrm{q}\left({\mathrm{t}}_{\mathrm{n}}\right)\right)$

where q_goni is the joint angle measured using a goniometer, q is the estimated joint angle according to Equation (21), ${\mathrm{t}}_{\mathrm{end}}$ is the motion duration, and N is the number of samples during the motion.

The value of the fitness function was used to determine the probability that each individual would be selected for the next generation during the selection process. Roulette-wheel selection was used, which selects the solutions with the lowest values of fitness function but also gives a chance to some weaker solutions to survive the selection process. A weaker solution may include some components that may prove useful after the crossover process. The roulette-wheel selection process was chosen to prevent convergence to a local minimum and to try to find a globally optimal solution. The crossover rate was chosen to be equal to 0.7, which means that 70% of new (child) individuals were created using the crossover operator. Other children were defined using the mutation operator or were elite individuals. Two individuals with the lowest values of the fitness function in the current generation (elite individuals) were selected and directly included in the next generation to accelerate the convergence of the GA to the optimal solution. To create new combinations of parameters (individuals) that differed from the existing ones, the mutation operator was used. Mutation provides genetic diversity, allows the GA to search a wider space, and can take individuals out of a local minimum and move them towards a global one. The mutation operator added a Gaussian distributed random value to a chromosome gene. If it fell outside the lower and upper limits for the gene, the new value of the gene was truncated. The imposed condition for stopping the GA was heuristically obtained since the change in the best value of the fitness function F was not greater than ${10}^{-6}$ for the previous 20 generations.

A PC platform with an Intel(R) Core(TM) i7-9700F CPU at 3.00 GHz and 64.0 GB RAM was used for the optimization procedure. The software used for the Hill’s and passive joint structure models was the Matlab2022 (license: Matlab Campus-Wide 2022) package Simulink. The Hill’s and passive joint structure model parameters were estimated using genetic algorithm functions ga within the Matlab Global Optimization toolbox. The ga function contained: an objective (fitness) function that is given in Equation (22); the number of variables set to 24 (the number of Hill’s and passive joint structure model parameters); ‘PopulationSize’ set to 20 individuals; ‘SelectionFcn’ set to ‘selectionroulette’; ‘EliteCount’ set to 2; and ‘CrossoverFraction’ set to 0.7. The values of the Hill’s and passive joint structure model parameters given in papers [19,20,21,22,23,24,25,26,27,28] were used to define the upper and lower bounds of the parameters. Half the value of the given parameters shown in the papers [19,20,21,22,23,24,25,26,27,28] was taken for the lower bound, while twice the value of the given parameters was taken for the upper bound. The CPU operating time was ~30 h for each participant and for each task.

2.4.4. Verification of the Training GA Outputs

The measured joint angle (q_goni) in the Training phase (described in Section 2.2) was used for verification of the GA estimation parameters from the Hill’s and passive joint structure models. The overall accuracy of the estimated parameters obtained by GA was obtained using the estimated joint angle (q) from the passive joint structure model for the values of the estimated parameters in the Training phase. The measured joint angle (q_goni) from the Training phase and the estimated joint angle (q) were compared using Dynamic Time Wrapping (DTW) [41]. DTW gives the sum of the Euclidian distance between corresponding points of estimated joint angle (q) and the measured joint angle (q_goni). The Training GA output verification is represented by the obtained DTW value divided by the number of samples.

2.4.5. Recursive Least Square Algorithm

Analogously to the approach presented in [42,43], the RLS algorithm was applied to the estimated joint torque ${\mathsf{\tau }}_{\mathrm{J}}$ and joint angle q to estimate the elbow joint stiffness. The total joint torque ${\mathsf{\tau }}_{\mathrm{J}}$ was expressed as in Equation (23):

(23)${\mathsf{\tau }}_{\mathrm{J}}=\mathsf{\Phi }\mathsf{\Pi }$

where $\mathsf{\Phi }=\left(1,\mathrm{q},{\mathrm{q}}^2,\dots ,{\mathrm{q}}^{\mathrm{n}}\right)$ is the $\mathrm{n}$ order torque regressor, which is the function of the joint angle, and $\mathsf{\Pi }={\left({\mathsf{\pi }}_1,{\mathsf{\pi }}_2,\dots {\mathsf{\pi }}_{\mathrm{n}}\right)}^{\mathrm{T}}$ is a vector of the state parameters to be estimated. The order $\mathrm{n}=3$ was chosen to capture the main features of the RLS mathematical model while the estimation was simultaneously denoised. The estimation of the state parameter $\mathsf{\Pi }$ is proposed in [44] and given by Equations (24)–(28):

(24)$ϵ\left[\mathrm{k}\right]={\mathsf{\tau }}_{\mathrm{j}}\left[\mathrm{k}\right]-\mathsf{\Phi }\left[\mathrm{k}\right]\mathsf{\Pi }\left[\mathrm{k}\right]$

(25)$\mathsf{\rho }\left[\mathrm{k}\right]={\mathsf{\Phi }}^{\mathrm{T}}\left[\mathrm{k}\right]\mathrm{P}\left[\mathrm{k}-1\right]\mathsf{\Phi }\left[\mathrm{k}\right]$

(26)$\mathrm{G}\left[\mathrm{k}\right]={\left(1+\mathsf{\rho }\left[\mathrm{k}\right]\right)}^{-1}\mathrm{P}\left[\mathrm{k}-1\right]\mathsf{\Phi }\left[\mathrm{k}\right]$

(27)$\mathsf{\Pi }\left[\mathrm{k}\right]=\mathsf{\Pi }\left[\mathrm{k}-1\right]+\mathrm{G}\left[\mathrm{k}\right]ϵ\left[\mathrm{k}\right]$

(28)$\mathrm{P}\left[\mathrm{k}\right]=\mathrm{P}\left[\mathrm{k}-1\right]-\mathrm{G}\left[\mathrm{k}\right]{\mathsf{\Phi }}^{\mathrm{T}}\left[\mathrm{k}\right]\mathrm{P}\left[\mathrm{k}-1\right]$

where $ϵ\left[\mathrm{k}\right]$ is the error of the measured (${\mathsf{\tau }}_{\mathrm{j}}\left[\mathrm{k}\right]$) and estimated torque values ($\mathsf{\Phi }\left[\mathrm{k}\right]\mathsf{\Pi }\left[\mathrm{k}\right]$) in the kth time sample, $\mathrm{P}\left[\mathrm{k}\right]$ is the covariance matrix, and $\mathrm{G}\left[\mathrm{k}\right]$ is the gain of the RLS algorithm. The initial values of $\mathsf{\Pi }$ and $\mathrm{P}$ are given as $\mathsf{\Pi }\left[0\right]={\left[000\right]}^{\mathrm{T}}$ and $\mathrm{P}\left[0\right]=0.01{\mathrm{I}}_{3\mathrm{x}3}$, where $\mathrm{I}$ is an identity matrix. The elbow joint stiffness (Stiffness) in each time sample $\mathrm{k}$ was calculated using Equation (29).

(29)$Stiffness\left[\mathrm{k}\right]=-\mathsf{\Phi }{\left[\mathrm{k}\right]}^{\mathrm{T}}\mathsf{\Pi }\left[\mathrm{k}\right]$

2.4.6. Stiffness Functional Scale

The data acquired in E1 and E2 were used for designing novel functional scales. First, the parameters of the Hill’s and passive joint structure models estimated in the Training phase were applied to the acquired data in E1 and E2. Second, Dynamic Time Wrapping (DTW) [41] was performed on a measured joint angle (q_goni) and the estimated joint angle (q) to quantify the estimation error. Finally, fitting curves (so-called functional scales) for different task frequencies (Stiffness (tempo) in E1) and different weights (Stiffness (weight) in E2) were defined.

3. Results

3.1. Identification of Hill’s and Passive Joint Structure Model Parameters

Table 2 presents the range (maximal and minimal values) of the parameter values of the Hill’s and passive joint structure models for elbow joints obtained using the GA algorithm. The muscle activation function of the biceps brachii long head muscle (a_AG(t)) triceps brachii long head (a_ANT(t)) and the measured elbow joint angle (q_goni) from the Training phase were used as inputs in the Hill’s and passive joint structure models. The range values (minimum and maximum values) of the Hill’s and passive joint structure models’ parameters for the elbow joint for each parameter were calculated as the boundary values (minimum and maximum) of all values estimated for all participants for each parameter.

3.2. Verification of Hill’s Model and Passive Joint Structure Parameter Estimation

Verification of the accuracy of the Hill’s model and passive joint structure parameter estimation was performed using the DTW algorithm (defined in Section 2.4.4) and is presented in Table 3. The values in Table 3 represent the average deviation (in degrees) of the measured (q_goni) and estimated (q) angles in the elbow joint (DTW divided by the number of samples) for 10 participants.

Figure 5 shows an example of the typical pattern for WMFT8_m movement. The mean value (bold line) and standard deviation (shaded area) of these signals were obtained by overlapping all WMFT8_m movements and their normalization with the time: (A) measured (q_goni) and estimated (q) elbow joint angles; (B) estimated stiffness value (Stiffness); and (C) estimated total joint torque value (${\mathsf{\tau }}_{\mathrm{J}}$).

The comparative analysis of stiffness accuracy in relation to the change in model parameter values is presented in Table 4. The values of the model parameters were modified for ±10%, ±20%, and ±30% regarding the parameter values estimated by GA. The accuracy of the stiffness was calculated using the DTW algorithm, comparing the estimated stiffness values using GA and stiffness calculated from modified parameters. The results show that stiffness sensitivity increases with the increase or decrease in each of the parameters. The stiffness is more sensitive to changes in the k₁, k₂, k₄, q₂, l_tc, l_ts, F_max, l_mopt, l_ms, and v_max parameters (DWT ≥ 0.1) (see Table 4), while for the others, the change in stiffness is negligible. It can be observed that in the DWT value for l_tc, l_ts is high, because the algorithm for stiffness diverged after parameter modification.

3.3. Elbow Joint Stiffness Estimation and Functional Scales

The comparative analysis of the stiffness value through different tempos in task E1 is presented in Figure 6. The comparative analysis of the mean value of the stiffness signal for participant 9 for each tempo (15 bpm, 30 bpm, 45 bpm, and 60 bpm) is presented in Figure 6A. The stiffness mean values were calculated by overlapping all WMFT8_m movements and their normalization with the time. Using the boxplots, the stiffness behavior for each tempo for all 10 participants is shown in Figure 6B. Each boxplot represents the median and the interquartile range of the area under the stiffness signal (Figure 5B) for all repetitions of each participant.

A stiffness fitting curve (a so-called functional scale) for different task tempos $Stiffness\left(\mathrm{tempo}\right)$ is given in Equation (30) and shown in Figure 7.

(30)$Stiffness\left(\mathrm{tempo}\right)=0.001975· {\mathrm{tempo}}^{1.709}$

The comparative analysis of the stiffness value through different loads in task E2 is presented in Figure 8. The comparative analysis of the mean value of the stiffness signal for participant 9 for each weight (0, 0.25, 0.5, 0.75, and 1 kg) is presented in Figure 8A. The stiffness mean values were obtained by overlapping all WMFT8_m movements and their normalization with the time. Using the boxplots, the stiffness behavior for each load for all 10 participants is shown Figure 8B. Each boxplot represents the median and the interquartile range of the area under the stiffness signal (Figure 5B) for all repetitions of each participant.

Stiffness fitting curves (so-called functional scales) for different “pulling” weights $Stiffness\left(\mathrm{weight}\right)$ are given in Equation (31) and shown in Figure 9.

(31)$Stiffness\left(\mathrm{weight}\right)=0.39· exp\left(1.014· \mathrm{weight}\right).$

4. Discussion

This paper presents a novel methodology for elbow joint stiffness estimation in people without sensorimotor disorders during different movement scenarios regarding movement tempo and the upper limb load. The overall research goal of the paper is to define a Stiffness functional scale (general Stiffness curve pattern) that is useful for the further neurorehabilitation follow-up of patients.

Considering that the musculoskeletal system of the upper limb is a nonlinear model, Hill’s muscle model and the passive joint structure model are convenient for use in such research [18,19]. One of the challenges when using Hill’s muscle model and passive joint structure is the estimation of the model parameters, which are the physical parameters of the muscles and tendons themselves, as well as the parameters that describe the characteristics of their behavior during movement. To the best of our knowledge, no report includes a complete list of parameter values for the antagonistic pair (the triceps brachii long head and biceps brachii long head muscle) and passive joint structure. Parameter values that describe the lengths of the above-mentioned muscles (l_mopt, l_ts, l_ms) are presented in [20,21,22,23,24,25,26,27,28]; additionally, F_max values are given in [20,21,23,24,27,28], ${\mathrm{M}}_{\mathrm{m}}$ values are given in [22,27], ${\mathrm{a}}_{\mathrm{f}}$ values are presented in [22], and ${\mathrm{v}}_{\max }$ values are given in [22,28]. In this paper, we proposed an experimental task corresponding to the clinical Wolf test, which aims to assess human/patient capabilities. The task comprises the repetition of a movement with different velocities and loads. To the best of our knowledge, there has been no similar scenario performed in the literature, and this fact makes comparison with previous research difficult. Additionally, it is worth remarking that exact equivalence among parameters is impossible due to the high variability in the human body. However, good matching is observed for the parameters l_mopt, l_ts, l_ms, F_max, ${\mathrm{M}}_{\mathrm{m}}$, ${\mathrm{a}}_{\mathrm{f}}$, and ${\mathrm{v}}_{\max }$, as presented in Table 5.

Cavallaro et al. [20] and Zhao et al. [45] used GA in their work to estimate a subset of Hill’s model parameters. Other parameters of the Hill model, as well as all parameters of the passive joint structure for the elbow joint, are not described in the literature, so one of the major contributions of this paper is the presentation of the complete estimation of all parameters of the Hill’s muscle model and the passive joint structure model using GA. The range values (minimum and maximum values) of all the parameters estimated in the Training phase of ten participants are presented in Table 2. Considering that during the experiment, the tempo of the WMFT8m movement was variable (see Section 2.2), different dynamics of the elbow joint movement were taken into account during parameter estimation. The quality of parameter estimation using GA was tested using the DTW algorithm by comparing acquired and estimated values in the elbow joint. The values of the DTW analysis show that the maximum deviation is 9.16 deg, which corresponds to 6.1% of the full elbow angle range values during one WMFT8m movement (see Table 3). The parameters estimated from the Training phase were applied in experiments E1 and E2 for all participants to estimate stiffness in the elbow joint. A maximum deviation of 12.16 deg, which corresponds to 8.1% of the full elbow angle range values during one WMFT8m movement, was obtained for E1 and E2 (see Table 3). Additionally, a visual representation of the pattern matching of the estimated and acquired elbow joint angle is presented in Figure 5A. Stiffness in the elbow joint for participant 9 in task E1 and task E2 is presented in Figure 6A and Figure 8A, respectively. The maximum values of the estimated stiffness for task E1 are in the range [1.75 Nm/rad–7 Nm/rad], while those for task E2 are in the range [2.2 Nm/rad–5.25 Nm/rad]. Obtained maximal values of joint stiffness are aligned with the values of elbow joint stiffness presented in papers [15] (from 3 to 20 Nm/rad), [46] (from 2.5 to 9 Nm/rad), [47] (from 9 to 14 Nm/rad), and [48] (~10 Nm/rad).

Due to the variability in biological mechanisms and experimental conditions, we qualitatively compared the stiffness evolution in our work with results from the relevant literature in cases where the human arm is moving with different velocities and where it is acting against different forces. Humans perform interaction tasks with ease and precision by modulating their mechanical stiffness. It is intuitively expected and proven in the literature that the co-contraction of antagonistic muscle groups (which is highly related to joint mechanical impedance) depends on velocity [49] or load [50]. However, for a number of challenges that include human motion analysis, rehabilitation, and human–robot interaction, comprehensive study and understanding of the mutual influence of velocity and load on human motion and stiffness planning are needed. Furthermore, considering the latest enabling technology—collaborative robots capable of mechanical impedance modulation—the impact that cobots can have on healthcare and industry highly depends on our capabilities to understand human behavior patterns in contact and non-contact tasks.

The stiffness estimation in tasks E1 and E2 is presented in Figure 6, Figure 7, Figure 8 and Figure 9. To the best of our knowledge, no report analyzes the stiffness in such a way in the existing literature. The obtained results for experiment E1 (Figure 6) show that the stiffness value increases with the increase in the movement tempo and that the same effect was observed in all participants. Therefore, it can be concluded that there is a universal relationship that connects the stiffness and the movement tempo, and this dependence can be represented by Equation (30). The corresponding curve is shown in Figure 7 and it represents the functional scale of the stiffness change regarding the movement tempo change in the population without sensorimotor disorders. The obtained functional scale can be used as a reference curve for assessing the recovery of patients with upper extremity weakness. Compared to the papers [46,49], which analyze the correlation between joint stiffness and joint velocity, we obtained quite similar results: the increase in speed affects the increase in stiffness. A similar stiffness-increasing trend with increasing weight in the E2 task is presented in Figure 8. It can be observed that in participants 5, 7, and 10, stiffness has an increasing trend for all increasing weights, including 1 kg. Upon comparing the stiffness–force relationship with the results in [8], it can be noticed that with increasing force, stiffness also increases, which complies with the conclusions of [8]. However, in most cases (participants 1, 2, 3, 4, 6, 8, and 9), stiffness has an increasing trend for all increasing weights except for 1 kg, where the stiffness value is close to the stiffness at 0.75 kg. A possible reason for this phenomenon is the saturation of stiffness in the elbow joint with an increase in weights greater than 0.75 kg. This saturation property of elbow joint stiffness can be explained by the compensatory ability of motor behavior [47]. Based on the obtained results for different weights, the functional scale of stiffness change regarding weight change can be defined and is given in Equation (31). The corresponding curve is given in Figure 9. Due to the previously mentioned saturation effect, the stiffness deviation at 1 kg is larger than for other weights. Considering that the 0.5 kg weight is conventionally used in task 8 of WMFT, the obtained functional scale can also be used as a reference for assessing the recovery of patients with upper extremity weakness [51].

5. Conclusions

This paper introduced a general assessment methodology for elbow joint stiffness estimation based on Hill’s muscle model and the passive joint structure model, and a genetic optimization procedure during the “Reach and retrieve” task that is conventionally used within neurorehabilitation testing. The model parameters were estimated using the short-term tempo ramp scenario, and then, verified in two independent scenarios, resulting in 6.1% and 8.1% errors for the training and test datasets, respectively. The estimated joint stiffness in the two long-term test scenarios (the first one with tempo variation and the second one with load variation) was used for creating elbow stiffness curves to represent the general pattern of motor behavior in the population without sensorimotor disorders. The proposed objectification criteria of motor behavior in the form of stiffness patterns could significantly contribute to the monitoring of motor recovery in patients with sensorimotor disorders. A comparison to the results obtained using the proposed model with existing clinical scales (for the selected task as well as for the overall functional assessment of the arm) would be of interest for future work, as would further automatization of the overall system using collaborative robot-based approaches. Understanding human stiffness patterns and stiffness monitoring is a prerequisite for the development of rehabilitation devices based on collaborative robots capable of imposing recommended stiffness patterns, on the one hand, and monitoring Cartesian stiffness in contact tasks, on the other hand. Future research will exploit the presented methodology in clinical studies for assessment of the rehabilitation process with a larger number of participants.

Author Contributions

Conceptualization, M.R., D.U., M.M.J., S.D.D., T.J.D.T. and K.J.; methodology, M.R., D.U., M.M.J., S.D.D., T.J.D.T. and K.J.; software, M.R., D.U., M.M.J. and M.T.; formal analysis, M.R., D.U., M.M.J., M.T. and K.J.; data acquisition, M.R., D.U., M.M.J., S.D.D., T.J.D.T. and M.T.; resources, M.R., D.U., M.M.J., S.D.D., T.J.D.T., M.T. and K.J.; data curation, M.R., D.U., M.M.J., S.D.D., T.J.D.T., M.T. and K.J.; writing—original draft preparation, M.R., D.U., M.M.J., S.D.D., M.T. and K.J.; writing—review and editing, M.R., D.U., M.M.J., S.D.D., T.J.D.T. and K.J.; visualization, M.R., D.U. and M.M.J.; project administration, K.J. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

This study was conducted according to the guidelines of the Declaration of Helsinki and approved by the ethics committee of the University of Belgrade—School of Electrical Engineering, Serbia.

Informed Consent Statement

Informed consent was obtained from all subjects involved in the study.

Data Availability Statement

All information is given on project official website http://fornextcobot.etf.bg.ac.rs/sr/naslovna/en/naslovna-english/.

Acknowledgments

We would like to thank to Ljubica Konstantinović from the Clinic for Rehabilitation and “Miroslav Zotović”, Belgrade, for their valuable comments regarding the study’s conceptualization and draft preparation. Additionally, we would like to thank all volunteers for their participation in the study.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Experimental setup: (A) initial position of the arm; (B) the position of the arm during the motion. 1—EMG electrodes located on triceps brachii long head muscle, 2—EMG electrodes located on biceps brachii long head muscle, 3 and 5—IMU units, 4—goniometer, 6—wrist fixation, 7—feedback display of the movement tempo.

Figure 2. Muscle activation estimation: (A) Processed EMG signal (blue line) and EMG envelope (red line); (B) neural activation [Forumla omitted. See PDF.]; (C) muscle activation [Forumla omitted. See PDF.] obtained using Equation (1).

Figure 3. A block diagram of the joint stiffness estimation.

Figure 4. Hill’s muscle model.

Figure 5. An example of WMFT8m movement pattern: mean value (bold line) and standard deviation (shaded area) of (A) measured and estimated elbow angle; (B) elbow joint stiffness; (C) total joint torque. The presented signals belong to participant 9 for task E1 and have a tempo of 30 bpm.

Figure 6. Comparative analysis of (A) elbow joint stiffness mean value for participant 9; (B) stiffness value for each tempo for all 10 participants in task E1.

Figure 7. Stiffness fitting curve (green curve) for different task tempos (E1). Mean value for all participants is represented by red dots and the corresponding standard deviation is represented by blue lines.

Figure 8. Comparative analysis of (A) elbow joint stiffness mean value for participant 9; (B) stiffness value for each tempo for all 10 participants in task E2.

Figure 9. Stiffness fitting curve (green curve) for different “pulling” weights task (E2). Mean value for all participants is presented by red dots and corresponding standard deviation is presented by blue lines.

Table 1

Selected body segment parameters (BSP).

No	Parameter	No	Parameter
1	Length, Hand	6	Circumference, Forearm
2	Length, Wrist to Knuckle	7	Circumference, Elbow
3	Length, Forearm	8	Width, Hand
4	Circumference, Fist	9	Width, Wrist
5	Circumference, Wrist	10	Width, Elbow

Table 2

Parameters of Hill’s model and passive joint structure model for the elbow joint.

Parameter	Mean Value ± Standard Deviation	Min Value	Max Value
k1 [Nm]	7.16 ± 1.88	4.54	9.05
k2 [${\mathrm{rad}}^{-1}$]	0.35 ± 0.05	0.25	0.43
k4 [Nm]	35.89 ± 13.24	15.00	60.5
k5 [${\mathrm{rad}}^{-1}$]	0.026 ± 0.006	0.02	0.04
q1 [$\mathrm{rad}$]	4.20 ± 0.77	3.11	5.49
q2 [$\mathrm{rad}$]	5.25 ± 1.29	3.22	6.95
Bp [Nms]	3.32 ± 0.81	1.55	4.05
Fmax [N]	2289.76 ± 786.18	1401.59	3913.00
w [n.u.]	0.78 ± 0.14	0.64	0.96
lmopt [m]	0.19 ± 0.05	0.10	0.27
af [n.u.]	0.39 ± 0.19	0.25	0.75
kte [${\mathrm{m}}^{-1}$]	2611.95 ± 862.11	1511.02	4042.59
ktl [N/m]	27,230.04 ± 8846.59	15,570.31	43,133.01
kt [N/m]	538,599.51 ± 196,458.08	317,771.71	923,909.99
ltc [m]	0.34 ± 0.093	0.23	0.47
lts [m]	0.33 ± 0.070	0.24	0.46
kml [N/m]	408.21 ± 137.80	244.96	710.60
kme [${\mathrm{m}}^{-1}$]	75.84 ± 23.35	46.05	113.92
km [N/m]	4551.31 ± 1590.63	2500.06	7352.47
lmc [m]	0.14 ± 0.03	0.10	0.19
lms [m]	0.19 ± 0.05	0.10	0.27
Bm [Ns/m]	0.19 ± 0.05	130.97	342.17
${\mathrm{M}}_{\mathrm{m}}\left[\mathrm{kg}\right]$	241.77 ± 75.83	0.33	0.45
${\mathrm{v}}_{\max }\left[\mathrm{m}/\mathrm{s}\right]$	0.39 ± 0.04	0.54	1.46

Table 3

Verification of the Training GA outputs using DTW—average deviation (in degrees).

		1	2	3	4	5	6	7	8	9	10
Training phase		8.47	4.46	8.02	1.95	8.67	9.16	4.89	2.56	2.77	1.83
E1	15 bpm	2.78	5.21	8.43	4.37	5.42	5.48	6.67	1.63	3.55	5.80
	30 bpm	4.45	1.98	2.91	4.05	5.16	5.83	1.78	5.25	2.96	2.16
	45 bpm	4.80	3.11	6.05	2.06	6.32	6.51	2.43	2.23	2.95	2.31
	60 bpm	5.31	5.74	10.79	6.56	11.85	7.03	3.13	4.76	3.65	0.99
E2	0 kg	4.55	7.80	8.47	2.52	4.35	12.16	11.43	5.10	4.75	6.94
	0.25 kg	4.76	6.55	6.64	4.91	5.33	3.53	10.15	6.29	2.23	3.02
	0.5 kg	4.45	1.98	2.91	4.05	5.16	5.83	1.78	5.25	2.96	2.16
	0.75 kg	3.70	7.22	11.22	3.20	3.93	8.65	5.49	1.69	2.93	11.14
	1 kg	8.13	5.44	10.79	4.38	5.35	8.82	5.32	6.96	2.30	3.26

Table 4

Verification of the accuracy of the generated stiffness patterns depending on the change in parameter values using DTW—average deviation (in Nm/rad).

Parameter	Percentage Variation
	−10%	+10%	−20%	20%	−30%	30%
k1	0.0622	0.0429	0.1784	0.1370	0.2854	0.2916
k2	0.1712	0.1341	0.3541	0.3405	0.4556	0.5504
k4	0.0331	0.0353	0.1078	0.0916	0.2543	0.1592
k5	0.0019	0.0019	0.0040	0.0038	0.0064	0.0059
q1	0.0024	0.0024	0.0052	0.0050	0.0084	0.0077
q2	0.1669	0.1934	0.3211	0.3747	0.4050	0.3827
Bp	0.0041	0.0040	0.0100	0.0094	0.0173	0.0155
Fmax	0.0239	0.0243	0.0538	0.0614	0.1028	0.1066
w	0.0142	0.0104	0.0446	0.0235	0.0994	0.0367
lmopt	0.0444	0.0287	0.0938	0.0548	0.1887	0.0899
af	0.0193	0.0159	0.0559	0.0377	0.0994	0.0609
kte	0.00035	0.00042	0.0005	0.00064	0.00057	0.00072
ktl	0.00028	0.00033	0.00037	0.00036	0.00047	0.00047
kt	0.00074	0.00071	0.0012	0.00087	0.0020	0.0011
ltc	5.0826	0.0038	7.1025	0.0038	5.8225	0.0038
lts	0.0038	4.9648	0.0038	6.9452	0.0038	6.2682
kml	0.00028	0.00031	0.00043	0.00044	0.0005	0.00051
kme	0.00032	0.00026	0.00046	0.00039	0.00056	0.00046
km	0.0074	0.0072	0.0186	0.0168	0.0322	0.0277
lmc	0.0082	0.0133	0.0110	0.0215	0.0114	0.0298
lms	0.0444	0.0287	0.0938	0.0548	0.1887	0.0899
Bm	0.00079	0.00079	0.0013	0.0012	0.0018	0.0018
${\mathrm{M}}_{\mathrm{m}}$	0.00005	0.00009	0.00007	0.00014	0.00009	0.00017
${\mathrm{v}}_{\max }$	0.0254	0.0215	0.0733	0.0517	0.1440	0.0846

Table 5

Comparative analysis of the Hill’s muscle model and passive joint structure model parameters.

Reference	lmopt [m]	lts [m]	lms [m]	Fmax [N]	Mm [kg]	af [n.u.]	vmax [m/s]
Cavallaro et al. [20]	0.15	0.19–0.23	n.a.	393–1000	n.a.	n.a.	n.a.
Garner et al. [21]	0.15	0.19–0.23	n.a.	462–629	n.a.	n.a.	n.a.
Winters et al. [22]	0.09–0.14	0.22	0.31–0.37	n.a.	0.06–0.13	0.35–0.45	0.49–0.76
Desplenter et al. [23]	0.11–0.17	n.a.	0.16–0.25	2875–2397	n.a.	n.a.	n.a.
Holzbaur et al. [24]	0.11–0.13	0.14–0.27	n.a.	624–799	n.a.	n.a.	n.a.
Murray et al. [25]	0.13	0.22–0.23	0.22–0.36	n.a.	n.a.	n.a.	n.a.
Yu et al. [26]	0.15	0.11–0.13	n.a.	n.a.	n.a.	n.a.	n.a.
Song et al. [27]	0.16–0.17	0.21–0.23	0.39–0.41	630–764	0.34–0.45	n.a.	n.a.
Bingshan et al. [28]	0.12–0.13	0.14–0.27	n.a.	624–799	n.a.	n.a.	0.96–1.04
This paper	0.1–0.27	0.24–0.46	0.1–0.27	1400–3900	0.33–0.45	0.25–0.75	0.54–1.46

Appendix A

The specific measurement data (Processed EMG signal and its envelope) for tasks E1 and E2 for participant 9 are presented in Figure A1 and Figure A2, respectively. The results for task E1 are presented in Figure A1. The EMG signals and envelopes are represented for a time sequence of 4 s. The different tempo of the motion is evident in the EMG signals and envelopes. For a tempo of 15 bpm during 4 s, the motion is performed once, while for a tempo of 60 bpm during 4 s, the motion is performed four times. Additionally, increasing the tempo affects the values of the Processed EMG signal and its envelope (see Figure A1). The maximal envelope value for a tempo of 15 bpm is ~0.02 mV, for a tempo of 30 bpm is ~0.03 mV, for a tempo of 45 bpm is 0.05 mV, and for a tempo of 60 bpm is 0.1 mV.

Figure A1. Muscle activation estimation: Processed EMG signal (blue line) and EMG envelope (red line) in task E1.

Figure A2. Muscle activation estimation: Processed EMG signal (blue line) and EMG envelope (red line) in task E2.

In task E2, increasing the load affects the values of the Processed EMG signal and its envelope (see Figure A2). The EMG signals and envelopes are represented for the time sequence of 4 s. Because task E2′s tempo is 30 bpm, the motion is repeated four times for 4 s for each load. The results show that the EMG signals and envelope increase with increased load. The maximal envelope value for a load of 0 kg is ~0.025 mV, for a load of 0.25 kg is ~0.04 mV, for a load of 0.5 kg is ~0.06 mV, for a load of 0.75 kg is ~0.07 mV, and for a load of 1 kg is ~0.1 mV.

References

1. Tran, V.D.; Dario, P.; Mazzoleni, S. Kinematic Measures for Upper Limb Robot-Assisted Therapy Following Stroke and Correlations with Clinical Outcome Measures: A review. Med. Eng. Phys.; 2018; 53, pp. 13-31. [DOI: https://dx.doi.org/10.1016/j.medengphy.2017.12.005]

2. Gladstone, D.J.; Danells, C.J.; Black, S.E. The Fugl-Meyer Assessment of Motor Recovery After Stroke: A Critical Review of its Measurement Properties. Neurorehabil. Neural Repair; 2002; 16, pp. 232-240. [DOI: https://dx.doi.org/10.1177/154596802401105171]

3. Bernhardt, J.; Hayward, K.S.; Kwakkel, G.; Ward, N.S.; Wolf, S.L.; Borschmann, K.; Krakauer, J.W.; Boyd, L.A.; Carmichael, S.T.; Corbett, D. et al. Agreed Definitions and A Shared Vision for New Standards in Stroke Recovery Research: The Stroke Recovery and Rehabilitation Roundtable Taskforce. Int. J. Stroke; 2017; 12, pp. 444-450. [DOI: https://dx.doi.org/10.1177/1747493017711816] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28697708]

4. Bernaldo de Quirós, M.; Douma, E.H.; van den Akker-Scheek, I.; Lamoth, C.J.C.; Maurits, N.M. Quantification of Movement in Stroke Patients under Free Living Conditions Using Wearable Sensors: A Systematic Review. Sensors; 2022; 22, 1050. [DOI: https://dx.doi.org/10.3390/s22031050] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35161796]

5. Kwakkel, G.; Lannin, N.A.; Borschmann, K.; English, C.; Ali, M.; Churilov, L.; Saposnik, J.W.; Winstein, C.; van Wegen, E.E.; Wolf, S.L. et al. Standardized Measurement of Sensorimotor Recovery in Stroke Trials: Consensus-Based Core Recommendations from the Stroke Recovery and Rehabilitation Roundtable. Int. J. Stroke; 2017; 31, pp. 784-792. [DOI: https://dx.doi.org/10.1177/1747493017711813]

6. Borzelli, D.; Pastorelli, S.; d’Avella, A.; Gastaldi, L. Virtual Stiffness: A Novel Biomechanical Approach to Estimate Limb Stiffness of a Multi-Muscle and Multi-Joint System. Sensors; 2023; 23, 673. [DOI: https://dx.doi.org/10.3390/s23020673] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36679467]

7. Li, K.; Zhang, J.; Wang, L.; Zhang, M.; Li, J.; Bao, S. A Review of The Key Technologies for sEMG-Based Human-Robot Interaction Systems. Biomed. Signal Process. Control; 2020; 62, 102074. [DOI: https://dx.doi.org/10.1016/j.bspc.2020.102074]

8. Li, K.; Liu, X.; Zhang, J.; Zhang, M.; Hou, Z. Continuous Motion and Time-Varying Stiffness Estimation of the Human Elbow Joint Based on SEMG. J. Mech. Med. Biol.; 2019; 19, 1950040. [DOI: https://dx.doi.org/10.1142/S0219519419500404]

9. He, B.; Huang, H.; Thomas, G.C.; Sentis, L. A Complex Stiffness Human Impedance Model with Customizable Exoskeleton Control. IEEE Trans. Neural Syst. Rehabil. Eng.; 2020; 28, pp. 2468-2477. [DOI: https://dx.doi.org/10.1109/TNSRE.2020.3027501]

10. Perreault, E.J.; Kirsch, R.F.; Crago, P.E. Voluntary Control of Static Endpoint Stiffness During Force Regulation Tasks. J. Neurophysiol.; 2002; 87, pp. 2808-2816. [DOI: https://dx.doi.org/10.1152/jn.2002.87.6.2808]

11. Perreault, E.; Kirsch, R.; Crago, P. Multijoint Dynamics and Postural Stability of the Human Arm. Exp. Brain Res.; 2004; 157, pp. 507-517. [DOI: https://dx.doi.org/10.1007/s00221-004-1864-7] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15112115]

12. Ajoudani, A.; Fang, C.; Tsagarakis, N.G.; Bicchi, A.A. A reduced-Complexity Description of Arm Endpoint Stiffness with Applications to Teleimpedance Control. Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2015); Hamburg, Germany, 28 September–2 October 2015; [DOI: https://dx.doi.org/10.1109/iros.2015.7353495]

13. Ajoudani, A.; Fang, C.; Tsagarakis, N.; Bicchi, A. Reduced-Complexity Representation of the Human Arm Active Endpoint Stiffness for Supervisory Control of Remote Manipulation. Int. J. Robot. Res.; 2017; 37, pp. 155-167. [DOI: https://dx.doi.org/10.1177/0278364917744035]

14. Fang, C.; Ajoudani, A.; Bicchi, A.; Tsagarakis, N.G. Online Model Based Estimation of Complete Joint Stiffness of Human Arm. IEEE Robot. Autom. Lett.; 2018; 3, pp. 84-91. [DOI: https://dx.doi.org/10.1109/LRA.2017.2731524]

15. Stanev, D.; Moustakas, K. Stiffness Modulation of Redundant Musculoskeletal Systems. J. Biomech.; 2019; 85, pp. 101-107. [DOI: https://dx.doi.org/10.1016/j.jbiomech.2019.01.017]

16. Wu, Y.; Zhao, F.; Kim, W.; Ajoudani, A. An Intuitive Formulation of the Human Arm Active Endpoint Stiffness. Sensors; 2020; 20, 5357. [DOI: https://dx.doi.org/10.3390/s20185357]

17. Liu, J.; Ren, Y.; Xu, D.; Kang, S.H.; Zhang, L.Q. EMG-Based Real-Time Linear-Nonlinear Cascade Regression Decoding of Shoulder, Elbow, and Wrist Movements in Able-Bodied Persons and Stroke Survivors. IEEE Trans. Biomed. Eng.; 2019; 67, pp. 1272-1281. [DOI: https://dx.doi.org/10.1109/TBME.2019.2935182] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31425016]

18. Jovanović, K.; Vranić, J.; Miljković, N. Hill's and Huxley's Muscle Models: Tools for Simulations in Biomechanics. Serbian J. Electr. Eng.; 2015; 12, pp. 53-67. [DOI: https://dx.doi.org/10.2298/SJEE1501053J]

19. Robleto, R. An Analysis of the Musculotendon Dynamics of Hill-Based Models. Master’s thesis; Texas Tech University Libraries: 2500 Broadway Lubbock, TX, USA, 1997.

20. Cavallaro, E.E.; Rosen, J.; Perry, J.C.; Burns, S. Real-Time Myoprocessors for a Neural Controlled Powered Exoskeleton Arm. IEEE Trans. Bio-Med. Eng.; 2006; 53, pp. 2387-2396. [DOI: https://dx.doi.org/10.1109/TBME.2006.880883]

21. Garner, B.A.; Pandy, M.G. Musculoskeletal Model of the Upper Limb Based on the Visible Human Male Dataset. Comput. Methods Biomech. Biomed. Eng.; 2001; 4, pp. 93-126. [DOI: https://dx.doi.org/10.1080/10255840008908000]

22. Winters, J.M.; Stark, L. Estimated Mechanical Properties of Synergistic Muscles Involved in Movements of a Variety of Human Joints. J. Biomech.; 1988; 21, pp. 1027-1041. [DOI: https://dx.doi.org/10.1016/0021-9290(88)90249-7] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/2577949]

23. Desplenter, T.; Trejos, A. Evaluating Muscle Activation Models for Elbow Motion Estimation. Sensors; 2018; 18, 1004. [DOI: https://dx.doi.org/10.3390/s18041004] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29597281]

24. Holzbaur, K.R.S.; Murray, W.M.; Delp, S.L. A Model of the Upper Extremity for Simulating Musculoskeletal Surgery and Analyzing Neuromuscular Control. Ann. Biomed. Eng.; 2005; 33, pp. 829-840. [DOI: https://dx.doi.org/10.1007/s10439-005-3320-7] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16078622]

25. Murray, W.M.; Buchanan, T.S.; Delp, S.L. The Isometric Functional Capacity of Muscles that Cross the Elbow. J. Biomech.; 2000; 33, pp. 943-952. [DOI: https://dx.doi.org/10.1016/S0021-9290(00)00051-8] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10828324]

26. Yu, T.F.; Wilson, A.J. A Passive Movement Method for Parameter Estimation of a Musculoskeletal Arm Model Incorporating a Modified Hill Muscle Model. Comput. Methods Progr. Biomed.; 2014; 114, pp. e46-e59. [DOI: https://dx.doi.org/10.1016/j.cmpb.2013.11.003] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24290234]

27. Song, D.; Lan, N.; Loeb, G.E.; Gordon, J. Model-Based Sensorimotor Integration for Multi-Joint Control: Development of a Virtual Arm Model. Ann. Biomed. Eng.; 2008; 36, pp. 1033-1048. [DOI: https://dx.doi.org/10.1007/s10439-008-9461-8] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18299994]

28. Bingshan, H.; Haoran, T.; Hongrun, L.; Xiangxiang, Z.; Jiantao, Y.; Hongliu, Y. An Improved EMG-Driven Neuromusculoskeletal Model for Elbow Joint Muscle Torque Estimation. Appl. Bionics Biomech.; 2021; 2021, 1985741. [DOI: https://dx.doi.org/10.1155/2021/1985741]

29. Huang, Y.; Chen, K.; Zhang, X.; Wang, K.; Ota, J. Motion Estimation of Elbow Joint From sEMG using Continuous Wavelet Transform and Back Propagation Neural Networks. Biomed. Signal Process. Control; 2021; 68, 102657. [DOI: https://dx.doi.org/10.1016/j.bspc.2021.102657]

30. Lemerle, S.; Grioli, G.; Bicchi, A.; Catalano, M.G. A Variable Stiffness Elbow Joint for Upper Limb Prosthesis. Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2019); Macau, China, 4–8 November 2019; pp. 7327-7334.

31. Pearsall, D.J.; Reid, G. The Study of Human Body Segment Parameters in Biomechanics. Sports Med.; 1994; 18, pp. 126-140. [DOI: https://dx.doi.org/10.2165/00007256-199418020-00005]

32. Wolf, S.L.; Lecraw, D.E.; Barton, L.A.; Jann, B.B. Forced Use of Hemiplegic Upper Extremities to Reverse the Effect of Learned Nonuse Among Chronic Stroke and Head-Injured Patients. Exp. Neurol.; 1989; 104, pp. 125-132. [DOI: https://dx.doi.org/10.1016/S0014-4886(89)80005-6]

33. Cestelein, B.; Cagnie, B. Parlevliet Danneels, T.L.; Cools, A. Optimal Normalization Tests for Muscle Activation of the Levator Scapulae, Pectoralis Minor, and Rhomboid Major: An Electromyography Study Using Maximum Voluntary Isometric Contractions. Arch. Phys. Med. Rehabil.; 2015; 96, pp. 1820-1827. [DOI: https://dx.doi.org/10.1016/j.apmr.2015.06.004]

34. Nieminen, H.; Taklal, E.-P.; Viikari-Juntura, E. Normalization of Electromyogram in the Neck-Shoulder Region. Eur. Appl. J. Physiol. Occup. Physiol.; 1993; 67, pp. 199-207. [DOI: https://dx.doi.org/10.1007/BF00864215]

35. Ekstrom, R.A.; Soderberg, G.L.; Donatelli, R.A. Normalization Procedures Using Maximum Voluntary Isometric Contractions for the Serratus Anterior and Trapezius Muscles During Surface EMG Analysis. J. Electromyogr. Kinesiol.; 2004; 15, pp. 418-428. [DOI: https://dx.doi.org/10.1016/j.jelekin.2004.09.006]

36. Halaki, M.; Ginn, K. Normalization of EMG Signals: To Normalize or Not to Normalize and What to Normalize to?. Computational Intelligence in Electromyography Analysis—A Perspective on Current Applications and Future Challenges; InTech: Rijeka, Croatia, 2012; Chapter 7 ISBN 978-953-51-0805-4

37. Manal, K.; Buchanan, T.S. A One-Parameter Neural Activation to Muscle Activation Model: Estimating Isometric Joint Moments from Electromyograms. J. Biomech.; 2003; 36, pp. 1197-1202. [DOI: https://dx.doi.org/10.1016/S0021-9290(03)00152-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12831746]

38. Hill, A.V. The Heat of Shortening and the Dynamic Constants of Muscle. Proc. R. Soc.; 1938; 126, pp. 136-195. [DOI: https://dx.doi.org/10.1098/rspb.1938.0050]

39. He, J. A Feedback Control Analysis of the Neuro-Musculo-Skeletal System of a Cat Hindlimb. Ph.D. Dissertation; The University of Maryland: College Park, MD, USA, 1988.

40. Audu, M.L.; Davy, D.T. The Influence of Muscle Model Complexity in Musculoskeletal Motion Modeling. J. Biomech. Eng.; 1985; 107, pp. 147-157. [DOI: https://dx.doi.org/10.1115/1.3138535]

41. Sakoe, H.; Seibi, C. Dynamic Programming Algorithm Optimization for Spoken Word Recognition. IEEE Trans. Acoust. Speech Signal Process.; 1978; 26, pp. 43-49. [DOI: https://dx.doi.org/10.1109/TASSP.1978.1163055]

42. Fagiolini, A.; Trumić, M.; Jovanović, K. An Input Observer-Based Stiffness Estimation Approach for Flexible Robot Joints. IEEE Robot. Autom. Lett.; 2020; 5, pp. 1843-1850. [DOI: https://dx.doi.org/10.1109/LRA.2020.2969952]

43. Trumić, M.; Grioli, G.; Jovanović, K.; Fagiolini, A. Force/Torque-Sensorless Joint Stiffness Estimation in Articulated Soft Robots. IEEE Robot. Autom. Lett.; 2022; 7, pp. 7036-7043. [DOI: https://dx.doi.org/10.1109/LRA.2022.3178467]

44. Ljung, L. System identification. Signal Analysis and Prediction; Procházka, A.; Uhlíř, J.; Rayner, P.W.J.; Kingsbury, N.G. Springer, Birkhäuser: Boston, MA, USA, 1998; pp. 163-173. [DOI: https://dx.doi.org/10.1007/978-1-4612-1768-8] ISBN 978-1-4612-1768-8

45. Zhao, Y.; Zhang, Z.; Li, Z.; Yang, Z.; Dehghani-Sanij, A.A.; Xie, S. An EMG-Driven Musculoskeletal Model for Estimating Continuous Wrist Motion. IEEE Trans. Neural Syst. Rehabil. Eng.; 2020; 28, pp. 3113-3120. [DOI: https://dx.doi.org/10.1109/TNSRE.2020.3038051] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33186119]

46. Bennett, D.J.; Hollerbach, J.M.; Xu, J.; Hunter, I.W. Time-Varying Stiffness of Human Elbow Joint During Cyclic Voluntary Movement. Exp. Brain Res.; 1992; 88, pp. 433-442. [DOI: https://dx.doi.org/10.1007/BF02259118]

47. Shin, W.; Handdeut, C.; Kim, J. Human elbow motor learning skills of varying loads: Proof of internal model generation using joint stiffness estimation. J. Biomech. Sci. Eng.; 2021; 16, pp. 21-88. [DOI: https://dx.doi.org/10.1299/jbse.21-00088]

48. Franklin, D.W.; Burdet, E.; Osu, R.; Kawato, M.; Milner, T.E. Functional Significance of Stiffness in Adaptation of Multijoint Arm Movements to Stable and Unstable Dynamics. Exp. Brain Res.; 2003; 151, pp. 145-157. [DOI: https://dx.doi.org/10.1007/s00221-003-1443-3] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12783150]

49. Suzuki, M.; Shiller, D.M.; Gribble, P.L.; Ostry, D.J. Relationship Between Cocontraction, Movement Kinematics and Phasic Muscle Activity in Single-Joint Arm Movement. Exp. Brain Res.; 2021; 140, pp. 171-181. [DOI: https://dx.doi.org/10.1007/s002210100797] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11521149]

50. De Luca, C.J.; Mambrito, B. Voluntary Control of Motor Units in Human Antagonist Muscles: Coactivation and Reciprocal Activation. J. Neurophysiol.; 1987; 58, pp. 525-542. [DOI: https://dx.doi.org/10.1152/jn.1987.58.3.525]

51. Schwarz, A.; Kanzler, C.M.; Lambercy, O.; Luft, A.R.; Veerbeek, J.M. Systematic Review on Kinematic Assessments of Upper Limb Movements After Stroke. Stroke; 2019; 50, pp. 718-727. [DOI: https://dx.doi.org/10.1161/STROKEAHA.118.023531]