Model-Based Estimation of Ankle Joint Stiffness

We address the estimation of biomechanical parameters with wearable measurement technologies. In particular, we focus on the estimation of sagittal plane ankle joint stiffness in dorsiflexion/plantar flexion. For this estimation, a novel nonlinear biomechanical model of the lower leg was formulated that is driven by electromyographic signals. The model incorporates a two-dimensional kinematic description in the sagittal plane for the calculation of muscle lever arms and torques. To reduce estimation errors due to model uncertainties, a filtering algorithm is necessary that employs segmental orientation sensor measurements. Because of the model’s inherent nonlinearities and nonsmooth dynamics, a square-root cubature Kalman filter was developed. The performance of the novel estimation approach was evaluated in silico and in an experimental procedure. The experimental study was conducted with body-worn sensors and a test-bench that was specifically designed to obtain reference angle and torque measurements for a single joint. Results show that the filter is able to reconstruct joint angle positions, velocities and torque, as well as, joint stiffness during experimental test bench movements.


Introduction
An accurate real-time estimation of human joint stiffness would be beneficial in many situations. For example, active prosthetics, orthotics, or exoskeletons could use information on quantitative joint stiffness to adjust their control strategies [1]. Currently, such systems are increasingly fitted with physically adjustable elements [2]. Therefore, the adjustment of, for example, a physically controllable interaction stiffness in robotic support is a possible area of application. Besides these movement support applications, time-resolved knowledge of human joint stiffness can be useful in motion analysis or as a detection system for pathologic movement states, such as spasticity [3].
The stiffness of human lower extremity joints changes naturally during locomotion. This phenomenon is typically called stiffness modulation and serves as an involuntary adaptation of the neuromusculoskeletal system to handle contact situations in a highly uncertain environment. Stiffness modulation is suggested to increase gait stability, or absorb shocks [4]. Although it can be argued that stiffness is the major component that is changed during human locomotion, other passive elements (such as viscous damping) also exist. Together with inertia, joint dynamic properties can be generalised to the concept of mechanical impedance. A mechanical impedance can be regarded as a frequency dependent and (for most human joints) nonlinear mechanical resistance to excitation of a force (torque) or a velocity (angular rate). Since its recognition in locomotion in the 1980s [5], mechanical impedance has been well researched and is suggested to play a major role in the stabilisation of unstable walking with ankle and hip dynamic models. The human ankle joint consists of an upper (articulatio talocruralis; TC) and a lower (articulatio subtalaris and articulatio talocalcaneonavicularis) part. This study focusses on the modelling of dorsiflexion and plantar flexion, which is mainly achieved by the TC [29][30][31]. Although the TC is generally considered as a one-degree-of-freedom joint disagreement exists about the rotation axis/axes. Early studies indicated two discrete axes of rotation [32,33], one at dorsiflexion and one at plantar flexion. However, a subsequent study identified a single, fixed axis of rotation [34]. More recent research considered a continuously shifting axis of rotation with foot movement [29,30,35]. Nonetheless, some studies support the approach of a single fixed axis [31,36]. Moreover, the modelling approach with a single hinge joint was successfully employed in several studies that yielded plausible results [37][38][39][40][41]. Consequently, in the present study we model the TC as a simple hinge joint.

Ankle Kinematic Model
The local foot coordinate system is set-up with respect to the knee reference system, where the latter is (later on) related to the global reference system located at the hip. The knee reference coordinate system is located at the contact points between condyles and tibia plateau at a knee angle of 0 • . Thus, the y-axis y knee points along the tibia, whereas x knee is anterograde at upright standing. z knee is chosen so as to span a right-handed coordinate system. A local coordinate system is defined for the rotation of the foot, which is body-fixed to the foot. However, this allows to rotate the foot with respect to the global coordinate system. Figure 1 shows the local foot coordinate system, which rotates about the z-axis z ankle . The inclination and orientation of z ankle were determined according to the results of [34].  Left: View at the coordinate system from the front; α TC tla is the angle between tibia axis and rotation axis. Right: View at the coordinate system from above; α TC f ml is the angle between mid-line of the foot and the rotation axis.
The origin of the foot coordinate system was placed between both malleoli. The local x-axis x f oot is perpendicular to z f oot and parallel to the foot sole. y f oot is aligned so as to span a right-handed coordinate system. The ankle joint angle ϕ ankle is defined by a rotation of the foot coordinate system around z f oot with respect to the global y-axis. Furthermore, the zero position ϕ ankle = 0 • is defined at upright standing (y knee is perpendicular to x ankle ). The angular range of the ankle joint was defined to lie on the interval −60 • < ϕ ankle < 20 • . The resulting transformation of a point p f oot in the local coordinate system of the foot is defined in homogeneous coordinates with help of the homogeneous transformation matrix T(ϕ ankle ) To set-up rotation matrices, x trans and y trans describe the origin of the foot coordinate system in the knee coordinate system. Moreover, the rotation angles are defined as β = α TC f ml − 90 • and γ = 90 • − α TC tla . For actuation of the ankle joint, kinematic considerations are made for three main muscles: the musculus tibialis anterior (TA) as a primary actuator for dorsiflexion, as well as the musculus soleus (SOL) and the musculus gastrocnemius (GAS) as the primary actuators for plantar flexion [42]. The reasons for this decision were: a) to have reliable sEMG signals in the filter, and b) to limit the dimension of the resulting state-space model (each muscle increases the filter order by 2.) Figure 2 is an overview of the muscle kinematics used to calculate the lever arms. Ankle joint model (green: TA, blue: SOL, red: GAS). Dotted lines illustrate the lever arms; the star indicates the variable centre of rotation that is obtained by projection into the sagittal plane. TA, musculus tibialis anterior; SOL, musculus soleus; GAS, musculus gastrocnemius.
The TA originates from the lateral condyle of the tibia and sets on along the inside of the foot. Figure 2 shows that the TA is guided by two medially aligned ligaments (retinaculum musculorum extensorum superius and inferius). The GAS (Note that the GAS model was scaled accordingly, to account for the medial and lateral part simultaneously) originates from its two heads at the condyles of the femur where, in case of the SOL, the origin is the posterior-proximal end of the tibia and fibula. Both muscles are joined in the Achilles tendon which inserts into the posterior surface of the heel bone (calcaneus). As described in [28], the muscles are approximated by line segments from origin to fixation point. The fixation of the TA by means of the retinacula is realised by two nodal points for determination of the effective origin and point of action. These two points were taken from the literature for the shank and for the foot [43,44], respectively ( Figure 2). From the position of the TC rotation axis and the coordinates of the nodal point, the length of the MTC l MT,i and the lever arms r ma,i can be calculated. The length l MT,i is a result of the sum of the single line segment. The lever arms are defined as the distance of the rotation axis to the effective direction of the muscle force. The effective direction lies on the straight line denote global coordinates of the projection of the origin and the contact point of the muscle into the sagittal plane. Hence, the lever arm of the muscle i results in In Equation (2),p global IC represents the global coordinates of the current rotation axis. Due to the projection into the sagittal plane, a resulting movement of rotation axis can be observed and is visualised in Figure 3a,b.p  Figure 3. Determination of the current centre-of-rotation in the sagittal plane. (a) Determination of the current centre-of-rotation by a small angle change ∆ϕ ankle .p IC results from the intersection of the perpendicular lines bisecting straight lines connecting p I,i (ϕ ankle ) and p I,i (ϕ ankle + ∆ϕ ankle ). (b) Change of current centre-of-rotation with ϕ ankle with projection into the sagittal plane.
Note that if the muscle has an effective origin or contact point, these are used to determine the effective direction. This is, for example, the case for the TA, where the effective origin is in the shank and the effective contact point is in the foot. In addition, the GAS has an effective origin at the condyles of the femur given the corresponding knee alignment. Moreover, the GAS is a biarticulate muscle acting on knee and ankle joint at the same time and, as such, the corresponding variables depend on the orientation of the ankle and knee joint angles.
Line segment lengths and lever arms of the modelled kinematics were calculated offline over the whole range of ankle angles ϕ ankle and knee angles ϕ knee at increments of ∆ϕ angle = 1 • and approximated by fourth-order polynomials. The muscle contraction speed v M can then be calculated from angular velocity and time-derivative of l MT (ϕ ankle ). For TA and SOL, the time-derivative of with muscle length l M and pennation angle α p . As for the GAS we obtain Note that in Equations (3) and (4) it is assumed that the tendon length l T is constant. However, whereas taking the Achilles tendon strain into account is important to determine joint stiffness, computational simplifications are necessary to make the model treatable.

Hip Kinematic Model
The hip joint is typically modelled as a 3 DOF ball joint [45] or a 1 DOF hinge joint in the sagittal plane [43]. For the purpose of the present model, only the influence of the sagittal plane orientation of the hip is considered. Hence, the hip joint is modelled as a hinge joint. The MTC length of the biarticular muscles that act on the knee and hip joint is influenced by the orientation of the hip joint. For further simplification, the influence of hip abduction/adduction and outer/inner rotation is considered to be small as compared to hip flexion/extension and is subsequently neglected. Biarticular muscle groups that are considered in the knee model presented by Misgeld et al. [28] are the musculus quadriceps femoris (QF) and the ischiocrural muscles (HAMS). The lengths of their corresponding MTCs l MT,QF and l MT,H AMS are influenced by the hip angle ϕ hip . The hip angle ϕ hip is defined as the angle between the pelvis and proximal extension of the thigh; it is zero in upright standing and has a range of −20 • < ϕ hip < 150 • , where ϕ hip > 0 • is defined to be flexion. The hip angle is considered in an analogous way for the GAS calculations.

Lower Extremity Dynamics
The dynamic model of the lower extremity follows the approach of Yamaguchi and Zajac [46]. In contrast to [46], in our model the ankle joint TC rotation axis is not perpendicular to the sagittal plane; moreover, other coordinate systems and angles are used. The individual determination of the centre-of-mass and centre-of-inertia is based on the relative data of de Leva [47]. The kinematics of the knee joint are based on [28], which assumes a fixed centre-of-rotation for the dynamics; this is in contrast to the knee kinematics, where the translation of the centre-of-rotation is still considered. The present model consists of three segments and three hinge joints ( Figure 4).
The derivation of the multibody dynamics follows the Euler-Lagrange equation approach. On the one hand, elastic and viscous joint properties are modelled as external forces. On the other, the potential energy consists of segmental mass m i in the corresponding centre-of-gravity, and the total kinetic energy results from the kinetic energy of the three modelled segments. Substitution of potential and kinetic energy into the Lagrangian functional and subsequent partial differentiation leads to the Euler-Lagrange equations M(q)q + D(q,q)q + k(q) = τ, where M(q) denotes the symmetric positive definite inertia matrix, q is the vector of generalised coordinates, and τ denotes the vector of external forces and torques acting on the lower extremity model. The matrix D(q,q) contains centrifugal and Coriolis terms and is defined as The last term in Equation (5) relates to potential energy and is defined as In order to determine the complete dynamics, the external joint torques are required, these external torques act on the dynamic system (5) through τ. However, for the hip joint no external torques are taken into account, because the hip joint is only considered for its influence on the length of the corresponding MTC that act on the knee joint. The corresponding external torque for the hip joint, τ 1 , is unknown. Consequently, the corresponding hip angle q 1 is considered as an external input signal. The equations of motion, given by Equation (5), have to be reduced with respect to the two unknowns q 2 and q 3 :M (q, q 1 )q +D(q,q, q 1 ,q 1 )q +k(q, q 1 ) =τ(q,q, q 1 ,q 1 ,q 1 ), which is an equivalent representation with respect to the two lower rows of Equation (5). The external joint torques τ j consist of actively generated joint torques τ T j , passively generated joint torques τ E j and τ V j as well as externally generated torques τ ext j like, for example, ground reaction forces. These torques are combined as The active torques for each joint are considered as the sum of all active contributions τ T i,j of the involved muscles i acting on the considered joint j where N j is the number of muscles that act on the joint j, r ma i,j is the lever arm of the corresponding muscle that acts on joint j and F T i is the tendon force that is calculated from the MTC dynamic models. These include an activation dynamics model and a modified extended Hill-model; these two dynamic models are not repeated here as they are extensively described in [28]. The inclusion of the passively generated joint torques is based on a simplification of the extended Hill model [28,48]. For that, the tendons are assumed to be stiff elements. Therefore, for each single muscle the parallel elastic and viscous element can be separated from the contractile element and combined into a single elastic and a single viscous element for each joint, that represents the combined viscoelastic torques around that joint. The viscous torques τ V j are described by [49] with damping constant K V j for the joint j. The modelling of elastic moment of the knee joint is based on [15] and is given as where the additional term τ * k accounts for the steepest increase of the elastic torque at full knee extension: The elastic torque at the ankle joint is modelled by the following form The nonlinear elastic torques for knee and ankle joint are shown in Figure 5. The dependency of the joint torques on the alignment of neighbouring joints is clearly shown and is due to the biarticular nature of the muscle groups involved. . Dynamic lower extremity model. System 0: global coordinate system; system 1: reference frame thigh; system 2: reference frame shank; system 3: reference frame foot (z 3 is the rotation axis of the TC).
Passive elastic joint torques for knee (left) and ankle (right) joint, based on parameters presented in [15].
Finally, ground reaction forces (GRF) have to be considered in case the lower extremity has, for example, contact with the ground. Assuming the point of action p GRF is known, the components of the GRF vector can easily be translated to model torques where o j denotes the joint centre and e z j is defined as the unit vector in the direction of the rotation axis of the joint j. The model is extended with muscle activation dynamics [28,50] and muscle contraction dynamics of an extended Hill type [28,51]. For the total of six muscles (RF, ST, GAS, VM, SOL, TA) that are included in the model, the state vectorx is used to describe activation and contraction dynamics. Denoting x 1 = q 2 = ϕ knee , x 2 = q 3 = ϕ ankle , x 3 =φ knee and x 4 =φ ankle , the dynamic model of the lower extremity can be given in the form of a state-space model aṡ In Equation (17), 16 denotes the state vector and u ∈ R 13 is the input vector, consisting of sEMG-signals, hip angle ϕ hip , components of GRF f GRF and three coordinates of point of action of the force p GRF . Thus, the full state vector is given by where a i is the activation of muscle i andF isom i is the isometric force of muscle i, normalised with the maximal isometric force of muscle i. The input vector is given by where s i is the rectified, normalised and filtered sEMG raw signal from muscle i, x CoP , y CoP , z CoP are the coordinates of the centre of pressure, and F GRF x , F GRF y , F GRF z are the components of the ground reaction force. Moreover,f (·) : R 12 × R 12 → R 12 denotes the nonlinear activation and contraction dynamics [28]. The output equation consists of joint angle measurements y = [x 1 , The state-space model given by Equation (17) can be reformulated in a compact form aṡ where f (·) : R 16 × R 13 → R 16 is a nonsmooth, nonlinear vector field and C ∈ R 2×16 is the output measurement matrix. Also considered are system noise n ∈ R 16 in Equation (20) and measurement noise m ∈ R 2 in Equation (21). The noise signals are assumed to be uncorrelated, mean free, white Gaussian processes with E[nn T ] = Q n and E[mm T ] = R m . Figure 6 is a block diagram of the lower extremity dynamic model. Note that the input sEMG-signals to the diagram are the bandpass and linear envelope filtered signals that are generated from the raw sEMG measurements (the process of sEMG signal filtering is described in Section 3.2). rma has an additional second GAS-component, since the GAS is spanning two joints Figure 6. Block diagram of the lower extremity dynamic model. The dashed boxes are unchanged compared to the model presented in [28], expcept for the additional muscles (SOL, TA). The bold boxes are fundamentally extended or new.

Methods and Material
This section presents the model-based filtering approach used to minimise the influence of modelling errors and measurement noise. In addition, the experimental setup and procedures used for the validation study are described. The filter of choice is cubature Kalman filter, since it increases numerical accuracy towards numerically problematic operations, such as the quadratic terms associated with the prediction of the covariance matrices [27].

Cubature Kalman Filter and Square-Cubature Kalman Filter
The cubature Kalman filter (CKF) was first presented in [27]. The CKF is different from, for example, the extended Kalman filter (EKF) or the unscented Kalman filter (UKF). A comparison between the CKF and UKF shows that the CKF is contained as a special case in the unscented transform [27]. However, in contrast to the UKF, the CKF does not use any negatively weighted sampling points and is, therefore, well suited for a numerically more stable implementation, the so-called square-root cubature Kalman filter (SCKF). The CKF is based on the nonlinear state-space model represented by a discretised version of Equations (20) and (21) x and on Bayes' equations. In Equations (22) and (23), k denotes sampling instances of discrete time, x k is the discrete state vector and y k is the discrete measurement vector. The associated noise covariances Q k and R k are the discretised versions of the noise covariances in Equations (20) and (21). f d (·) : R n × R m → R n and h d (·) : R n → R p are nonlinear nonsmooth vector fields. For the prediction step the predictive probability density results in where D k = {u i , y i } k i=1 denotes the input measurement pairs up to time k and p is the probability density. The posterior density is given by where c k denotes the normalising constant, given by The observation density p(y k+1 |x k+1 ) is determined by Equation (23). It should be noted that the multi-dimensional integrals of Equations (24) and (26) are not always solvable in complete form. However, for processes with a Gaussian probability distribution, Bayes' equations can be brought to a recursive form, where only the expected value and the covariance are to be determined. The solution of the multidimensional integral is reduced to the form of nonlinear function · Gaussian distribution. For the expected value and the covariances, the following equations are obtained where N (x;x, P) denotes a Gaussian distribution with mean valuex and variance P. The idea of the CKF is to approximate the integrals given in Equations (27)- (31) for the nonlinear vector fields with sets of weighted points. For these points ξ i and the weightings ω i , we ask for an arbitrary function g(·) : R n → R n and an argument x ∈ R n R n where I is the unity matrix of dimension n. In the CKF, the points are determined with a third-degree spherical-radial cubature rule. Then the points are given as where [1] i denotes the ith point that is obtained by permutation and sign change of the components of the n-dimensional generator (1, 0, . . . , 0) T for the set of points The integrals given in Equations (27)- (31) can be solved with the Equations (32)- (34) for Gaussian distributions. The resulting CKF algorithm is according to the UKF algorithm [52], with the exception of the CKF cubature sampling points X i which are chosen in a different way as Comparison of the sigma points of the UKF and the cubature points of the CKF shows that the CKF is a special case of the UKF [27]. However, in contrast to the unscented transformation, the CKF gives a mathematical justification for an accurate sampling of the weighted Gaussian integrals. To avoid numerical difficulties associated with limited numerical precision, the CKF is formulated in its SCKF version that avoids computation of the squareroot of the covariance matrices (positive definiteness of covariance matrices is ensured). In the following, for convenience, the SCKF is presented in a compact form.
The SCKF is initialised with S 0|0 that denotes the root factor of the error covariance P k|k = S k|k S T k|k . In the filter prediction phase, the cubature points are determined and propagated through the state equation and the predicted state is estimated The prediction of the error covariance and accordingly its root factor S k+1|k with P k+1|k = S k+1|k S T k+1|k follows a triangularisation(S = tria(A) denotes the triangularisation of the matrix A, with e.g., a QR-decomposition) with Q k = S Q,k S T Q,k and moreover The first steps of the filter correction phase can be formulated analogous to the filter prediction phase with R k = S R,k S T R,k and The cross-covariance is estimated explicitly (not as a root factor) where The Kalman gain matrix and the innovation are calculated as (The operator/denotes right-hand side matrix divisions. Thus, B/A is short for the solution of xA = B with a least squares approach).
where finally state and error covariance root factor update is given bŷ

Experimental Setup
The experimental setup consists of body-worn sensor systems and a custom built test-bench. Two synchronised sensor systems were used that were attached to the body of a test person. A surface electromyography (EMG) measurement was conducted with the MyoMuscle system (Noraxon Inc., Scottsdale, AZ, USA) following the SENIAM (Surface ElectroMyoGraphy for the Non-Invasive Assessment of Muscles) guidelines to place the electrodes [53]. Segmental orientation was measured with nine-degrees-of-freedom (9DOF) magnetic-angular-rate-gravity (MARG) sensors of the MyoMotion system (Noraxon Inc., Scottsdale, AZ, USA). Both systems consist of wireless sensor nodes attached to the test person, that allow the user to move freely during stiffness determination tests. The sensor systems allow for synchronised recordings with a standard computer (MyoResearch software, Noraxon Inc., Scottsdale, AZ, USA). sEMG signals are sampled at a sampling frequency of 1500 Hz (16 bit analog-digital conversion) and low-pass filtered (corner frequency 500 Hz), as well as high-pass filtered (corner frequency 10 Hz). The segmental orientation angles with respect to an inertially fixed reference system are provided with a sampling frequency of 100 Hz. The experimental setup is shown in Figure 7 and consists of a personal computer with data acquisition and real-time operating system (DS1104, dSpace GmbH, Paderborn, Germany) and test-bench integrated sensors. These sensors are mounted to the rotational axis of the test-bench and are used to determine the torque (DR-2477, Lorenz Messtechnik GmbH, Alfdorf, Germany) and the rotational angle (HDmag MHAD 50, Baumer Group, Fraunfeld, Switzerland). For the determination of a joint torque, a known axial reference inertia is employed in combination with a nonlinear dynamic model. The parameters of this nonlinear dynamic model were determined in an experimental identification procedure, described in [54]. Thus, by subtracting all of the test-bench torques from the torque sensor signal, the human joint torque that acts externally on the test-bench can be extracted and used as a reference for the filter validation. Finally, the neuromusculoskeletal model and the SCKF were implemented in Matlab (The MathWorks, Natick, MA, USA).

Simulation and Experimental Results
The model-based stiffness estimator was validated in simulations and in an experiment including five test persons. However, for the in silico tests the model was parameterised with values from literature [15,43,44,55]. Further parameters were taken from [28]. A table with model parameters and sources is given in the Appendix D, Tables A3 and A4.

Simulation Study
Simulation tests were conducted in static and dynamic situations. In the static case, the focus of investigation lies in the varying length and lever arms of single muscles that are investigated depending on the kinematics of the lower extremity model. Of interest are the biarticulate muscles, also with regard to changes in the neighbouring joints. In case of the ankle joint, this is the GAS. Figure 8 shows the lever arm of the GAS with respect to the ankle joint.
In addition to the lever arms that were calculated via the distance to the joint centre, lever arms that were calculated with the principle of virtual displacements are shown, where for the torque the lever arm of the muscle i at the joint j can be stated as The lever arm calculated from the principle of virtual displacements shows an offset when compared to the lever arms that were calculated from the model. On the one hand, these differences might be due to anatomical limitations in the attachment of the muscles and/or the reduction to the two-dimensional case. On the other hand, the fact that the lever arm is calculated via a rotation with a constant generalised coordinate axis might be the reason for the differences in the observed lever arm. However, the main shape of the two lever arms curves are similar. Other than the GAS, the ankle joint contains the uniarticular muscles SOL and TA for which the muscle and the lever arms are shown in Figure 9. As before, the model-calculated lever arms are compared to lever arms computed with the principle of virtual displacement. For the other muscles acting on the knee, analogous results are obtained. Moreover, the muscle length and lever arm values are in accordance with published values [56][57][58]. In the Appendix A, Figure A1 compares the muscle length to data from the literature. We conclude that a hinge model approach is sufficient for ankle joint modelling. After static simulation tests, the SCKF was verified in dynamic in silico tests. For that, the model of Equations (20) and (21) was discretised at a sampling time of T s with the resulting discrete state-space equations where n k and m k are the corresponding discretised versions of the system and measurement noise processes. In the simulation model, in comparison to the filter model, additive white Gaussian noise is added to the states and to the measurements. Moreover, the external hip velocity and acceleration are obtained from the hip position by means of differentiation (a low pass with relatively high crossover frequency is used to limit to the influence of high frequency noise). This filter is thereby implemented by using the following transfer function To test the stiffness estimation in an appropriate dynamic situation, all three angles (hip, knee and ankle) are initialised with an initial state of 0. Note that the hip flexion-extension angle is important for the full-order model for knee and ankle. The muscles that are considered in the lower extremity model are then activated in different patterns over a time window of 5 s, followed by a resting phase of 5 s. Figure 10 presents an example of the SCKF for the full-order nonlinear model. During the experiment, the hip angle is changed in a sinusoidal manner with an amplitude of 20 • and a frequency of π 2 rad/s. Filter internal states are used to calculate an overall ankle and knee torque, which is compared to the torque of the lower extremity model. As can be seen in Figure 10, the filter shows good state estimation performance, with slightly larger errors in the angular rates compared to the angles. Angular errors (RMS) remain below 0.4 • and angular rate errors (RMS) remain below 0.7 rad/s. Ankle and knee torques are estimated quite accurately, also during dynamic situations (e.g., during muscle activation).

Experimental Study
The substantial extension of the lower extremity dynamic model is the ankle model. To validate the ankle joint stiffness estimator under real-world experimental conditions, the test-bench described in Section 3.2 was used to provide reference measurements for ankle joint stiffness and angle. Before the measurement, Ag/AgCl surface electrodes were placed on the corresponding muscles and body-worn sensors of the MyoMuscle system were attached to a test person. In addition, the orientation of foot, shank and thigh was determined by applying 9DOF MARG sensor nodes of the MyoMotion system with straps to the respective segments. The lower extremity SCKF filter model had to be adapted to the experimental setup. Due to a connection of the ankle joint rotational axis with the rotational axis of the test-bench, there are no ground reaction forces. Furthermore, the torque sensor, used in the test-bench, measures additional external torques from the reference mass and viscous friction. These external torques were determined by employing a nonlinear dynamic model of the test-bench and subtracted from the sensor signal [28,54]. Finally, the dynamic model of the SCKF was reduced by assuming constant knee and hip angle values q 2 = ϕ knee and q 1 = ϕ hip , respectively. This simplification is assumed to be valid, since the test persons were instructed to sit upright and not to move their knee or their hip. Thus, in the modelq 1 =q 1 =q 2 =q 2 = 0 is assumed with subsequent simplifications in the model dynamics.
After model parameterisation and preparation of the experimental setup, a study with five male test persons was conducted (age [years]: 29 ± 5.7, height [cm]: 181.6 ± 1.82, weight [kg]: 79.2 ± 8.58). During the experiments, the test persons were instructed to conduct a number of different extensions/flexions after an initial maximum voluntary contraction test of muscles with relevance to the ankle joint. sEMG signals were then normalised to maximum voluntary contraction before they were applied to the filter. Characteristic experimental movements included plantar flexion and dorsiflexion, as well as coactivation of flexors and extensors. An overview of the experimental test protocol is given in the Appendix B (Table A1). An example of such an experiment with test person ID004 is shown in Figure 11; as can be seen, the SCKF is able to reconstruct the ankle torque based on surface EMG and segmental orientation measurements. Ankle joint torques of up to 12 Nm were generated during the measurements. Moreover, the joint torque estimation is accurate for muscle coactivation (phase (e) in Figure 11) as well as a dynamic cyclic movement (phase (f) in Figure 11). In addition to the model-based estimation of the ankle joint torque, the quasi-stiffness κ ankle was computed from the model states where r ma ankle = [r ma 1,ankle , r ma 2,ankle , . . . , r ma N,ankle ] T and F T ankle = [F T 1 , F T 2 , . . . , F T N ] T denote the lever arms and the forces of the muscles that were considered in the model, respectively. In the resulting stiffness in Figure 11, peaks on the ankle joint stiffness κ ankle can be observed that are of unphysiological value and dynamics. We assume that the noisy sEMG signals which are given to the contraction dynamics are responsible for the high dynamics in the observed stiffness values. If an average mean stiffness value is determined, stiffness values show good agreement with published values [59,60]. For a quantitative investigation of the estimation quality, a statistical analysis was conducted. The results ( Figure 12) show the deviation of the estimated filter torque from the reference torque of the test-bench at maximum plantar flexion. The average RMS error is 1.27 Nm. In addition to these results, in the Appendix C, Table A2 lists the estimation errors obtained for the different test persons.

Discussion and Conclusions
The goal of this study was to develop a novel model-based estimator for the determination of ankle joint stiffness, employing body-worn sensors only. For this, an existing sagittal plane dynamic model of the human knee was extended to describe a full sagittal plane model of the lower extremity. The dynamic model is thus driven by a measured hip angle (and corresponding time derivatives), where ankle and knee torques are estimated based on sEMG and segmental orientation data. To minimise the influence of measurement noise and model errors that would occur in a forward simulation, a nonlinear state-observer was developed. The filter of choice was the square-root version of the CKF because of its ability to handle nonsmooth, highly nonlinear dynamic systems, and no necessity to process Cholesky decompositions of the covariance matrix. The SCKF was able to reconstruct knee and ankle joint torque/quasi-stiffness in simulations, as well as ankle joint torque/quasi-stiffness in an experimental series with five test persons. However, in different ankle joint dynamic experiments a decreasing estimation performance was observed. Since this decrease in the performance was non-consistent with multiple measurements within the group of test persons, a general observer problem for measurements with higher dynamics can be excluded. An improved test bench is required to further investigate observed effects. In the current setup, the rotation axis of the test-bench is fixed, whereas the axis of the TC is moving. This might lead to involuntary inversion/eversion movements that could be responsible for the individually observed errors. Moreover, the estimation of the ankle joint stiffness shows unphysiologically high peaks, whereas the average mean of the estimated ankle joint stiffness lies within the physiological range. An obvious explanation for this effect is a noisy sEMG measurement, including movement artifacts. An additional explanation might the uncertainty in the muscle model time constants that lead to activation mismatches. To improve the performance of the model one option would be segmental parameterisation of a certain subset of model parameters. A combination of segmental parameterisation and calibration procedures, as described in [14,20], might improve the quality of our model. Moreover, a comparison of an extended model with elastic tendons to our model (concentrated visco-elastic elements) should be investigated. We are currently investigating these items, as well as application of the full-order observer for combined knee/ankle stiffness estimation.  Table A1. Overview of the experimental protocol.