Unifying system identification and biomechanical formulations for the estimation of muscle, tendon and joint stiffness during human movement

In vivo joint stiffness estimation during time-varying conditions remains an open challenge. Multiple communities, e.g. system identification and biomechanics, have tackled the problem from different perspectives and using different methods, each of which entailing advantages and limitations, often complementary. System identification formulations provide data-driven estimates of stiffness at the joint level, while biomechanics often relies on musculoskeletal models to estimate stiffness at multiple levels, i.e. joint, muscle, and tendon. Collaboration across these two scientific communities seems to be a logical step toward a reliable multi-level understanding of joint stiffness. However, differences at the theoretical, computational, and experimental levels have limited inter-community interaction. In this article we present a roadmap to achieve a unified framework for the estimation of time-varying stiffness in the composite human neuromusculoskeletal system during movement. We present our perspective on future developments to obtain data-driven system identification and musculoskeletal models that are compatible at the theoretical, computational, and experimental levels. Moreover, we propose a novel combined closed-loop paradigm, in which reference estimates of joint stiffness via system identification are decomposed into underlying muscle and tendon contribution via high-density-electromyography-driven musculoskeletal modeling. We highlight the need for aligning experimental requirements to be able to compare both joint stiffness formulations. Unifying both biomechanics’ and system identification’s formulations is a necessary step for truly generalizing stiffness estimation across individuals, movement conditions, training and impairment levels. From an application point of view, this is central for enabling patient-specific neurorehabilitation therapies, as well as biomimetic control of assistive robotic technologies. The roadmap we propose could serve as an inspiration for future collaborations across broadly different scientific communities to truly understand joint stiffness bio- and neuromechanics. Video Abstract: Unifying system identification and biomechanical formulations for the estimation of muscle, tendon and joint stiffness during human movement


Introduction
Human movement emerges from the coordinated interplay between neural, muscular and skeletal structures, interacting with the environment [1,2]. Despite extensive knowledge on anatomy and function of neural and musculoskeletal structures, there is a knowledge gap in the causal interplay, i.e. cause-effect relationships, between the neural and the biomechanical levels of movement [3][4][5]. The traditionally studied joint angles and net joint torques, if considered in isolation, cannot completely characterize how we respond to changes in the surrounding environment and how we interact with it [6]. For instance, for a given joint angle, the same net torque acting around a joint may be generated either by completely relaxing the muscles spanning the joint or by contracting the agonist and antagonist muscle groups simultaneously. However, in the latter case, it would require a larger external force/torque to change the joint angle by a fixed amount. To characterize the difference across the presented scenarios, joint impedance needs to be considered [7].
Human joint impedance is a fundamental neuromechanical property, closed-loop controlled by the central nervous system (CNS) and contributed by musculoskeletal properties [8]. It enables efficient and stable interaction with the environment and navigation across different types of terrain [9]. Joint impedance represents a joint's elastic-, viscous-and inertia-dependent dynamic behavior, thereby relating the torque τ ϕ (t) response to an externally imposed angular displacement perturbation ϕ(t) [10]. This article focuses on the elastic-dependent component of joint impedance, i.e. joint stiffness, which is the position-dependent component.
Reliable joint stiffness estimation during movement may help elucidate human motor control strategies and, consequently, guide the rationale for robust control of biomimetic prosthetic or exoskeletal robots [11]. Moreover, joint stiffness estimation in clinical settings may help deconstruct the biomechanical causes of common movement disorders (e.g. stiff-knee gait, or crouch gait), thereby improving neurorehabilitation and clinical decision-making procedures [12]. However, the joint's complex anatomy and underlying physiological and neuromechanical processes make joint stiffness estimation a challenging task. Different scientific communities, e.g. system identification and biomechanics, have faced the challenge by looking at the problem from different perspectives and at different levels, e.g. net joint, muscle, or tendon levels. Moreover, researchers within each community use different vocabulary, theoretical concepts, computational methods, and experimental protocols, each entailing strengths and weaknesses that are often complementary.
On the one hand, the system identification community largely uses a top-down approach, i.e. a data-driven estimation of joint impedance and joint stiffness is built from and validated on experimentally recorded perturbation-based data [13]. Since the experimental data is collected at the joint level, the contributions of all components spanning the joint are captured. However, a detailed physiological description is missing and, consequently, it is difficult to infer causal relationships between muscle-tendon unit (MTU) stiffness and overall joint stiffness. On the other hand, the biomechanics community often employs musculoskeletal models following a bottom-up approach, i.e. the constituting elements of a joint are modeled using basic principles of physics and physiology [14]. Joint stiffness, as well as the stiffness of the joint's constituting elements, are predicted from the musculoskeletal model's neuromechanical states, i.e. no external perturbations are required. However, a priori assumptions on the modeled elements must be made, and these elements are described by parameters that need to be calibrated across subjects using ex vivo measurements.
Despite the potential for merging both methods, the lack of a common framework to estimate joint stiffness has kept communication and cross-fertilization between communities to a low level. The aim of this perspective article is to propose a new paradigm to unify and generalize joint stiffness formulations during movement across the system identification and biomechanics communities. This will enable the definition of a common language and working framework, within which new methodologies can be developed for causally linking the neural and musculoskeletal mechanisms underlying joint stiffness. Although this paper's focus is joint stiffness, the paradigm is generalizable to joint impedance as well as to other estimation techniques, e.g. elastography.
This article is organized as follows. Section 2 includes a brief overview of joint stiffness from a physiological point of view and state-of-the-art methods in joint stiffness estimation. Section 3 describes the fundamentals of two selected methodologies from the system identification and the biomechanics communities. Section 4 presents our proposed solution to unify these techniques at three different levels, i.e. theoretical, computational, and experimental. Section 5 discusses open challenges and future steps. Finally, section 6 summarizes the article's take-home messages.

Joint stiffness determinants and estimation methodologies
This section provides a brief overview of the physiological mechanisms and anatomical structures contributing to joint stiffness. The objectives of this section are to provide the necessary physiological background to understand the remainder of the article and to define a common terminology across communities. Furthermore, an overview of the literature on joint stiffness estimation from the system identification and biomechanics communities is given. Key methodologies are selected and sorted in table 1 according to their required level of a priori knowledge.

Main determinants of joint stiffness 2.1.1. Musculoskeletal determinants
Joint stiffness is defined by the elastic properties of its constituent musculoskeletal elements including muscles, tendons, skeletal bones, ligaments, cartilages, as well as other connective and epithelial tissues (figure 1) [15].
The geometrical arrangement of fibers within a given muscle volume (i.e. parallel vs pennate muscles) dictates both fiber length range and maximal force-generating properties, thereby influencing muscle stiffness [16]. Muscle elastic properties also arise from both passive and active elements across subcellular to organ scales [17]. Muscle passive elements are those generating resistive mechanical force in response to induced stretches [18] and not being recruited by nervous system excitatory signals. Muscle active elements are intended as those being recruited by the nervous system, via afferent and efferent neural pathways, for regulating muscle contraction and mechanical force generation [18,19]. At the intersection of passive and active mechanisms there is the phenomenon called short-range stiffness (SRS), which is the instantaneous increase in force in response to an initial stretch due to the deformation of cross-bridges, i.e. activation-dependent links within the muscle's myofilaments [20]. The short-range comprises small deformations (up to 3% of the muscle fiber length) and a small amount of time, i.e. approximately within the first 50 ms after the stretch [21]. Muscles operating outside of this short-range regime function in the dynamic-range.

Neural determinants
At the spinal level, the alpha motor neurons integrate several afferent and efferent neural signals and generate the net neural drive recruiting muscle active elements, ultimately modulating force and stiffness. A clear efferent pathway contributing to the descending neural drive is the corticospinal tract, which transmits voluntary commands from the motor cortex down to the alpha motor neurons. An important afferent contribution derives from sensory elements in the MTU: the Golgi-tendon organs measure tendon tension, while the muscle spindles measure the muscular fiber stretch. Following perceived tension or stretch stimuli, the sensors send, respectively, inhibitory and excitatory contributions to the alpha motor neurons via monoand poly-synaptic pathways. Multiple other efferent pathways from the supraspinal centers play a role in the regulation of the joint stiffness [22].
Efferent and afferent regulation of net joint stiffness involves, respectively, feedforward mechanisms, e.g. co-contraction, affecting the neural drive prior to the application of a stretch, and feedback mechanisms, e.g. the reflexive response from the muscle spindles, affecting the neural drive as a response to the applied stretch [23]. Feedforward and feedback regulation mechanisms are in accordance with the performed task, the environment [24], and the properties of external perturbations [25]. Furthermore, the geometrical configuration of the joint and MTUs influences the modulation of the net joint stiffness by the CNS [23].

Toward a common terminology
The term intrinsic is used differently by the system identification and biomechanics communities. The system identification community often uses the term intrinsic impedance to refer to the combined mechanical properties of the muscle active elements regulated by the feedforward pathway (section 2.1.2) as well as by the joint constituent's passive components [26]. The biomechanics community, instead, often uses the term intrinsic to refer to material properties of musculoskeletal structures [27]. This creates potential for misunderstanding across disciplines about to what extent activation-dependent factors are considered to contribute to intrinsic properties.
To avoid misunderstandings across disciplines, we recommend the disuse of the term intrinsic, especially for collaborative studies. Instead, we propose the following classification: passive versus active joint stiffness. Passive joint stiffness represents the net joint stiffness contribution of passive skeletal and muscular elements. Active joint stiffness represents the net joint stiffness contribution from active muscular components. Additionally, a distinction between feedforward and feedback contributions to active joint stiffness may be made (section 2.1.2).

Key methodologies to obtain joint stiffness during human movement
Joint stiffness can be modeled via methodologies that assume different amounts of a priori knowledge. These methodologies can be placed on a spectrum that goes from white-box to black-box through different shades of gray (table 1). White-box models use purely mechanistic methodologies to describe joint stiffness [28]. A priori knowledge is used to obtain mathematical equations predicting the dynamics of all mechanisms contributing to joint stiffness. Black-box models use purely data-driven approaches (e.g. neural networks), where the functional relationships between recorded inputs and outputs are estimated without any reference to the underlying mechanisms [29]. Pure white-box models cannot be obtained, as a complete understanding of all details describing complex muscle geometrical properties and contraction mechanisms would be required. On the other hand, pure black-box models can be obtained, but lack a physical interpretation of the results, as they simply seek for mathematical relationships between the two given signals. To unravel the causal interplay between the neural and musculoskeletal mechanisms underlying joint stiffness, methodologies for joint stiffness estimation should ideally combine mechanistic and data-driven approaches to represent joint stiffness by a well-tuned gray-box model.
The system identification community typically employs dark-gray-box models to represent joint impedance, which is modeled using some level of a priori knowledge. Data-driven approaches are then used to retrieve the parameters of the model. Contrarily, the biomechanics community employs light-gray-box models to describe the elements and mechanisms contributing to joint stiffness. Sensor fusion is used to fill in parameter values that are not known to the model using information from multiple sensors, typically ex vivo measurements, electromyography (EMG), kinematic and dynamic data. Table 1 sorts representative articles along the black to white spectrum. Methodologies from both communities are classified based on the following aspects, in decreasing order of relevance: • Model architecture: The amount of distinct physiological pathways contributing to joint stiffness within the model. Within system identification methods, for instance, a model with a single, lumped pathway [30] is darker than a cascaded model with a dedicated reflexive pathway. Within musculoskeletal models, the amount of, and detail within, building blocks, i.e. muscle activation and contraction dynamics, determines the method's shade of gray. For instance, a model that uses experimental data and optimization routines to obtain muscle-tendon forces is darker than a model including explicit equations for muscle activation dynamics (e.g. a muscle twitch model), fiber contraction dynamics and series-elastic tendon architecture (e.g. inclusion of an elastic tendon instead of a stiff one). Moreover, a model that does not include muscle architecture features such as pennation angle dependencies and/or the kinematic relationship between changes in moment arms with changes in joint angles is darker than a model that includes these features. • System dynamics representation: The representation of the system dynamics within the model, either nonparametric or parametric. Nonparametric models are darker as system dynamics are represented by functions, infinite-dimensional in nature and with very limited assumptions on the underlying dynamics, e.g. impulse response function (IRF), or frequency response function (FRF) [13]. Parametric models are whiter, as analytical expressions with a predefined structure, order and assumptions are used, e.g. second-order inertia-damper-spring (IBK) model, or Hill-type muscle models.
Articles describing methodologies falling under the same classification criteria define a state-of-the-art cluster. Within a cluster, only representative articles are included in the table. State-of-the-art clusters are Table 1. Key state-of-the-art system identification and musculoskeletal modeling articles for joint stiffness estimation, sorted in a black to white spectrum. Articles with an asterisk * do not model joint stiffness explicitly.

Gray level
Model details Author Lumped nonpar. Ensemble TV IRF with nonparam. time variation Ludvig 2012 [31] Lumped joint impedance modeled by a TV IRF built from perturbation angle to torque response. Stiffness is defined as the integral of the IRF.
2nd order diff. Equation with nonparam. time variation Xu 1998 [32] Lumped joint impedance modeled by a second order inertia, spring, damper (IBK) difference equation from angle to torque. The IBK parameters of the model change nonparametrically over time.
2nd order diff. Equation with quasi-param. time variation van de Ruit 2020 [39] Lumped joint impedance modeled by a second order IBK difference equation from angle to torque. The statistical properties of the parameters in time are expressed by a kernel, which assumes smooth behavior.
2nd order diff. Equation with time-periodic time variation Guarin 2017 [40] Lumped joint impedance modeled by a Box-Jenkins model from angle to torque. Variation in dynamics are assumed time-periodic over the cycles. Time variation is represented by polynomials and continuity at the transition between cycles is not imposed.
Parallel of two TV IRFs with smooth, nonparam. time variation Guarin 2018 [30] Parallel of intrinsic joint impedance (TV IRF: from angle to torque) and reflexive joint impedance (TV IRF: from EMG to torque). Time variation is assumed to be smooth, described as cubic B-splines.
Parallel of high-pass filter and Hammerstein system with delay Jalaleddini 2017 [41] Parallel-cascade joint impedance model from angle to torque. The intrinsic pathway is described nonparametrically as an LTI IRF, whereas the reflex path is described by a cascade of a delay, differentiator and a Hammerstein model. A parametric state-space model describes the dynamic part of the latter.
Parallel of 2nd order diff. Equation and param. function with Ludvig 2007 [42] smooth, nonparam. time variation Parallel-cascade joint impedance model from angle to torque. The intrinsic pathway is modeled by a second order IBK difference equation, whereas the reflex path is modeled by a cascade of 40 ms delay, differentiator, half-wave rectifier and IRF. Low-pass filters enforce slow, smooth time variation.
SRS + stat. opt. Static optimization and muscle SRS model Hu 2011 [43] Endpoint stiffness is predicted from the active SRS of the modeled MTUs. Muscle forces are obtained via static optimization.
Inverse EMG-driven sim. and muscle SRS model Pfeifer 2012 [44] Active joint stiffness is predicted from the SRS of modeled MTUs. Co-contraction is captured via EMG, which is followed by static optimization to obtain muscle forces.
Inverse sim. and combination of SRS and passive force-length Zonnino 2019 [27] relationship Joint stiffness is estimated from active muscle SRS and passive force-length curve. The solution space of activations is studied for a given angle and torque target. Dependencies of moment arms and pennation angle with respect to joint angle are considered.

EMG-Hill
Forward EMG-driven sim. and muscle stiffness from force-length Sartori 2015 [45] relationship Joint stiffness is predicted from the length-dependent component of the modeled MTUs' force. Co-contraction is captured via EMG, which provides the excitations to perform forward simulations of a standard Hill-type muscle model.
A constitutive model is used to compute the passive and active tissue properties. Detailed musculoskeletal geometry is obtained from medical images. incorporated in the table if they contain at least one journal article and one article where the proposed methodology was applied to experimental data. For the system identification community, only methodologies capable of tracking the variation of joint stiffness as continuous functions of time are considered, i.e. methods for time-varying (TV) and parameter-varying models. Such methodologies are therefore applicable for the identification of joint stiffness during human movement. Table 1 provides an overview of representative methodologies for capturing joint stiffness over time. This section provides the authors' selected methodologies in the context of system identification and musculoskeletal modeling separately. Section 4 outlines how the selected methodologies can be combined for achieving the proposed unified framework.

Selected system identification methodology
System identification methodologies require a trade-off between the a priori knowledge to be included in the model and the amount of data required. For the proposed framework, we choose a system identification methodology resulting in a slowly TV lumped nonparametric representation of joint impedance. We select a dark-gray box methodology to limit the number of assumptions made from the system identification side and provide a reliable estimate of joint stiffness. The drawback is that an ensemble of realizations is required, increasing the experimental time.
Computationally, the short data segment (SDS) method (such as Ludvig et al [31]) is a desirable option to estimate the nonparametric model dynamics, which slowly vary with time. The SDS method has been widely used in the literature to estimate joint impedance and stiffness from human experimental data [33][34][35][36]. An overview of pros and cons relative to the SDS method is depicted in figure 3.

Theoretical level description
Joint impedance can be represented by a linear time-varying (LTV) model, capturing the variations of joint impedance with changes in the operational point [37,38]. The TV FRF H(t, jω) is a common TV model to represent joint impedance nonparametrically. The TV FRF linearly relates the imposed angular displacement perturbation signal, ϕ(t), to the torque response signal, τ ϕ (t) [38]. The relation is expressed in the frequency domain. The TV FRF is two-dimensional: it considers both the system dynamics, as dependence on the frequency jω, and time variation, as dependence on time t. The relationship between ϕ(t) and τ ϕ (t) is described as: where F −1 and F are the inverse and direct Fourier transform operators, respectively. The value of H(t, jω) for ω = 0 is the static component of the TV FRF and corresponds to joint stiffness as a function of time [47]. As a special case, by rewriting joint impedance parametrically with a commonly used IBK model in the linear time-invariant case, the following relation is obtained: from which the static gain for ω = 0 corresponds to the stiffness K. In the TV case, a similar comparison can be done between H(t, jω) and a parametric TV IBK model, in this case obtaining K(t): stiffness as a function of time.
In equation (1), a single, lumped pathway is used to model the dynamics through the TV FRF. The main constraints of modeling joint impedance by a TV FRF are that the measurements should be sufficiently long to extract joint stiffness and the fact that the input and output are related linearly. On the other hand, the main advantage is that, since no direct distinctions between the mechanisms contributing to joint impedance are made in the model, there are limited a priori assumptions on the determinants of joint stiffness. Furthermore, the model structure is valid to represent joint impedance under isometric and dynamic experimental conditions [37].

Computational level description
To compute joint stiffness, we propose to estimate the TV FRF in equation (1) via the SDS computational method. The SDS method estimates the TV FRF via a correlation-based technique, as described in [37]. The correlation-based technique is analogous to the singular-value-decomposition based linear-least squares proposed in [38], although it is computationally more efficient [31]. Since a nonparametric structure is used to represent the time variation, a data ensemble of realizations with the same TV conditions is required to perform the estimation. This corresponds to assuming that joint impedance is linear time-periodic. To reduce the number of realizations required for the SDS method, the variation of the parameters in time is assumed to be locally time-invariant [31]. Local time-invariance entails that, in a short segment of length N, joint impedance can be considered constant. Such assumption significantly reduces the number of ensemble measurements required, although it is only valid for slow TV behaviors.
In equation (1), the system dynamics are modeled by means of the TV FRF, while in the SDS the TV IRF is used instead. Since the TV IRF is the time-domain equivalent of the TV FRF [39], the same considerations hold.

Experimental level description
The TV FRF model is estimated from an ensemble of input-output realizations obtained under the same experimental conditions. The data ensemble can be obtained experimentally by using a robotic manipulator to apply small amplitude perturbations to a joint [34]. The joint angle can be recorded using an encoder integrated into the manipulator's motor, whereas torque could be measured with a load cell placed in the proximity of the motor. The experimental protocol can include isometric or dynamic conditions [34]. In isometric conditions, the joint interacts with a stiff position-controlled manipulator, which imposes angular perturbations around a fixed joint position. To induce time variation in joint stiffness, subjects are instructed to change the voluntary torque by tracking a predefined trajectory. In dynamic conditions, the joint interacts with a compliant torque-controlled manipulator imposing torque perturbations, while subjects vary the joint angle. To disentangle joint impedance from the compliant manipulator environment, the SDS algorithm as presented in [31] can be adapted to be applicable to dynamic experimental conditions, following the formulations in [34]. Currently, most experiments using the SDS method have been performed with subjects holding a fixed posture (e.g. sitting) [34,35]. Moreover, perturbations were applied to a single degree of freedom (DoF) of a single joint.

Selected musculoskeletal modeling methodology
We propose high-density electromyography (HD-EMG)-driven musculoskeletal modeling techniques [48] for the estimation of muscle, tendon, and joint forces as a function of measured HD-EMG signals and joint angles [49][50][51]. This methodology was recently extended to estimate muscle, tendon, MTU-equivalent, as well as joint stiffness as proposed by Sartori et al [45] and generalized to account for pennation angles in [52]. In this section we provide a further refined formulation. An overview of pros and cons of the proposed method is depicted in figure 3.

Theoretical level description
Net joint stiffness results from the stiffness of joint's constitutive elements (section 2.1). A fundamental step in any mechanistic modeling technique is to choose which biological structures to model a priori. In the remainder of this sub-section, we consider MTUs to be the main contributors to joint stiffness, and therefore the key biological structure to be modeled [16]. Nevertheless, modeling of passive skeletal elements such as ligaments should also be considered for joint stiffness estimation [53].
MTUs could be modeled using the standard Hill-type muscle model as a basis for additional development. Hill-type muscle models are preferred over the more detailed Huxley-type [54,55] (previous research showed both model types yield similar results [56]) or finite element models [57,58]. This is because they are computationally efficient, thereby making them a good candidate for achieving computational tractability and physiological plausibility simultaneously, something crucial for enabling real-time robotics and rehabilitation applications [59,60]. The Hill-type muscle model often used in musculoskeletal simulations comprises two distinct blocks for activation and contraction dynamics [49][50][51].
In the activation dynamics block, muscle activation is directly extracted from experimental EMG recordings [17,50]. An advanced activation block would combine HD-EMG recordings and decomposition techniques to extract muscle motor units spike trains to drive the model's MTUs [48,61,62]. The main advantage of decomposing motor units spike trains is the ability to extract and preserve the natural frequency content of motor unit firings across any type of contraction. Consequently, slow type motor units could be linked to a slow twitch model and fast type motor units to a fast twitch model [63], enabling a more robust estimation of force, and therefore stiffness, across a large repertoire of tasks. This is different to typically used bipolar EMG recordings in combination with linear filtering with an a priori chosen cutoff frequency. A priori chosen cutoff frequencies would be suitable only for movements involving muscle contractions at certain (matching) frequencies (i.e. for which ad hoc filtering cutoff frequencies should be chosen upon) [Chapter 3 of [19]]. Instead, using HD-EMG and decomposition techniques it would be possible to selectively remove noise from the EMG, thus recovering the motor unit pool's net neural drive, i.e. emerging from efferent and afferent pathways [64,65].
The contraction dynamics block should include the following features: • A contractile element (CE) representing the muscle's activation dependent force generating mechanisms. These mechanisms should comprise both the well established active force-length and force-velocity relationships [66], and the less understood history-dependent muscle properties, e.g. lengthening-induced force enhancement and shortening-induced force depression [67,68]. Modeling history-dependent muscle properties might have implications in daily-life movements [69] and in the context of wearable robotics, where exoskeletal robots may alter the natural stretch-shortening cycle of a muscle. With the proposed CE, isometric, concentric, and eccentric muscle contractions would be well characterized. • A parallel element (PE) representing the passive tissue properties of the muscle. Ideally, a subject-and muscle-specific passive force-length relationship should be constructed [70]. • A SRS element accounting for instantaneous increase in muscle force due to short stretches of muscle fibers [71]. This feature is paramount to estimate muscle stiffness using perturbation-based data. • A series elastic element (SEE) representing the tendon dynamics. A recent study assessed the use of surrogate tendon models, e.g. using the finite element method, to extract subject specific tendon forces and strains [72].
The aforementioned model would allow for the computation of muscle, tendon, as well as the equivalent MTU force using EMG and joint angles as inputs. In contrast to inverse dynamics (ID) formulations, our proposed method would not require availability of experimental ground reaction forces (GRFs) once the model is calibrated (i.e. GRFs would be needed to derive reference joint torques via ID for calibrating the HD-EMG-driven model's MTU parameters once per subject). This would have benefits in robotics and rehabilitation applications, where measurement of GRFs would no longer be required, thereby simplifying experimental setups and enhancing portability. To obtain stiffness, we propose the following equations.
The stiffness of a muscle along the MTU's line of action, K m eq , could be expressed as follows (adapted from [73]): where F m eq and l m eq are the force and length, respectively, of the muscle fiber along the direction of the tendon's line of action, F m and l m are the force and length, respectively, of the muscle fiber along its axis. The variables a,l m ,v m and α represent the muscle fiber's activation, normalized length, normalized velocity, and pennation angle, respectively.
Tendon stiffness (K t ) is computed as: where F t (ϵ) is the non-linear tendon force as a function of tendon strain ϵ [ϵ = (l t − l s )/l s , where l s is the tendon slack length and l t is the tendon length]. The total MTU stiffness is the series arrangement of K m eq and K t : K mtu = Then, the net joint stiffness (K J ) can be computed from the stiffness of its constituting MTUs: where F mtu i , r i and K mtu i represent the force, moment arm and stiffness, respectively, of the ith MTU spanning the joint, and θ is the joint angle. The first term in equation (5) captures the stiffness of all the MTUs projected at the joint level via their moment arms. The second term in equation (5) represents joint stiffness contributions resulting from angle-dependent changes in moment arms. Figure 2 displays the different parameters of interest of the musculoskeletal model.

Computational level description
OpenSim, an open source software to perform multibody dynamics simulations [74,75], contains a wide range of generic musculoskeletal models obtained from cadavers. The geometric and inertial properties of MTU length, l mtu , could be obtained using OpenSim or directly measured using ultrasonography. A computational viable way to compute moment arms could be that of fitting cubic splines to MTU length nominal data across the joints' range of motion and using the tendon excursion method (based on the principle of virtual work) to estimate moment arms [76]. This would be computationally efficient because the evaluation of the spline is in closed-form. Additionally, the piece-wise polynomial fit would allow for simple computation of ∂r i /∂θ (used in equation (5)) and jacobians. However, this approach could be heavy on computer memory as the splines' coefficients (for each MTU and DoF) need to be stored, i.e. with memory storage requirements increasing with the number of DoFs spanned by the MTUs.
The open source 'Calibrated Electromyography-Informed NeuroMusculoSkeletal' modeling toolbox (CEINMS) [77] is chosen as a basis to implement the proposed MTU model to obtain muscle and MTU forces (used in equations (3) and (5)) from EMGs and MTU kinematics in a forward dynamics fashion. A muscle formulation based on a Hill-type muscle model comprises only one differential equation per modeled MTU that can be efficiently computed in real time using a root solver [78]. In contrast, Huxley-type and finite element models would require the simultaneous integration of multiple differential equations, which nowadays is not feasible in real time computations [46,54,79,80].
Fiber activation is obtained from experimental HD-EMG recordings, and fiber length and velocity are estimated from the joint kinematics. The muscle force is projected to the tendon taking the pennation angle into account: F mtu = F t = F m cos α. The pennation angle, α is computed assuming that the width of the Hill-type muscle is constant. Computational speed is assured because all aforementioned steps are solved via closed-form equations that do not require iterative methods. Therefore, CEINMS can estimate muscle and joint forces and stiffness in open-loop, i.e. reference joint torques to solve an optimization problem are not necessary [51]. However, a prior calibration stage requires experimental joint torques to tune MTU parameters to create a subject-specific musculoskeletal model. This calibration could be achieved via optimization routines [77] or machine learning algorithms [81]. Global optimizers such as simulated annealing [82,83] or covariance matrix adaptation [84] are preferred choices due to their ability to avoid local minima, but they can be time consuming. If computational speed is a priority, gradient-based optimizers, such as the interior point optimizer [85] should be considered. Machine learning algorithms obtain promising results when it comes to model personalization [81], but they require large experimental data sets to assure true learning.
The computation of ∂F m /∂l m is required in equation (3). This partial derivative could be obtained via the a priori computation of a multidimensional spline that maps the whole solution space of muscle force as a function of activation and normalized velocity and length, allowing for a fast computation of the partial derivative with respect to fiber length, i.e. muscle stiffness. This makes the method feasible for real-time computation of joint stiffness. Nevertheless, it relies on the computer's memory as many coefficients must be stored. Another option to compute these derivatives are automatic differentiation methods [86].

Experimental level description
Kinematic, kinetic and EMG data are required to obtain joint stiffness using our proposed HD-EMG-driven model-based approach [45,87]. To obtain kinematic data, i.e. joint angles, MTU lengths, and moment arms, of multiple DoFs and MTUs simultaneously in a computationally and experimentally viable way, we recommend the use of motion capture analyses, inertial measurement units, or angle encoders from orthoses, in combination with multibody dynamics software such as OpenSim. Kinetic data, i.e. joint torques, are needed as reference data to tune the model's parameters in the calibration phase. Reference joint torques could be obtained via ID and GRFs, or via a load cell in an orthosis or insole. Once the model's parameters are calibrated, the reference joint torque is no longer required as the model runs in open-loop using EMGs and joint angles as inputs. The activity of the muscles of interest could be recorded using surface HD-EMG electrodes. Further decomposition of the acquired EMGs into motor unit spike trains would allow measuring activity of muscle firing at different frequencies depending on the task [48].

A window of opportunity for a combined cross-disciplinary method
Sections 3.1 and 3.2 briefly described the fundamentals of the selected system identification method and musculoskeletal modeling technique, respectively. Each of them has strengths and limitations that are complementary. Figure 3 summarizes these advantages and drawbacks. This compatibility across communities provides an excellent opportunity to merge these methods in a combined framework to estimate, and consequently better understand, joint stiffness.

Roadmap for a unified framework
In section 3 we described the selected system identification method and the proposed musculoskeletal model to estimate joint stiffness during movement. In this section we present our vision for a unified framework that combines both methodologies. From a theoretical level, we discuss the mathematical languages used in system identification and model-based methods as well as how these can be put in relation to each other. From a computational level, we propose a novel closed-loop algorithm that combines both methods to extract multiscale stiffness estimates, i.e. from joint to muscle, to tendon scales. From an experimental level, we discuss the experimental protocols to be used to generate the movement data needed by the proposed closed-loop framework and the experimental requirements needed to assure the validity of both methods within the same experiment. . Schematic diagram of the HD-EMG-driven musculoskeletal model and LTV nonparametric system identification approaches for the estimation of (ankle) joint stiffness. (Middle) The neuromechanical states affecting joint impedance (i.e. the operational point) include the joint angle, the properties of the applied angular displacement, the joint torque and the HD-EMGs (characterizing the joint activation). Joint impedance relates the torque response τ ϕ (t) to an externally imposed angular displacement perturbation ϕ(t). (Bottom) LTV nonparametric system identification analyses the relationship between ϕ(t) and τ ϕ (t). An impedance model is fit to the signals by minimizing the error between τ ϕ (t) and the predicted torque responseτ ϕ . An estimate of net joint stiffnessK J (t) is obtained from the model. (Top) The HD-EMG-driven musculoskeletal model predicts the net joint stiffnessK J (t) and joint torque τ (t) from the measured joint angle and HD-EMG, using models of the MTU activation, kinematics and dynamics and of the resultant joint dynamics. (Dashed lines) Proposed closed-loop tracking approach to adjust the parameters of the HD-EMG-driven musculoskeletal model by comparing the stiffness estimates between the two approaches, and the predicted torque to the measured torque. The closed-loop tracking can be used for calibration of the MTU dynamics model or for (minimal) adjustment of the HD-EMG signals. Figure 4 illustrates the bottom-up approach followed in biomechanics as well as the top-down approach followed by system identification communities. In musculoskeletal modeling, sensor fusion (e.g. fusing EMGs and joint angles) is used to drive a computational model that synthesizes the kinematics and dynamics of the MTU and joint [88]. In this, joint stiffness is decoded from the stiffness of the model's constituting elements [27,45]. In the SDS system identification method, a model of joint impedance is built using perturbation-based measurements and joint stiffness is extracted from joint impedance. It is non-trivial to compare joint stiffness estimations from two completely different methodologies, and in order to use them in a single framework it is paramount that both techniques are estimating the same underlying neuromechanical property.

Theoretical level description
The mathematics used in both techniques is typically based on different formulations. On the one hand, the SDS system identification method obtains a TV FRF from experimental data, representing the transfer function in response to an angular displacement perturbation, i.e. joint impedance. The static component of the transfer function is computed by taking the DC gain of the TV FRF for ω = 0, representing the component of the transfer function that does not depend on velocity or higher-order terms. Since the TV FRF in equation (1) is computed between position and torque signals, the static component corresponds to joint stiffness, as illustrated in section 3.1.
On the other hand, our proposed musculoskeletal model obtains joint stiffness from the projected stiffness of its constituting elements equation (5). The stiffness is computed by taking the partial derivative of force with respect to length (i.e. physical notion of stiffness) as described previously, equations (3) and (4). In both cases, stiffness represents the change in force/torque due to an infinitesimal change in displacement.
However, it is computed via different functions. In the SDS, stiffness is computed from the TV FRF, a function relating dynamically angular perturbation and torque response, while in the HD-EMG-driven musculoskeletal modeling technique, it is computed from a multidimensional mapping function relating muscle length, velocity, pennation angle and activation to muscle force (multidimensional spline), equation (3).
Additionally, each method operates at a different level. The SDS method uses experimental data collected at the joint level, i.e. joint torque and angle. Hence, the TV FRF and stiffness estimates are directly at the joint level. Musculoskeletal models compute the stiffness at muscle and tendon levels, and then project the stiffness to higher levels, i.e. MTU and joint. Assumptions on the number of modeled elements and approximations in the computation of variables such as the moment arms are required in this multiscale definition of joint stiffness. These dissimilarities regarding mathematical formulations and scale suggest that special attention must be paid when comparing joint stiffness resulting from both methodologies.
Moreover, multiple mechanisms and structures influence joint stiffness, and the properties of the applied perturbation signal affect the mechanisms that are excited and that produce an observable change in the response torque. For example, the properties of the perturbation signal (e.g. acceleration, velocity, duration), can affect the reflexive responses that are elicited and can induce a co-contraction and/or predictive response, i.e. the subjects performing the experiment might modulate their neuromechanical states to counteract the applied perturbations. Furthermore, the amplitude and speed of the perturbation signal influence whether the muscular stretch response is in the short-or in the dynamic-range (section 2.1.1). In this context, a musculoskeletal model should be set up in relation to the design choices made at the system identification level. That is, if the system identification approach (algorithm and experimental protocol) has been designed to probe joint stiffness response in the short-range, then the associated musculoskeletal model should be equipped with SRS modules in parallel with muscle CEs (section 3.2.1).
The most generalizable musculoskeletal (light-gray-box) model would contain analytical formulations of each 'perturbation-dependent mechanism' (i.e. reflexive stiffness, short-and dynamic-range stiffness, etc) and combine formulations to match the neuromuscular phenomena triggered by the perturbation signal. Modeling the different neuromuscular phenomena is not required for dark-gray-box system identification methods because they are data-driven methods, i.e. the data already includes the contribution of all excited neuromuscular phenomena (e.g. SRS, reflexive active stiffness, etc). In this scenario, in order to achieve consistency between stiffness estimates from system identification and musculoskeletal modeling methods: (a) The musculoskeletal model structure should a priori take into account what neuromuscular phenomena are triggered by the perturbation signal and explicitly model them. (b) The stiffness estimate comparison between model-based and system identification methods can only be done using data from the same perturbation-driven experimental conditions.

Computational level description
Here we propose a coupled closed-loop computational framework that overcomes shortcomings and highlights advantages of each method (figure 3). A closed-loop formulation in which data-driven joint stiffness estimates obtained via system identification are tracked by model-based predictions of joint stiffness is proposed and illustrated in figure 4. The proposed system identification method is applied to perturbation-based data to obtain an estimate of joint stiffness over time. The estimate, obtained from experimental data that can be validated, provides a reference value to be tracked by the model-based prediction [89]. Similarly to the torque-tracking closed-loop EMG-driven biomechanical modeling paradigm presented by Sartori et al [89], we propose a modeling paradigm based on simultaneous closed-loop tracking of stiffness and torque.
The purpose of the presented closed-loop control scheme (figure 4) is twofold: (1) to refine TV model inputs such as HD-EMGs and (2) to tune, during the calibration phase, physiological MTU parameters that directly depend on joint stiffness (section 3.2.2, i.e. slope of muscle/tendon passive force-length curves). HD-EMG-derived muscle activations can be minimally adjusted to minimize the combined joint stiffness and torque tracking error. Moreover, muscle activations of MTUs included in the model from which it is not feasible to acquire experimental measurements e.g. deep muscles, can be synthesized following this formulation [89]. During the calibration phase (section 3.2, Computational level description), minimizing the tracking error at the stiffness level, in addition to the current torque tracking formulation, might enable a more physiologically consistent identification of MTU parameters.
Whether an estimated parameter can be used to infer a physiological parameter could be checked by determining the identifiability of physiological parameters set forth. Such identifiability can be extracted from a sensitivity analysis, performed either analytically or numerically [90,91], to compute the covariance of the model parameters. The sensitivity analysis can provide information on the confidence in the model Figure 5. Example of how the proposed paradigm would work in a simple experimental protocol involving sinusoidal ankle plantar-and dorsiflexions. Experimentally recorded signals from a perturbation-based protocol, i.e. joint torque and joint angle, are used to obtain reference joint stiffness values via the SDS method. Joint torque, angle and HD-EMG signals are used to obtain joint stiffness predictions via the HD-EMG-driven musculoskeletal model. These predictions are iteratively optimized by minimally adjusting HD-EMG-derived muscle activations. Once the tracking error is minimized, MTU, muscle, and tendon contributions to the net joint stiffness can be obtained thus enabling a multiscale understanding of joint stiffness. and show whether the model parameters can be determined independently from each other or interdependences between parameters must be taken into account.
The proposed paradigm will provide reliable joint stiffness estimates while opening a window toward smaller scales estimates, e.g. muscles and tendons ( figure 5). Moreover, having a musculoskeletal model that tracks not only joint torque but also joint stiffness would further reduce the solution space of muscle forces, advancing the biomechanics community toward solving the muscle redundancy problem [92,93]. Finally, the ability to track a reliable estimate of joint stiffness error might push forward the field of musculoskeletal optimizations via the development of more complete objective functions that take into account joint stability in addition to traditional optimization criteria, e.g. energy expenditure [94].
In addition to tracking joint stiffness and torque, the proposed framework could be further extended to simultaneously track other biomechanical variables as well. For instance, ultrasonography recordings could be used in real time to measure, and consequently track, muscle fascicle kinematics in vivo [95].

Experimental level description
System identification and HD-EMG-driven musculoskeletal modeling methods possess different experimental requirements. In this section we present how a generic experimental setup that would allow the simultaneous application of both techniques.
The most simple experimental setup that allows studying joint stiffness for a single DoF contains the following two elements: a specialized manipulator capable of applying perturbations while measuring the joint's angle and torque; and EMG measurements (bipolar or high-density) of the main muscles spanning the joint of interest. To account for muscles spanning multiple joints, e.g. gastrocnemius medialis, more detailed biomechanical models require joint angles of adjacent joints as well. This can be achieved by preferably using a motion capture system or inertial measurement units. In the absence of these devices, an electro-goniometer could be used to measure static joint angles. An example of this setup can be seen in [87].
A variety of tasks and experimental conditions can be studied with this setup, from static to dynamic tasks. Depending on the type of perturbation that is applied, different physiological mechanisms can be excited, e.g. reflexes or SRS. The collected data will always allow obtaining joint stiffness via both techniques. From the system identification community, the selected method uses a LTV model to represent joint impedance. To satisfy this assumption, small amplitude perturbations should be applied.
Once the proposed paradigm is validated in a controlled setup and an isolated DoF, more complex experiments involving multiple DoFs and muscle contraction mechanisms, e.g. eccentric contractions, can be designed.

Open challenges and future steps
In this article, we presented a roadmap to combine data-driven musculoskeletal modeling and system identification techniques for the estimation of joint stiffness, a neuromechanical property that dictates how humans interact with the environment. We highlighted the potential benefits of combining the efforts of two communities that rarely communicate and cross-fertilize, despite being so complementary. Whilst system identification provides reliable stiffness estimates at the joint level, musculoskeletal modeling provides stiffness predictions at muscle-tendon levels too. We therefore proposed key developments in selected musculoskeletal modeling and system identification methodologies as well as an integrative closed-loop formulation in which joint stiffness profiles estimated via system identification are tracked by an optimized HD-EMG-driven musculoskeletal model.
The proposed framework should help understand the potentials of such a cross-disciplinary collaboration and serve as a launching step for merging the theoretical, computational and experimental divergences across communities. Previous studies already extended CEs to include modules accounting for history-dependent muscle properties (figure 2), e.g. stretch-induced force enhancement and shortening-induced force depression [67,68], and therefore constitute a solid foundation for future implementations. Moreover, an extended Hill-type muscle model including an SRS module already exists [71], and the effect of SRS on joint kinematics and kinetics was evaluated. Nevertheless, to the best of our knowledge, the effect on joint stiffness of the inclusion of an SRS module in an extended muscle model has not been explicitly investigated yet. Methods to obtain the net neural drive of the motor neuron pool to drive a musculoskeletal model have also been implemented [48,62,91,96]. Our proposed musculoskeletal model would include features of all aforementioned methods in a single formulation and be combined with the suggested system identification method.
The proposed musculoskeletal model advances the current state of the art in joint stiffness estimation via bottom-up techniques. However, further research is required to address current shortcomings in our physiological understanding of the mechanisms regulating joint stiffness and experimental limitations to properly record data to estimate joint stiffness. For instance, the model comprises a set of parameters that need to be calibrated using experimental data. Whether the underlying interpretation of the model parameters is correct depends on the correctness of the model's physiological representation. Given that the state-of-the-art has still limited understanding of the underlying physiological mechanisms regulating joint stiffness (section 2.1), it is not possible to ensure convergence toward the true physiological values. Nevertheless, by validating the calibrated model against unseen reference torque and stiffness estimates across a variety of tasks, we gain confidence in the model.
As is typical for system identification of complex systems, a dark-gray-box method was selected to limit the underlying assumptions and provide a reliable comparison method. With the advancing of the knowledge on a multiscale joint stiffness model, whiter models can be used, explicitly modeling the changes of joint stiffness with respect to the operational point.
For simplicity, the proposed closed-loop approach was here illustrated for a single DoF ( figure 4). Nevertheless, the presented method could be potentially extended to multiple DoFs if experimental limitations (regarding, for instance, simultaneous perturbation of multiple joints) are overcome. This would enable a multiscale understanding of joint stiffness across multiple DoFs during functional tasks such as locomotion.
Finally, several future paths could be traced from the proposed framework. The concepts explored in this article could be transferred to other communities or methodologies. System identification techniques have already been applied in combination with ultrasonography [97]. Moreover, despite being still premature, elastography is showing promising results to perform system identification of the stiffness at a muscular level [98].

Conclusion
Reliable estimation of joint stiffness during human movement requires establishing a novel cross-disciplinary framework. We highlighted the strengths and weaknesses of the chosen system identification and musculoskeletal modeling methods, i.e. SDS method and HD-EMG-driven musculoskeletal modeling, respectively, and consequently detailed a roadmap to combine these techniques from different scientific communities. Even though the focus was on joint stiffness across system identification and biomechanics communities, the proposed roadmap may serve to favor cross-fertilization across scientific communities spanning from medical imaging to robotics, from numerical modeling to sensor fusion and system identification.

Data availability statement
No new data were created or analysed in this study.