## 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 [6–8], 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. [16–18], 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 *m _{i}* as the mass of segment

*i*,

*l*as the distance measured between the two consecutive joints,

_{i}*r*gives the distance from the distal joint to the segment's center of mass,

_{i}*k*and

_{mi}*b*are the joint's stiffness and damping parameters associated with the muscle element, and

_{mi}*k*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.

_{ti}where *θ*_{21} = *θ*_{2} − *θ*_{1} and $[i\u0302,j\u0302]$ are the basis for “a” reference frame as shown in Fig. 1(b).

where *I _{i}* 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.

where *g* represents the acceleration due to gravity and *ψ _{i}* =

*θ*−

_{i}*φ*.

_{i}*L*=

*T*−

*V*. Considering the dynamic properties of the system [21], its motion is given by

where *F _{j}* represents one of

*p*external forces applied to the system at location

*c*.

_{j}*θ*

_{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 $\phi \u02d91$ and $\phi \u02d92$ as functions of up to the third derivative of

*θ*. Finally, $\phi 1,\u2009\phi 2,\u2009\phi \u02d91,\u2009\phi \u02d92$ 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:

### 2.1 Parameter Estimation.

#### 2.1.1 Estimation of Constant Parameters.

where **Z** is a measurement vector, $\lambda \u0302$ 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).

*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 $\lambda \u0302=H+Z$ such that it minimizes the Euclidean norm of

*[23] by using the Moore–Penrose pseudo-inverse [24]*

**ρ**#### 2.1.2 Recursive Methods.

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

**λ**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;

**Z**

*is the measurement vector;*

_{k}**H**

*is the configuration matrix;*

_{k}**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 (**H*** _{k}*) and measurement vector

**Z**

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

_{k}## 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 (*m*_{1}, *m*_{2}, *l*_{1}, *r*_{1}, and *l*_{2}) 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.

### 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\xd710\u22123$. Measurement uncertainty was considered higher due to the noise introduced to angular measurements making $Ri,i=1\xd710\u22121$. The initial parameter estimate ($\lambda \u03020$) 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\xd7103$.

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 *R*^{2} 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 *R*^{2} 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 (*k _{t}*

_{2},

*k*

_{m}_{2}, and

*b*

_{m}_{2}) are close to the ones used for simulation. This is not so for those corresponding to the ankle joint (

*k*

_{t}_{1},

*k*

_{m}_{1}, and

*b*

_{m}_{1}). 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. [16–18], 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.

*k*= 100 N/m, and the stiffness of the ligament is

_{m}*k*= 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

_{l}*r*= 0.1 m and the moment arm of the ligament as

_{mt}*r*= 0.001 m. When the rotational stiffness is calculated for both the muscle–tendon unit and the ligament using

_{l}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 [33–35]. 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

*ν*

_{1}=

*k*

_{t}_{1}

*φ*

_{1}. After substituting this relationship into Eq. (7), we obtain that

*ϕ*

_{1}

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

*ϕ*

_{2}, since

*ν*

_{2}=

*k*

_{t}_{2}

*ϕ*

_{2}. In this way, we find that

*ψ*

_{1}=

*θ*

_{1}−

*φ*

_{1}and

*ψ*

_{2}=

*θ*

_{2}−

*φ*

_{2}, Eqs. (11) and (12) become

*θ*

_{1},

*θ*

_{2}, and their derivatives.