## Abstract

Postural stability is important in everyday life as falls can cause severe injuries. Risk of injuries is higher in the elderly whose balance is often impaired. Modeling postural stability and the parameters that govern it is important to understand the balance mechanism and allow for the development of fall prevention strategies. Several mathematical models have been proposed to represent postural stability of bipeds. These models differ on the number of degrees-of-freedom (DOF) of the skeletal structure, force generation function for the muscle models, and capability to change their behavior as a function of the task. This work proposes a nonlinear model that captures fall recovery using a hip–ankle strategy. The muscle actuation is modeled as a third-order Poynting–Thomson's (PT) mechanical system where muscles and tendons are represented as lumped parameters actuating the aforementioned joints. Both a regression technique and a Kalman Filter (KF) are used to estimate the muscle–tendon parameters of the model. With a good model, the direct estimation of these parameters would allow clinicians to improve postural stability in the elderly, monitor the deterioration of the physical condition in individuals affected by neuro-degenerative diseases, and develop rehabilitation appropriate processes.

## 1 Introduction

Balance impairment affects an ever-growing percentage of the elderly population in the developed world. Injuries that are consequence of a fall include, but are not limited to, bruises, contusions, and bone fractures. When severe injuries occur, immobilization can limit functional autonomy for an extended period of time thus hampering the cardio-pulmonary system. This is often conducive to health complications or even death. Thus, fall prevention is paramount. It starts with the understanding of the mechanisms that humans use to control their balance and maintain a stable erect posture.

Balance is often modeled as a function of the number of degrees-of-freedom (DOF) of the system's kinematics, the order of the muscle models, and the capacity of the dynamic model to change its behavior in a task dependent manner, the latter applying to both the control parameters and the coefficients of the state equations. The number of kinematic DOFs depends on the number of joints that have a significant excursion from their equilibrium position. The actuation of the joints is often represented as a combination of elastic and viscous elements, representing the storage and dissipation of elastic energy within tendons and muscles. If the system is represented as an inverted pendulum, the number of DOFs is equivalent to the number of joints that the model has only if the actuators are first-order (i.e., a Kelvin–Voigt system with elastic and viscous elements in parallel). For higher order actuation systems, there can be DOFs that are internal to the actuators which manifest within the joints' kinematics as a superposition of multiple signals.

A dynamic system actuated by a second-order Kelvin–Voigt (KV) muscle–tendon model, where a single stiffness and damping elements are in parallel, is not adequate to represent tendons and contractile elements independently [1]. It was observed that a KV model is not able to properly describe the in vivo behavior of the muscle [2], and a third-order system was proposed in the form of a Poynting–Thomson (PT) model. In the PT model, an additional elastic element is positioned in series with a KV model to represent the tendon. While the PT model is an approximation limited to a linearized, but physiologically consistent, region of the force–length–velocity curve, it can provide a realistic account of the “contraction dynamics” which arises from the interaction between the stiffness of the tendon in series with the contractile muscle element. This model seems to provide the minimum increase in complexity that can account for the storage of elastic energy in the tendon. If the stiffness of the contractile element is sufficiently high, the tendon will allow for the storage of elastic energy that can be transformed in kinetic energy without being dissipated directly by the damping element (which is constrained in its movement by the elastic element in parallel with it).

Joint stiffness has often being proposed as a figure of merit to assess postural stability [3]. The recovery from fall to quiet standing has been modeled using the skeletal structure [4] actuated by first-order muscle models. The stiffness and damping parameters of these models can be estimated using a variety of techniques either in the frequency [5], time [68], and time–frequency domains [9,10]. On the other hand, recent literature has demonstrated that when modeling an actuator as second-order (and thus the whole musculoskeletal system as third-order), exponential stability might not be achievable within the physiological range of the parameters [11]. Using models with higher complexity could highlight different control strategies involving the synchronized co-ordination of multiple degrees-of-freedom instead of the mere increase in lower limb joint stiffness, which is more energy demanding [12].

The present work uses the hold-and-release paradigm [13], combined with a time-based system identification technique [14] to identify the parameters of a multiple DOF third-order skeletal-muscle model actuated at the hip and ankle. We modeled the nonlinear equations of motion of a double inverted pendulum and used both a linear regression and a Kalman filter (KF) to identify the stiffness associated with both the tendons and the muscle, considered as separate lumped parameters, as well as the damping of the muscles. This work is different from the previous literature [15] as it represents both the behavior at the hip and ankle, estimating a third-order nonlinear skeletal-muscle system with only one perturbation. Note that this study does not include the forces at the joints due to the joint ligaments that stabilize the joint. While this is a very important problem as described in Refs. [1618], we assumed that the ligaments maintain the joint stable during the physical activity. It can be demonstrated that while the ligaments provide a substantial force to prevent the relative translation between bones, their contribution to rotational stiffness around the joint is minimal. We will address this in the discussion.

In this work, we present a system identification technique that can identify the control parameters used for the recovery from a fall. The system that we modeled is nonlinear and high-order. This technique could be used as a diagnostic tool that can help clinicians to separately identify the properties of the tendons and muscles and thus help determine if an eventual impairment is due more to tissue atrophy or poor neural control.

## 2 Third-Order Model for a Balancing Human

Human movement along the sagittal plane can be modeled by means of a double inverted pendulum such as the one shown in Fig. 1(a). In the figure, the bottom segment is used to represent the legs while the top segment stands for the head, arms, and trunk. This model does not account for the flexing of the knee but is enough to study the hip control strategy as reported in Ref. [13]. The two segments are connected by muscles and tendons of different visco-elastic properties [11,19]. Additionally, the leg segment is considered to be attached to the ankle joint and rigidly connected to the ground. Finally, the legs and the feet are also connected using visco-elastic elements. In an attempt to determine the numerical value of these parameters while maintaining balance, the dynamics equations governing the model will be derived. They will then be used to write a nonlinear model suitable for estimation.

Based on Fig. 1(b) and applying the Euler–Lagrange methodology [20], we define mi as the mass of segment i, li as the distance measured between the two consecutive joints, ri gives the distance from the distal joint to the segment's center of mass, kmi and bmi are the joint's stiffness and damping parameters associated with the muscle element, and kti is the joint's stiffness associated with the tendon which is modeled as a purely elastic element. Because there are no bi-articular muscles spanning from below the ankle to above the hip, we are not including any cross-joint stiffness or damping.

The position of each mass in the system ($rmi$) is determined as a function of the generalized coordinates $q=[θ1φ1θ2φ2]T$ as
$rm1=−r1 sin θ1î+r1 cos θ1ĵ$
(1)

$rm2=(−l1 sin(θ1)+r2 sin(θ21))î+(l1 cos(θ1)+r2 cos(θ21))ĵ$
(2)

where θ21 = θ2 − θ1 and $[î,ĵ]$ are the basis for “a” reference frame as shown in Fig. 1(b).

The system's total kinetic energy (T) is obtained after deriving the velocity of each mass and is given by
$T=12(m1r12θ˙12+I1θ˙12+I2θ˙212 +m2(l12θ˙12+r22θ˙212−2l1r2θ˙1θ˙21 cos(θ2)))$
(3)

where Ii is the moment of inertia of the segment with respect to the center of gravity of link i and is a function of its geometry.

The system's potential energy (V) must include the effect of gravity and the energy stored by the elastic elements so that
$V=g(m1r1 cos(θ1)+m2(l1 cos(θ1)+r2 cos(θ21))) +12(kt1φ12+kt2φ22)+12(km1ψ12+km2ψ22)$
(4)

where g represents the acceleration due to gravity and ψi = θi − φi.

Finally, we consider the energy dissipated by the dampers [21] as
$D=12(bm1ψ˙12+bm2ψ˙22)$
(5)
The Lagrangian is defined as L = T − V. Considering the dynamic properties of the system [21], its motion is given by
$ddt∂L∂q˙−∂L∂qk+∂D∂q˙=∑j=1pFj·∂cj∂qk$
(6)

where Fj represents one of p external forces applied to the system at location cj.

Under no additional external forces ($∑j=1pFj·∂cj/∂q=0)$, the following relationships are found:
$0=(α+2β cos(θ2))θ¨1−(γ+β cos(θ2))θ¨2+ν1 +β sin(θ2)(θ˙22−2θ˙1θ˙2)−g(δ sin(θ1)−ε sin(θ21))$
(7)

$0=kt1φ1−ν1$
(8)

$0=γθ¨21+β(sin(θ2)θ˙12−cos(θ2)θ¨1)−gε sin(θ21)+ν2$
(9)

$0=kt2φ2−ν2$
(10)
where variables ν1 and ν2 have been defined as
$ν1=km1ψ1+bm1ψ˙1$
(11)

$ν2=km2ψ2+bm2ψ˙2$
(12)
and the model's inertial parameters are contained in the following expressions:
$α=I1+I2+m1r12+m2(l12+r22)β=m2l1r2γ=I2+m2r22δ=m1r1+m2l1ε=m2r2$
(13)
For the estimation of the visco-elastic parameters, it is convenient to write the previous relationships as functions of only θ1 and θ2 as they are directly measurable. This can be achieved by solving for ν1 and ν2 in Eqs. (8) and (10) and substituting in Eqs. (7) and (9). This will yield two expressions which can be used to solve for φ1 and φ2 in terms of θ1, θ2, and their first and second derivatives. Differentiating these expressions yields $φ˙1$ and $φ˙2$ as functions of up to the third derivative of θ. Finally, $φ1, φ2, φ˙1, φ˙2$ are substituted into Eqs. (11) and (12) so as to write them as functions of θ1, θ2 and their derivatives. More details on these operations are given in the Appendix. In this manner, the equations of motion can be written in the following form:
$Mq…+Cq¨+Sq˙+G=τ$
(14)

$q=[θ1θ2]T$
(15)
and
$M=[(α+2β cos(θ2))bm1−(γ+β cos(θ2))bm1−(γ+β cos(θ2))bm2γbm2]$
(16)

$C=[C11C12C21C22]$
(17)

$C11=(α+2β cos(θ2))k1−4β sin(θ2)θ˙2bm1$
(18)

$C12=−(γ+β cos(θ2))k1+β sin(θ2)(3θ˙21+θ˙1)bm1$
(19)

$C21=−(γ+β cos(θ2))k2+β(θ˙2+2θ˙1)sin(θ2)bm2$
(20)

$C22=γk2$
(21)

$S=[S11S12S21S22]$
(22)

$S11=−2β sin(θ2)θ˙2k1−2β cos(θ2)θ˙22bm1+g(δ cos(θ1)+ε cos(θ21))bm1$
(23)

$S12=β sin(θ2)θ˙2k1+(β cos(θ2)θ˙22+gε cos(θ21))bm1$
(24)

$S21=β sin(θ2)θ˙1k2+(β cos(θ2)θ˙1θ˙2+gε cos(θ21))bm2$
(25)

$S22=−gε cos(θ21)bm2$
(26)

$G=[−g(δ sin(θ1)−ε sin(θ21))k1gε sin(θ21)k2]$
(27)

$τ=−[θ1km1kt1+θ˙1bm1kt1θ2km2kt2+θ˙2bm2kt2]$
(28)

$k1=km1+kt1$
(29)

$k2=km2+kt2$
(30)

### 2.1 Parameter Estimation.

The dynamics of the third-order model (14) can be expressed as a linear function of the visco-elastic parameters such that
$Z=Hλ$
(31)

$=[z1z2]TH=[A11θ1θ˙1000000A24θ2θ˙2]$
(32)
where $λ=[λ1λ2…λ6]T$ is composed of combinations of these parameters such that
$kt1=λ3kt2=λ6bm1=kt12λ1kt1−λ2km1=bm1λ2kt1km2=kt2λ5kt2λ4−λ5bm2=kt2km2λ5$
and
$z1=γθ⃛2−αθ⃛1+β sin(θ2)(4θ¨1θ˙2−θ¨2(θ˙2+2θ˙21)) +β cos(θ2)(θ⃛21−θ⃛1+θ˙22(θ˙21−θ˙1)) +g(δ cos(θ1)θ˙1−ε cos(θ21)θ˙21)$
(33)

$z2=−γθ⃛21+β cos(θ2)(θ⃛1−θ˙12θ˙2) −β sin(θ2)θ¨1(θ˙2+2θ˙1)+gε cos(θ21)θ˙21$
(34)

$A11=−β cos(θ2)(θ¨21−θ¨1)+β sin(θ2)(θ˙22−2θ˙1θ˙2) +g(ε sin(θ21)−δ sin θ1)+αθ¨1−γθ¨2$
(35)

$A24=β(sin(θ2)θ˙12−cos θ2θ¨1)−gε sin(θ21)+γθ¨21$
(36)

#### 2.1.1 Estimation of Constant Parameters.

Consider the following linear system:
$Z=Hλ̂+ρ$
(37)

where Z is a measurement vector, $λ̂$ is a vector containing the estimated system's linear parameters, ρ contains the estimation error, and H is known as the configuration matrix. Each row of Z and H represents one equation which, when aggregated, can be used to solve the linear system defined by Eq. (31).

When the parameters in λ are constant, it is possible to find their value using a least squares approach. Ideally, there should be at least enough linearly independent measurements such that H is invertible. In practice, it is desirable to create an overdetermined system [22]. It is possible to find a vector $λ̂=H+Z$ such that it minimizes the Euclidean norm of ρ [23] by using the Moore–Penrose pseudo-inverse [24]
$H+=(HTH)−1HT$
(38)

#### 2.1.2 Recursive Methods.

When the parameter vector λ evolves through time, with known dynamics, a Kalman filter is more suitable for estimating its values than the pseudo-inverse [25]. The Kalman filter equations are written as
$λ̂k+=Aλ̂k−+Bμk$
(39)

$Pk+=APk−AT+Q$
(40)

$Kk=Pk−HkT(HkPk−HkT+R)−1$
(41)

$λ̂k=λ̂k−+Kk(Zk−Hkλ̂k−)$
(42)

$Pk=(I−KkHk)Pk−$
(43)

where A and B determine how the models parameters (λ) change over time; P, Q, and R are the estimation, process, and measurement variance and covariance matrices, respectively; Zk is the measurement vector; Hk is the configuration matrix; K is the Kalman gain; and I is an identity matrix of suitable size.

## 3 Visco-Elastic Parameter Estimation Using Kinematic Data

We estimated the apparent visco-elastic parameters of six human subjects (age = 27 ± 8, weight = 77.3 ± 8.7 kg, and height = 1.63 ± 0.14 m). All subjects were willing participants, had no prior history of balance issues, and were considered healthy.

The experiment followed the hold-and-release paradigm [13]. In which the subjects lean forward with the hip at a neutral position, with her/his arms crossed and held close to the chest, and both feet on the ground. After a random amount of time, the subject is suddenly released. This disturbance forces the subject to initiate a recovery from a potential fall. The experiment procedure is illustrated in Fig. 2.

To keep track of the subjects' oscillations, she/he was equipped with markers at the torso, ankle, and knee joints; from which joint angles about the sagittal plane could be computed. The experiment was recorded using an RGB camera with an average framerate of 22.5 Hz. The video was processed using the open source software kinovea which was developed for its applications in sport science [26,27]. The marker position was filtered using a low-pass Butterworth with a cutoff frequency of 10 Hz prior to obtaining joint angles. Angular velocity and acceleration were estimated using central differences so as to avoid phase shifts [22].

Angular position, velocity, acceleration, and jerk are used to update the configuration matrix (Hk) and measurement vector Zk at each time-step by substituting them into Eqs. (31)(36). These data were aggregated, and the visco-elastic parameters for each subject were obtained using a least-squares algorithm constrained to return positive parameter values. Note that this assumes that the parameters remain constant throughout the motion thus the actuators remain passive. Table 1 shows the estimated parameters for the six subjects as well as its average value and standard deviation rounded to the nearest integer.

## 4 Simulation Study

In order to test the parameter estimation procedure, a total of 24 distinct simulations were created to emulate the hold-and-release experiment. Each of these simulations represent a unique individual for whom mass, height, and visco-elastic parameters were generated randomly. For all subjects, their geometric parameters (m1, m2, l1, r1, and l2) were calculated from the anthropometric data presented by Winter [28]. Finally, an initial angular value for the segment representing the legs (θ1) was used as to mimic the disturbance applied to the subject while undergoing hold-and-release tests [13].

For each random subject, the third-order dynamic system presented in Eq. (14) was solved numerically using a fourth-order Runge–Kutta approximation with a fixed step of 10 ms. A total of 7 s of data were generated.

Once the simulated angular values were obtained, a random, normally distributed Gaussian noise was added to them. This noise was added in order to resemble the measurements obtained from marker based human motion tracking systems.

### 4.1 Simulations Based on Measured Visco-Elastic Parameters.

Twelve of the 24 simulated subjects presented visco-elastic parameters with the same average value and standard deviation observed in the human experiments described in Sec. 3. Relevant information can be found in Table 2.

### 4.2 Data Processing.

In order to create the configuration matrix (H) in Eq. (32) up to the third derivative of angular position must be obtained. Following the same procedure utilized with the experimental data, these derivatives were obtained numerically after a zero phase filter, low-pass Butterworth filter with a cutoff frequency of 10 Hz, was applied. An example of the noisy and filtered signals can be seen in Fig. 3.

### 4.3 Parameter Estimation.

After processing the simulated measurement data, the visco-elastic parameters for each of the 12 simulated subjects were estimated using (i) a least squares estimation by means of the configuration matrix pseudo-inverse (Pinv) and (ii) a KF. The KF equation parameters were chosen as follows. It was assumed that the visco-elastic parameters remained constant making A in Eq. (40) a suitable identity matrix and μk equal to zero. The model uncertainty Q was chosen as a small valued diagonal matrix where $Qi,i=1×10−3$. Measurement uncertainty was considered higher due to the noise introduced to angular measurements making $Ri,i=1×10−1$. The initial parameter estimate ($λ̂0$) was chosen randomly for each simulated subject. Finally, the initial estimate covariance matrix was chosen to reflect the uncertainty on the initial parameters with $Pi,i=1×103$.

Figure 4 shows the estimation results for all 12 simulated experiments. The markers on the six leftmost plots show the parameters obtained using the pseudo-inverse method while those on the six rightmost plots show those obtained using the KF. Each subplot has the values used for simulation on the abscissa and the estimated on the ordinate axis. The dashed line represents the ideal correspondence, that is, perfectly estimated parameters should fall on this line. For these estimated parameters, the coefficient of determination R2 against the data that originated them ranges from 0.57 to 0.99 (average 0.85 ± 0.16).

### 4.4 Simulated Subjects With Abnormal Visco-Elastic Parameters.

The remaining 12 simulated subjects presented larger visco-elastic parameters than those observed in the human experiments. Each parameter value used was obtained randomly from a normal distribution with an average and standard deviation as shown in Table 3. That is, 95% of the simulated values will be within two standard deviations of the mean.

For the simulated subjects, the estimation results, for both the Pinv and the KF, are also given in Table 3. Note that the Pinv and KF results are very similar, as they both give optimal estimations in the least-squares sense [25].

Figure 6 shows the estimated joint stiffness and damping for 12 different simulated subjects with abnormal visco-elastic parameters. The dashed line represents the ideal correspondence, that is, perfectly estimated parameters should fall on this line. This figure seems to indicate that the proposed parameter estimation technique is suitable for a wide range of combinations of the visco-elastic parameters. For the estimated parameters shown in Fig. 6, the coefficient of determination R2 against the data that originated them ranges from 0.38 to 0.97 (average 0.81 ± 0.19)

While the overall residual error (ρ) is small for the simulated experiments, we observe from Fig. 6 that not all parameters are correctly estimated. As an example, take the trial represented by the left facing triangle in this figure. Estimated values for the trunk segment (kt2, km2, and bm2) are close to the ones used for simulation. This is not so for those corresponding to the ankle joint (kt1, km1, and bm1). Note that the muscle presents a high stiffness (around 4000 N·m/rad) and low damping coefficients (around 500 Nm s/rad), indicating a large oscillation frequency. It is likely that this high oscillatory frequency of the ankle is lost due to the effect of the low pass filter, leading to incorrect parameter estimation.

## 5 Discussion

In this paper, we have analyzed a two-joint third-order coupled model of human standing actuated at each link. The coefficients of Eq. (14) are nonconstant and change nonlinearly as the function of the joint angles and their derivatives with respect to time.

The torque applied at each joint is generated by two Poynting–Thomson muscle–tendon systems where an equivalent tendon (represented by a linear spring) is in series to an equivalent muscle (represented by a parallel of a spring and a damper). Note that despite utilizing a linear model for the actuation, the equation of motion utilized for our model is highly nonlinear and has not been linearized due to the large variation of the angles in the act of recovering from a fall. We have proposed a method that allows for the estimation of the visco-elastic parameters actuating the joints. We have presented a least-square method that accurately estimates constant parameters and a Kalman filter that, with additional tuning, could allow for the estimation of the parameters when they are varying as a function of time.

The presentation of such a technique is important as it could help to quickly diagnose how the control parameters change with age and neurological diseases. Less complex models, where the movement of the hip is ignored, have highlighted that the stiffness of the ankle (modeled as a Kelvin–Voigt actuator) increases with age [6]. Furthermore, it is also important to analyze the effect that different impairments could have on the control parameters. Impairments are usually associated with a higher magnitude of the control parameters when compared to unimpaired individuals. This aspect has been reported widely in the literature for a variety of neurological impairments such as after a stroke [29], Parkinson's disease [30], and multiple scleroses [31] just to name a few. In order to take the large variation of the parameters which could occur due the aforementioned factors, we have performed the Monte Carlo analysis reported in Fig. 4 where the parameter to be identified was changed within a large range.

Increased stiffness was assumed as an attempt to increase the margin of stability so that oscillation around the ankle would have a smaller amplitude. Another way to easily control the dynamics of second-order systems is to increase the damping so to make the system exponentially stable. On the other hand, such a simple model requires a very high metabolic energy for the co-contraction necessary to increase the stiffness and could not explain behaviors where the stiffness of the ankle was lower than the one necessary to maintain postural stability, and yet the system was not falling [3]. Increasing the complexity of the system to a single-joint third-order model, provided mathematical proof that if specific ratios where maintained between the stiffness of the muscle and the tendon, the system could not be exponentially stable but could be stabilized asymptotically around a well prescribed orbit. This provides new insight into the mechanism by which this previously observed self-sustained oscillation could be obtained [12].

The estimation method proposed here will allow us to experimentally test the control strategy at both the hip and ankle and evaluate the interaction between the two joints. We speculate that by an appropriate modulation of the two joint movements, the two strategies illustrated above could be combined. For example, the ankle could work as a second-order system where the margin of stability is increased by increasing the stiffness while the hip strategy might be closer to a third-order system where a self-oscillation could provide prompt reactions while still allowing for a large displacement of the trunk. Indeed, we can observe in Table 1 that the muscle and tendon stiffness at the ankle are often very different from each other. Therefore, the behavior of the actuation system is dominated by the lower stiffness thus acting very closely to a second-order mass-spring system. On the other hand, the stiffness of the hip muscle actuator is comparable to the stiffness of the tendon and thus cannot be approximated as a second-order system.

The present work does not include the forces at the joints due to the joint ligaments that stabilize the joint. While this is a very important problem as described in Refs. [1618], we assumed that the ligaments maintain the joint stable during the physical activity. It can be demonstrated that while the ligaments provide a substantial force to prevent the relative translation between bones, their contribution to rotational stiffness around the joint is minimal.

As an example, we can assume that the cumulative translational stiffness of the muscle–tendon unit (MTU) is km = 100 N/m, and the stiffness of the ligament is kl = 1000 N/m. These are values that are assumed just for ease of the calculation but are of appropriate order of magnitude. For more precise values see Ref. [32]. On the other hand, the moment arm of the translational stiffness with respect to the point of rotation is much larger for the muscle. Comparing orders of magnitudes, we can loosely assume the moment arm of the MTU as rmt = 0.1 m and the moment arm of the ligament as rl = 0.001 m. When the rotational stiffness is calculated for both the muscle–tendon unit and the ligament using
$kθ=kmtrmt2+klrl2$
(44)

the square of the moment arm for the ligament will be many orders of magnitude smaller than the one for the MTU. Therefore, the resulting ligament rotational stiffness will be negligible compared to the rotational stiffness generated by the MTU. In the proposed example, the rotational stiffness for the ligaments bundle is about 1% of the MTU's rotational stiffness even though the translational stiffness of the ligament is much higher (Fig. 7).

Finally, while it could be argued that nonparametric techniques do not need to assume the order of the actuators' dynamics, they require the application of external random noise perturbations [3335]. Thus, there is a need for an external force/displacement generator that is difficult to use in the clinical setting and forces the system to work in un-natural conditions. The technique proposed in this work does not require an apparatus for perturbing the system. A perturbation with high power can be delivered by simply asking the subject to lean on the experimenter extended arm and then quickly releasing the support. This way a gravitational torque will act as a destabilizing perturbation on the subject that finds him/herself with a misaligned center of gravity with respect to the center of pressure. This can be quickly deployed in the clinical setting, together with a low cost motion tracking systems as those proposed, for example, in Refs. [6] and [7].

## Acknowledgment

The authors would like to thank the C. Christopher Cooney and Brian Jackman endowed Professorship.

## Funding Data

• Consejo Nacional de Ciencia y Tecnología (Catedra 2016-972; Graduate Scholarship Program; Funder ID: 10.13039/501100003141).

• Gannon University (C. Christopher Cooney and Brian Jackman Endowed Professorship; Funder ID: 10.13039/100006938).

• Universidad Autónoma de San Luis Potosí (C17-FAI-06-37.37; Funder ID: 10.13039/501100005324).

## Appendix

### Appendix: Derivation of the Third-Order Model

By solving Eq. (8), we find that ν1 = kt1φ1. After substituting this relationship into Eq. (7), we obtain that
$0=(α+2β cos(θ2))θ¨1−(γ+β cos(θ2))θ¨2+kt1 +β sin(θ2)(θ˙22−2θ˙1θ˙2)−g(δ sin(θ1)−ε sin(θ21))$
from where it is possible to solve for ϕ1
$φ1=1kt1(θ¨2(γ+β cos(θ2))−θ¨1(α+2β cos(θ2)) +g(δ sin(θ1)−ε sin(θ21))+β sin(θ2)(2θ˙1−θ˙2)θ˙2)$
Its time derivative ($ϕ˙1$) can be written as follows:
$φ˙1=1kt1(θ⃛2(γ+β cos(θ2))−θ⃛1(α+2β cos(θ2)) −g(ε cos(θ21)θ˙21+δ cos(θ1)θ˙1)−2β sin(θ2)(θ¨2+θ˙2)θ¨21 +β(cos(θ2)θ˙22+sin(θ2)(2θ˙1−θ˙2)θ¨2))$

Note that the third derivative of θ appears at this stage.

A similar procedure can be done to obtain ϕ2, since ν2 = kt2ϕ2. In this way, we find that
$φ2=1kt2(β(cos(θ2)θ¨1−sin(θ2)θ˙12)−γθ¨21+gε sin(θ21))φ˙2=−1kt2(β(sin(θ2)θ¨1(2θ˙1+θ˙2)+cos(θ2)(θ˙12θ˙2−θ⃛1)) +γθ⃛21−gε cos(θ21)θ˙21)$
Finally, knowing that ψ1 = θ1 − φ1 and ψ2 = θ2 − φ2, Eqs. (11) and (12) become
$kt1φ1=km1(θ1−φ1)+bm1(θ˙1−φ˙1)kt2φ2=km2(θ2−φ2)+bm2(θ˙2−φ˙2)$
which can then be written in terms of θ1, θ2, and their derivatives.

## References

References
1.
Tsuji
,
T.
,
Morasso
,
P. G.
,
Goto
,
K.
, and
Ito
,
K.
,
1995
, “
Human Hand Impedance Characteristics During Maintained Posture
,”
Biol. Cybern.
,
72
(
6
), pp.
475
485
.
2.
Kistemaker
,
D. A.
, and
Rozendaal
,
L. A.
,
2011
, “
In Vivo Dynamics of the Musculoskeletal System Cannot Be Adequately Described Using a Stiffness-Damping-Inertia Model
,”
PLoS One
,
6
(
5
), p.
e19568
.
3.
Morasso
,
P. G.
, and
Sanguineti
,
V.
,
2002
, “
Ankle Muscle Stiffness Alone Cannot Stabilize Balance During Quiet Standing
,”
J. Neurophysiol.
,
88
(
4
), pp.
2157
2162
.
4.
Chavez-Romero
,
R.
,
Cardenas
,
A.
,
Maya
,
M.
,
Vernaza
,
K. M.
, and
Piovesan
,
D.
,
2014
, “
Experimental Validation of Vision-Based System for the Characterization of Human Standing
,”
Latin American Congress of Automatic Control
(
IFAC-CLCA
), Cancún, México, Oct. 14–17, pp. 301–306.https://www.researchgate.net/publication/269690393_Experimental_Validation_of_Vision-Based_System_for_the_Characterization_of_Human_Standing
5.
de Gooijer-van de Groep
,
K. L.
,
de Vlugt
,
E.
,
de Groot
,
J. H.
,
van der Heijden-Maessen
,
H. C. M.
,
Wielheesen
,
D. H. M.
,
van Wijlen-Hempel
,
R. S.
,
Arendzen
,
J. H.
, and
Meskers
,
C. G. M.
,
2013
, “
Differentiation Between Non-Neural and Neural Contributors to Ankle Joint Stiffness in Cerebral Palsy
,”
J. Neuroeng. Rehabil.
,
10
(
1
), p.
81
.
6.
Chavez-Romero
,
R.
,
Cardenas
,
A.
,
Rendon-Mancha
,
J. M.
,
Vernaza
,
K. M.
, and
Piovesan
,
D.
,
2015
, “
Inexpensive Vision-Based System for the Direct Measurement of Ankle Stiffness During Quiet Standing
,”
ASME J. Med. Devices
,
9
(
4
), p.
041011
.
7.
Segura
,
M. E.
,
,
E.
,
Maya
,
M.
,
Cardenas
,
A.
, and
Piovesan
,
D.
,
2016
, “
Analysis of Recoverable Falls Via Microsoft Kinect: Identification of Third-Order Ankle Dynamics
,”
ASME J. Dyn. Syst. Meas. Control
,
138
(
9
), p.
091006
.
8.
Segura
,
M. E.
,
,
E.
,
Cardenas
,
A.
, and
Piovesan
,
D.
,
2015
, “
Time-Based Identification of Human Ankle Impedance Via Microsoft Kinect
,”
IEEE Signal Processing in Medicine and Biology Symposium
(
SPMB
), Philadelphia, PA, Dec. 12, pp.
1
5
.
9.
Piovesan
,
D.
,
Pierobon
,
A.
,
DiZio
,
P.
, and
Lackner
,
J. R.
,
2012
, “
Measuring Multi-Joint Stiffness During Single Movements: Numerical Validation of a Novel Time-Frequency Approach
,”
PloS One
,
7
(
3
), p.
e33086
.
10.
Piovesan
,
D.
,
Pierobon
,
A.
,
DiZio
,
P.
, and
Lackner
,
J. R.
,
2013
, “
Experimental Measure of Arm Stiffness During Single Reaching Movements With a Time-Frequency Analysis
,”
J. Neurophysiol.
,
110
(
10
), pp.
2484
2496
.
11.
Piovesan
,
D.
,
Pierobon
,
A.
, and
Mussa Ivaldi
,
F. A.
,
2013
, “
Critical Damping Conditions for Third Order Muscle Models: Implications for Force Control
,”
ASME J. Biomech. Eng.
,
135
(
10
), p.
101010
.
12.
Loram
,
I. D.
, and
Lakie
,
M.
,
2002
, “
Direct Measurement of Human Ankle Stiffness During Quiet Standing: The Intrinsic Mechanical Stiffness Is Insufficient for Stability
,”
J. Physiol.
,
545
(
3
), pp.
1041
1053
.
13.
Bortolami
,
S. B.
,
DiZio
,
P.
,
Rabin
,
E.
, and
Lackner
,
J.
,
2003
, “
Analysis of Human Postural Responses to Recoverable Falls
,”
Exp. Brain Res.
,
151
(
3
), pp.
387
404
.
14.
Roberts
,
T. J.
, and
Azizi
,
E.
,
2011
, “
Flexible Mechanisms: The Diverse Roles of Biological Springs in Vertebrate Movement
,”
J. Exp. Biol.
,
214
(
3
), pp.
353
361
.
15.
,
L. E.
,
Chavez-Romero
,
R.
,
Maya
,
M.
,
Cardenas
,
A.
, and
Piovesan
,
D.
,
2015
, “
Combining Genetic Algorithms and Extended Kalman Filter to Estimate Ankle's Muscle-Tendon Parameters
,”
ASME
Paper No. DSCC2015-9781.
16.
Caruntu
,
D. I.
, and
Hefzy
,
M. S.
,
2004
, “
3-D Anatomically Based Dynamic Modeling of the Human Knee to Include Tibio-Femoral and Patello-Femoral Joints
,”
ASME J. Biomech. Eng.
,
126
(
1
), pp.
44
53
.
17.
Caruntu
,
D. I.
,
Moreno
,
R.
, and
Freeman
,
R.
,
2017
, “
Knee Muscle and Ligament Forces During Drop Landing Exercise
,”
ASME
Paper No. IMECE2017-70982.
18.
Salinas
,
J. M.
, and
Caruntu
,
D. I.
,
2018
, “
2-D Inverse Dynamics Knee Model: Aligning Anatomical Knee Model With Knee Extension Kinematic Data Using Ligament Forces
,”
ASME
Paper No. DETC2018-85386.
19.
Piovesan
,
D.
,
Kennett
,
C.
,
Chavez-Romero
,
R.
,
Panza
,
M. J.
, and
Càrdenas
,
A.
,
2015
, “
Stiffness Boundary Conditions for Critical Damping in Balance Recovery
,”
ASME
Paper No. IMECE2015-50564.
20.
Meirovitch
,
L.
,
2010
,
Methods of Analytical Dynamics
, Wiley, Hoboken, NJ.
21.
Ogata
,
K.
,
1998
,
System Dynamics
, Vol.
3
,
Prentice Hall
,
Upper Saddle River, NJ
.
22.
Khalil
,
W.
, and
Dombre
,
E.
,
2004
,
Modeling, Identification & Control of Robots
,
Kogan Page Science
,
London
.
23.
Mooring
,
B. W.
,
Roth
,
Z. S.
, and
Driels
,
M. R.
,
1991
,
Fundamentals of Manipulator Calibration
,
Wiley
,
New York
.
24.
Penrose
,
R.
,
1955
, “
A Generalized Inverse for Matrices
,”
Math. Proc. Cambridge Philos. Soc.
,
51
(
3
), pp.
406
413
.
25.
Simon
,
D.
,
2006
,
Optimal State Estimation: Kalman, H Infinity, and Nonlinear Approaches
,
Wiley
, Hoboken, NJ.
26.
Charmant
,
J.
,
2018
, “Kinovea Version 0.8.15,” Bordeaux, France.
27.
El-Raheem
,
R. M. A.
,
Kamel
,
R. M.
, and
Ali
,
M. F.
,
2015
, “
Reliability of Using Kinovea Program in Measuring Dominant Wrist Joint Range of Motion
,”
Trends Appl. Sci. Res.
,
10
(
4
), pp. 224–230.
28.
Winter
,
D. A.
,
2009
,
Biomechanics and Motor Control of Human Movement
,
Wiley
, Hoboken, NJ.
29.
Lamontagne
,
A.
,
Malouin
,
F.
, and
Richards
,
C. L.
,
2000
, “
Contribution of Passive Stiffness to Ankle Plantarflexor Moment During Gait After Stroke
,”
Arch. Phys. Med. Rehabil.
,
81
(
3
), pp.
351
358
.
30.
Marusiak
,
J.
,
Kisiel-Sajewicz
,
K.
,
,
A.
, and
,
A.
,
2010
, “
Higher Muscle Passive Stiffness in Parkinson's Disease Patients Than in Controls Measured by Myotonometry
,”
Arch. Phys. Med. Rehabil.
,
91
(
5
), pp.
800
802
.
31.
Sinkjær
,
T.
,
Andersen
,
J. B.
,
Nielsen
,
J. F.
, and
Nielsen
,
J.
,
1996
, “
Impaired Stretch Reflex and Joint Torque Modulation During Spastic Gait in Multiple Sclerosis Patients
,”
J. Neurol.
,
243
(
8
), pp.
566
574
.
32.
Schmitz
,
A.
, and
Piovesan
,
D.
,
2016
, “
Development of an Open-Source, Discrete Element Knee Model
,”
IEEE Trans. Biomed. Eng.
,
63
(
10
), pp.
2056
2067
.
33.
Perreault
,
E. J.
,
Kirsch
,
R. F.
, and
Acosta
,
A. M.
,
1999
, “
Multiple-Input, Multiple-Output System Identification for Characterization of Limb Stiffness Dynamics
,”
Biol. Cybern.
,
80
(
5
), pp.
327
337
.
34.
Mirbagheri
,
M.
,
Barbeau
,
H.
, and
Kearney
,
R.
,
2000
, “
Intrinsic and Reflex Contributions to Human Ankle Stiffness: Variation With Activation Level and Position
,”
Exp. Brain Res.
,
135
(
4
), pp.
423
436
.
35.
Jalaleddini
,
K.
,
Golkar
,
M. A.
, and
Kearney
,
R. E.
,
2017
, “
Measurement of Dynamic Joint Stiffness From Multiple Short Data Segments
,”
IEEE Trans. Neural Syst. Rehabil. Eng.
,
25
(
7
), pp.
925
934
.