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

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.


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.

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 Sections 2.4.1 and 2.4.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. 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.

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 (rcm), 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 Sections 2.4.1 and 2.4.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.

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, Sensors 2023, 23, 1709 6 of 22 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 u(t) ( Figure 2). Muscle activation a(t) was calculated using the model with the lowest data point optimization time [23], as suggested by Manal et Buchanan [37], given in Equation (1): where A = −20 was taken to achieve the range of muscle activation function [0,1].

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 u t ( Figure 2). Muscle activation a t was calculated using the model with the lowest data point optimization time [23], as suggested by Manal et Buchanan [37], given in Equation (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.  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. 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 (a AG (t)) and triceps brachii long head (a 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 ( . 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. chii long head muscle (a (t)) and triceps brachii long head (a (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 ( ). 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 (qgoni, 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.

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 (lm), while the coupling force F will increase the length of the muscle.

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.
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 ( ). 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 (qgoni, 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.

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 (lm), while the coupling force F will increase the length of the muscle.  where F max is the maximum isometric force (the maximum force in the static state at MVC), a(t) is the muscle activation function, v m is the actual muscle velocity, v max is the maximum unloaded shortening speed, a f is the Hill's model parameter, l m is the actual muscle fiber length, l mopt is the optimal muscle fiber length, and w is the width of the force-length relationship. The total force F PE of the PE element is the sum of the forces delivered by the stiffness (F p ) and damping (F d ) force (Equation (5)).
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): where l m is the actual muscle fiber length, l ms is the passive muscle slack length, l mc is the length at which the spring becomes linear and is equal to the optimal muscle length l mopt , k m and k ml are parallel spring constants, and k 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): where B m is the parallel damping constant and v 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): l t < l ts , k tl k te (exp(k te (l t − l ts )) − 1), l ts ≤ l t ≤ l tc F tc + k t (l t − l tc ), l t < l tc (9) F tc = k tl k te (exp(k te (l tc − l ts )) − 1) where k t and k tl represent tendon spring constants, k te is the exponential shape parameter, l t is the tendon length, l ts is the tendon slack length, and l 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 l m , muscle speed v m , and tendon length l t are given in Equations (11)-(13): l t = l ms +l ts + r p q 0 − l m (13) where M m is the muscle mass, r p is the radius of the elbow joint calculated as a half value of the BSP parameter P10 from Table 1, and q 0 is the elbow joint angle during the full extension (q 0 = 0 deg). The joint torque τ t generated from the muscle and transmitted through the tendon F SE via r p is given in Equation (14).

Passive Joint Structures
Passive joint structure implies stability in the musculoskeletal model [39,40]. The passive joint torque τ p , the damping torque τ d , and the torque originating from gravity τ G are given in t Equations (15)- (17): where m is the mass of the forearm-hand segment, B p is the joint damping constant, g is the gravitational constant, r cm is the distance between the joint axis and the center of mass of the forearm-hand, q 1 and q 2 represent joint angles where the torque begins to increase, k 1 and k 4 are joint torque constants, and k 2 are k 5 are dimensionless parameters. The parameters k 1 , k 2 , k 4 , and 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 τ J is the sum of the torques originating from the passive joint structure and the torque induced by the antagonist's muscle τ tANT subtracted by the torque produced by the agonist's muscle τ tAG (Equation (18)).
According to Newton's Second Law for rotation, Equation (16) qdt + q 0 (21) where M is the moment of inertia of the forearm-hand segment.

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 Sections 2.4.1 and 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: where q goni is the joint angle measured using a goniometer, q is the estimated joint angle according to Equation (21), t 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. Roulettewheel 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.

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.

Recursive Least Square Algorithm
Analogously to the approach presented in [42,43], the RLS algorithm was applied to the estimated joint torque τ J and joint angle q to estimate the elbow joint stiffness. The total joint torque τ J was expressed as in Equation (23): where Φ = 1, q, q 2 , . . . , q n is the n order torque regressor, which is the function of the joint angle, and Π = (π 1 , π 2 , . . . π n ) T is a vector of the state parameters to be estimated. The order 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 Π is proposed in [44] and given by Equations (24)- (28):  (29).

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

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 (τ 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 1 , k 2 , k 4 , q 2 , 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.  Figure 5 shows an example of the typical pattern for WMFT8m movement. The mea value (bold line) and standard deviation (shaded area) of these signals were obtained b overlapping all WMFT8m movements and their normalization with the time: (A) measure (qgoni) and estimated (q) elbow joint angles; (B) estimated stiffness value (Stiffness); and (C estimated total joint torque value (τ ). The comparative analysis of stiffness accuracy in relation to the change in model pa rameter values is presented in Table 4. The values of the model parameters were modifie for ±10%, ±20%, and ±30% regarding the parameter values estimated by GA. The accurac of the stiffness was calculated using the DTW algorithm, comparing the estimated stiff ness 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 k1, k2, k4, q2, ltc, lts, Fmax, lmopt, lms, and vm

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 tempo is given in Equation (30) and shown in Figure 7. A stiffness fitting curve (a so-called functional scale) for different task tempos Sti f f ness (tempo) is given in Equation (30) and shown in Figure 7. 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 WMFT8m 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 weight are given in Equation (31) and shown in Figure 9. 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 Sti f f ness(weight) are given in Equation (31) and shown in Figure 9.

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

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], M m values are given in [22,27], a f values are presented in [22], and 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 , M m , a f , and 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 Figures 6A and 8A [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 humanrobot 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 noncontact tasks.
The stiffness estimation in tasks E1 and E2 is presented in Figures 6-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].

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.