# On Laterally Perturbed Human Stance: Experiment, Model, and Control

Authors: Dan Suissa; Michael Günther; Amir Shapiro; Itshak Melzer; Syn Schmitt
Journal: Applied Bionics and Biomechanics
 (2018)

---

## Abstract

Understanding human balance is a key issue in many research areas. One goal is to suggest analytical models for the human balance. Specifically, we are interested in the stability of a subject when a lateral perturbation is being applied. Therefore, we conducted an experiment, laterally perturbing five subjects on a mobile platform. We observed that the recorded motion is divided into two parts. The legs act together as a first, the head-arms-trunk segment as a second rigid body with pelvis, and the ankle as hinge joints. Hence, we suggest using a planar double-inverted pendulum model for the analysis. We try to reproduce the human reaction utilizing torque control, applied at the ankle and pelvis. The fitting was realized by least square and nonlinear unconstrained optimization on training sets. Our model is not only able to fit to the human reaction, but also to predict it on test sets. We were able to extract and review key features of balance, like torque coupling and delays as outcomes of the aforementioned optimization process. Furthermore, the delays are well within the ranges typically for such compensatory motions, composed of reflex and higher level motor control.


---

## Body

## 1. Introduction

How the human motor apparatus finds solutions to stabilize upright stance in the face of the body designed as a highly unstable multilinked-inverted pendulum has attracted the curiosity of modern scientists since at least 150 years [1], with quantitative analysis being intensified about 100 years ago by Murray et al. [2]. General neuromuscular control concepts have been derived therefrom by Asatryan and Feldman [3]; Feldman [4–6].Until today, analyses of human upright stance have been mainly performed in the sagittal plane. There has always been a separation in analyses of quiet and disturbed human stance. Contrary to quiet stance, disturbed stance has earlier been identified to be at least a double-inverted pendulum (DIP) problem by Nashner and McCollum [7]; Horak and Nashner [8]; Yang et al. [9]; Kuo [10]; Kuo et al. [11]; Runge et al. [12]; Kuo [13]; Fujisawa et al. [14]; Kiemel et al. [15]; Roth et al. [16] in the sagittal plane. Since about ten years ago, however, most researchers bore in mind a single-inverted pendulum (SIP) model by Geursen et al. [17] as an accompanying paradigm Winter et al. [18]; Morasso and Schieppati [19]; Morasso and Sanguineti [20]; Gage et al. [21] for their examinations of quiet stance in the sagittal plane. When examining quiet stance, even DIP models by Geursen et al. [17] equipped with just one additional degree of freedom (DOF) were not given attention until about a decade ago. Then, modern findings by Creath et al. [22]; Hsu et al. [23]; Saffer et al. [24]; Pinter et al. [25]; Günther et al. [26–28]; Yamamoto et al. [29] had finally pinpointed that quiet stance is in fact a multijoint phenomenon for which it takes at least triple-inverted pendulum (TIP) models to gain an understanding of its dynamics in the sagittal plane by Günther et al. [30]; Günther and Wagner [31]. To the best of our knowledge, only one study by Günther et al. [32] has examined rotations around the long axis of the body so far. And few studies by Bonnet et al. [33, 34] focused on the lateral (or frontal, resp.) plane, for which Day et al. [35] provided an analytic benchmark.It can be expected that stance dynamics in the lateral plane, be it externally disturbed or undisturbed (quiet), should show some similarity to those in the sagittal plane, with the primary difference that stance width is clearly higher than foot length, and the knees are rather compliant in sagittal but stiff in lateral direction. This leads us to the assumption that the DIP would be a well-suited model for lateral stance, accepting that the TIP is a well-suited “minimal model” for the sagittal plane, and that the knee joint is “missing” in the lateral direction. Anyone, gently pushed aside at the shoulder, or as in our perturbations case, undergoing a lateral shift of the ground, would probably respond with just a compensatory bending movement in the lower spine and hip, maybe combined with some load shift between both legs. The similarity to disturbed stance in the sagittal plane would be a DIP-like movement plus a shift in center of pressure (COP). Thus, aiming at an enhanced understanding of control strategies for lateral disturbances, we developed the aforementioned DIP model for simulating responses to such disturbances. We did this purposefully without modeling physiological structures for neural signal excitation and force generation (muscle-tendon units).Instead, we restricted ourselves to formulating abstract joint torque generators, as simple as possible. In particular, we based our control for both joints on an approach that can be seen as a “classical engineering method”: proportional-differential (PD) controllers. This abstract and simple approach allows a line of reasoning and inferences just reversal to an a priori implementation of physiologically based signal-processing and force-generating: initially neglecting all physiological properties (of control) but allowing some of them as a possible result of an optimization process can help to reveal their significance for the motor task under consideration. Accordingly, in this study, we chose the outlined approach to try and identify principal strategies and feedback loop properties supportive for stabilizing laterally disturbed human stance. The chosen approach particularly allows to disclose interdependences of control parameters like beneficial ratios (compare, e.g., Yang et al. [9]; Alexandrov et al. [36, 37]) and their relation to mechanical body design or beneficial time relations between sensor signals, actuator inputs, and actuator outputs. Based thereon, conclusions are possible relating to basic properties that may be imprinted into physiological and anatomical structures due to fundamental requirements for realizing the motor task (i.e., to the structures’ function): the mechanical body design (e.g., mechanical coupling through biarticular muscles) and the physical constraints (e.g., reaction times) the body is exposed to during the intended motor task.


## 2. Material and Methods

### 2.1. Inverse Kinematic Model

As can be seen in Figure1, the planar mechanical model consists of three rigid bodies. The feet are represented by a box ((0.2 × 0.2 × 0.02)[m], 1 kg). They are to be considered, a graphical component of the model, aside from lifting the ankle joint to a height of hf=0.02m. Furthermore, the DIP is composed of two rigid bodies: “Leg,” lumping both legs into one and “HAT,” depicting the head, arms, and trunk of the human body. The segments are defined between the centers (Ca, Ch, and Cs) by their masses, lengths, center of mass (COM) positions and inertia. Part of the anthropometric data were measured (i.e., total body height, segment lengths, and total body mass), while the segment masses and COM positions were calculated according to Winter [38]. Moments of inertia (2D) of the model segments with respect to their centers’ of mass are given in the following equation (infinitely thin, but rigid rod):
(1)JLeg=mLeg⋅lLeg2⋅112⋅I,JHAT=mHAT⋅lHAT2⋅112⋅I,where I is the identity matrix. The COMs are located at distance 0.45⋅lLeg (leg) and 0.5⋅lHAT (HAT) from the hip joint point, respectively (Figure 1).Figure 1
The marker positions (as in the experiment) at ankle, hip, and shoulder and their corresponding centers and the DIP model with: feet (heighthf), ankle joint (defined by Ca and angle: Θa), Leg segment between Ca and Ch (defined by: lLeg, COMLeg which is at dCOMHAT from Ch and by mLeg, and inertia JLeg), hip joint (Ch, Θh), and the HAT segment between Ch and Cs (lHAT, COMHAT at dCOMHAT from Ch, and by mHAT, JHAT).The kinematics of the model are as follows: the “HAT” segment rotates around the upper end of the “Leg” segment (Ch) and the “Leg” segment around the center of the feet (Ca). The feet are laterally translated from their initial position (=perturbation). Thus, the mechanical model’s three DOFs are the two angles Θa and Θh and the translational DOF x of the foot. Note that we did not constrained the DOFs which occur in the human body, for example, by ligaments. We assumed that Θa and Θh, as well as their first two derivatives and the acceleration of the perturbation p¨ (acting on Ca) may be sensed. These system states (and derivatives) are combined to form the input to the controllers that generate the torques applied to the ankle and hip. For further reading on the construction of the mechanical model utilizing the experimental data (i.e., marker positions), the physical properties of the model, and the implementation in SimMechanics, we refer to Appendix B.


### 2.2. Control

The DIP is perturbed by a motion signal (p,ṗ,p¨T) along x, deviating the feet laterally from their initial position. For our simulations, the perturbation signal p is experimentally obtained (lateral translation of Ca) and filtered by a PT2 low-pass filter at 6 Hz (as all experimentally obtained data). Note that the filter also provides the two derivatives ṗ,p¨ (for filtering, see Appendix A.4).As an abstraction of muscle force acting in humans ankle (A) and hip (H), joints are controlled by joint torques (MA,MH). Different control approaches are applied in order to represent central nervous system (CNS) strategies (note that I and IIa/b are reference strategies):
(I)
There is no control, simulating “no reaction” of the CNS to the perturbation (i.e., free falling/collapse of the body in response to the perturbation).
(II)
A “P(200)D(20)” controller with set point 0°, for both joints: (a) without delay and (b) with20 ms delay. This “technical” controller is representing the intuitive “control to 0°” strategy, similar to a SIP control strategy. Note that delays are introduced due to finite signal propagation velocities in the nervous system.
(III)
A dynamic spring-damper system is assumed by (2), based on that human posture is controlled by muscle-tendon complexes and sensor feedback, including proprioception and information about the perturbation:
(2)Mads←delay1M1=−k1⋅y1−d1⋅ẏ1+q1⋅p¨,Mhds←delay2M2=−k2⋅y2−d2⋅ẏ2+q2⋅p¨.The parametersX=delay1/2,k1/2,d1/2,q1/2 are acquired by unconstrained nonlinear optimization (see Appendix A.5), minimizing the squared error between given angular trajectories (i.e., experimental data) and those obtained in the simulation of the model (Eang):
(3)minXEang=minX∑iΘA,iexp−ΘA,isim2+ΘH,iexp−ΘH,isim2,where i is the sampling points. Following approaches for y1/2 for the controller III (2) are implemented to compare (a) uncoupled torques, (b) torques coupled by coupling angles one-to-one, and (c) torques coupled by coupling angles with the mixing ratios β:
(4a)IIIa:y1=Θa,y2=Θh;(4b)IIIb:y1=y2=Θa+Θh;(4c)IIIc:y1=β11⋅Θa+β12⋅Θh,y2=β21⋅Θa+β22⋅Θh.In approach (IIIc), the weight factors (i.e., mixing ratios)β equal the inertia tensor of the mechanical model (fitted to the current trials participant, not to be confused with the inertia of the model segments JLeg/HAT), numerically calculated by the inverse dynamics approach, which is explained in the next Section 2.2.1.


#### 2.2.1. Weights by Inverse Dynamics

Performing an inverse dynamics approach, the measured angular trajectories (Θaexpt,Θhexpt) were considered moving set points, of high gain (hg) PD controllers (P=2000, D=20, and N=100), to be retraced. The resulting torques Mahg and Mhhg are the bases for the computation of the weight factors.Alternatively starting from the Newtonian perspective, the joint torques may be simply approximated by a dependency of angular accelerations withM=I∗α where I is a moment of inertia and α an angular acceleration. With this ansatz
(5)MaLSQ=β11⋅Θ¨a+β12⋅Θ¨h+β31⋅p¨sim,MhLSQ=β21⋅Θ¨a+β22⋅Θ¨h+β32⋅p¨sim,we assume an inertia tensor β of dimension 3×2. The upper 2×2 part of β (given in kg⋅m2), corresponding to the contributions from the angular accelerations in the ankle and hip joints, approximates the 2D inertia tensor of the DIP (parametrized according to the participant). It is used as weighting factor between the ankle and hip angle in the states y1/2 of the controller (IIIc), see (4c). The lower part including β31 and β32 (given in kg⋅m) is not used further and thus ignored in the following.β is found from minimizing the squared and summed error of both Newtonian torque approximations Etot=Ea2+Eh2, over all time steps ts (LSQ):
(6)Ea=∑tsMahg−MaLSQ,Eh=∑tsMhhg−MhLSQ.In (4c) (IIIc), β11,12,21,22 (kg⋅m2) are multiplied with Θi (deg); thus, yi has the unit (kg⋅m2⋅deg). Since in (2) the left hand side has to be a torque (Nm=kg⋅m2/s2), the unit for ki is given by (1/deg⋅s2), for di is given by (1/deg⋅s), and for qi is given by (kg⋅m), since qi is multiplied with an acceleration.


### 2.3. Experiment

Experiments were conducted to (i) calculate the parameters of our assumed controllers, (ii) check the robustness of these controllers, and (iii) verify validity of the model. Five university students, four (1, 3–5) healthy and one (2) with occasionally light knee, hip, and back pain (i.e., rather “stiff” back and hip joint), of very similar physique participated in the experiment: all male, 26 ± 1.5 yrs, 69.4 ± 5.03 kg, 1.756 ± 0.03 m.The risk to fall for the healthy young adults is minimal, when perturbing stance with the BaMPer system; thus, no additional safety measures had to be taken. All participants gave informed consent. The experiments took place, all on the same date, in the “Schwartz Movement Analysis and Rehabilitation Laboratory,” Faculty of Health Sciences, Ben-Gurion University of the Negev, Beer-Sheva, Israel. The subjects took part on a purely voluntary basis. Studies on the BaMPer system were approved by the Helsinki Committee of Barzilai University Medical Center, Ashkelon, Israel (ClinicalTrials.gov registration numberNCT01439451).The perturbations were applied using the BalanceTutor© by MediTouch Ltd. which is based on the BaMPer system developed by Shapiro and Melzer [39]. Motion analysis was provided by the Ariel Performance Analysis Software (APAS version 10.100.0.1) running on a lab PC, utilizing two camcorders (Cannon NTC ZR100) with a frame rate of 59.94 fps. Six reflective LED markers (1 cm diameter) were placed on the participants’ front: one on each foot, placed directly above the ankle joint, two at the hip, on the most protruding points of each sides of pelvic bone, and another two at the shoulder, placed at the most protruding point of the clavicle. Additional markers were placed on the moving part of the treadmill in order to capture the perturbation, as well as a global reference marker on the lab wall.The participants were individually led into the lab and instructed to stand barefoot, upright, in hip-wide stance (on individually placed ground markers) on the BalanceTutor. They were asked to fold their hands behind their backs and keep them and the feet at their initial positions, throughout the whole experiment, if possible. Participants knew that they would receive lateral perturbations at their base of support, repeatedly during a period of a couple of minutes, executed through a shift of the laterally movable treadmill, to either the left or right.The set of at least 12 perturbations was of increasing magnitude (each to the left and right, for 10/14/18 cm distance with acceleration capped to either150cm/s2 or 200cm/s2) and was executed with no prior notice. Thus, perturbations were random in direction and time of onset, to mimic unexpected loss and recovery of balance. If arm movement or stepping was recorded, the trial was repeated once.


### 2.4. Average Controller and Prediction of Trials

From the set of valid experimental trials, we built a training (Str) and two test sets (Ste, Ste4). Training is done by optimizing the controller, individually for each trial and by then building the average over all obtained parameters. The final controller is then verified by predicting on the test sets. Due to strong qualitative similarities, we did not distinguish between the different perturbation magnitudes and directions, as well as the (in magnitude) different responses of the participants.


### 2.5. Overview: Data/Workflow

Figure2 depicts the data and workflow of the whole project.Figure 2
Data and workflow of the project. This sketch illustrates what was done, in which order and with the help of which tools. Large subscript A/H are used if the trajectory is meant, and small subscript a/h are used if only one value (e.g., sensor signals at a certain point in time).


## 2.1. Inverse Kinematic Model

As can be seen in Figure1, the planar mechanical model consists of three rigid bodies. The feet are represented by a box ((0.2 × 0.2 × 0.02)[m], 1 kg). They are to be considered, a graphical component of the model, aside from lifting the ankle joint to a height of hf=0.02m. Furthermore, the DIP is composed of two rigid bodies: “Leg,” lumping both legs into one and “HAT,” depicting the head, arms, and trunk of the human body. The segments are defined between the centers (Ca, Ch, and Cs) by their masses, lengths, center of mass (COM) positions and inertia. Part of the anthropometric data were measured (i.e., total body height, segment lengths, and total body mass), while the segment masses and COM positions were calculated according to Winter [38]. Moments of inertia (2D) of the model segments with respect to their centers’ of mass are given in the following equation (infinitely thin, but rigid rod):
(1)JLeg=mLeg⋅lLeg2⋅112⋅I,JHAT=mHAT⋅lHAT2⋅112⋅I,where I is the identity matrix. The COMs are located at distance 0.45⋅lLeg (leg) and 0.5⋅lHAT (HAT) from the hip joint point, respectively (Figure 1).Figure 1
The marker positions (as in the experiment) at ankle, hip, and shoulder and their corresponding centers and the DIP model with: feet (heighthf), ankle joint (defined by Ca and angle: Θa), Leg segment between Ca and Ch (defined by: lLeg, COMLeg which is at dCOMHAT from Ch and by mLeg, and inertia JLeg), hip joint (Ch, Θh), and the HAT segment between Ch and Cs (lHAT, COMHAT at dCOMHAT from Ch, and by mHAT, JHAT).The kinematics of the model are as follows: the “HAT” segment rotates around the upper end of the “Leg” segment (Ch) and the “Leg” segment around the center of the feet (Ca). The feet are laterally translated from their initial position (=perturbation). Thus, the mechanical model’s three DOFs are the two angles Θa and Θh and the translational DOF x of the foot. Note that we did not constrained the DOFs which occur in the human body, for example, by ligaments. We assumed that Θa and Θh, as well as their first two derivatives and the acceleration of the perturbation p¨ (acting on Ca) may be sensed. These system states (and derivatives) are combined to form the input to the controllers that generate the torques applied to the ankle and hip. For further reading on the construction of the mechanical model utilizing the experimental data (i.e., marker positions), the physical properties of the model, and the implementation in SimMechanics, we refer to Appendix B.


## 2.2. Control

The DIP is perturbed by a motion signal (p,ṗ,p¨T) along x, deviating the feet laterally from their initial position. For our simulations, the perturbation signal p is experimentally obtained (lateral translation of Ca) and filtered by a PT2 low-pass filter at 6 Hz (as all experimentally obtained data). Note that the filter also provides the two derivatives ṗ,p¨ (for filtering, see Appendix A.4).As an abstraction of muscle force acting in humans ankle (A) and hip (H), joints are controlled by joint torques (MA,MH). Different control approaches are applied in order to represent central nervous system (CNS) strategies (note that I and IIa/b are reference strategies):
(I)
There is no control, simulating “no reaction” of the CNS to the perturbation (i.e., free falling/collapse of the body in response to the perturbation).
(II)
A “P(200)D(20)” controller with set point 0°, for both joints: (a) without delay and (b) with20 ms delay. This “technical” controller is representing the intuitive “control to 0°” strategy, similar to a SIP control strategy. Note that delays are introduced due to finite signal propagation velocities in the nervous system.
(III)
A dynamic spring-damper system is assumed by (2), based on that human posture is controlled by muscle-tendon complexes and sensor feedback, including proprioception and information about the perturbation:
(2)Mads←delay1M1=−k1⋅y1−d1⋅ẏ1+q1⋅p¨,Mhds←delay2M2=−k2⋅y2−d2⋅ẏ2+q2⋅p¨.The parametersX=delay1/2,k1/2,d1/2,q1/2 are acquired by unconstrained nonlinear optimization (see Appendix A.5), minimizing the squared error between given angular trajectories (i.e., experimental data) and those obtained in the simulation of the model (Eang):
(3)minXEang=minX∑iΘA,iexp−ΘA,isim2+ΘH,iexp−ΘH,isim2,where i is the sampling points. Following approaches for y1/2 for the controller III (2) are implemented to compare (a) uncoupled torques, (b) torques coupled by coupling angles one-to-one, and (c) torques coupled by coupling angles with the mixing ratios β:
(4a)IIIa:y1=Θa,y2=Θh;(4b)IIIb:y1=y2=Θa+Θh;(4c)IIIc:y1=β11⋅Θa+β12⋅Θh,y2=β21⋅Θa+β22⋅Θh.In approach (IIIc), the weight factors (i.e., mixing ratios)β equal the inertia tensor of the mechanical model (fitted to the current trials participant, not to be confused with the inertia of the model segments JLeg/HAT), numerically calculated by the inverse dynamics approach, which is explained in the next Section 2.2.1.


### 2.2.1. Weights by Inverse Dynamics

Performing an inverse dynamics approach, the measured angular trajectories (Θaexpt,Θhexpt) were considered moving set points, of high gain (hg) PD controllers (P=2000, D=20, and N=100), to be retraced. The resulting torques Mahg and Mhhg are the bases for the computation of the weight factors.Alternatively starting from the Newtonian perspective, the joint torques may be simply approximated by a dependency of angular accelerations withM=I∗α where I is a moment of inertia and α an angular acceleration. With this ansatz
(5)MaLSQ=β11⋅Θ¨a+β12⋅Θ¨h+β31⋅p¨sim,MhLSQ=β21⋅Θ¨a+β22⋅Θ¨h+β32⋅p¨sim,we assume an inertia tensor β of dimension 3×2. The upper 2×2 part of β (given in kg⋅m2), corresponding to the contributions from the angular accelerations in the ankle and hip joints, approximates the 2D inertia tensor of the DIP (parametrized according to the participant). It is used as weighting factor between the ankle and hip angle in the states y1/2 of the controller (IIIc), see (4c). The lower part including β31 and β32 (given in kg⋅m) is not used further and thus ignored in the following.β is found from minimizing the squared and summed error of both Newtonian torque approximations Etot=Ea2+Eh2, over all time steps ts (LSQ):
(6)Ea=∑tsMahg−MaLSQ,Eh=∑tsMhhg−MhLSQ.In (4c) (IIIc), β11,12,21,22 (kg⋅m2) are multiplied with Θi (deg); thus, yi has the unit (kg⋅m2⋅deg). Since in (2) the left hand side has to be a torque (Nm=kg⋅m2/s2), the unit for ki is given by (1/deg⋅s2), for di is given by (1/deg⋅s), and for qi is given by (kg⋅m), since qi is multiplied with an acceleration.


## 2.2.1. Weights by Inverse Dynamics

Performing an inverse dynamics approach, the measured angular trajectories (Θaexpt,Θhexpt) were considered moving set points, of high gain (hg) PD controllers (P=2000, D=20, and N=100), to be retraced. The resulting torques Mahg and Mhhg are the bases for the computation of the weight factors.Alternatively starting from the Newtonian perspective, the joint torques may be simply approximated by a dependency of angular accelerations withM=I∗α where I is a moment of inertia and α an angular acceleration. With this ansatz
(5)MaLSQ=β11⋅Θ¨a+β12⋅Θ¨h+β31⋅p¨sim,MhLSQ=β21⋅Θ¨a+β22⋅Θ¨h+β32⋅p¨sim,we assume an inertia tensor β of dimension 3×2. The upper 2×2 part of β (given in kg⋅m2), corresponding to the contributions from the angular accelerations in the ankle and hip joints, approximates the 2D inertia tensor of the DIP (parametrized according to the participant). It is used as weighting factor between the ankle and hip angle in the states y1/2 of the controller (IIIc), see (4c). The lower part including β31 and β32 (given in kg⋅m) is not used further and thus ignored in the following.β is found from minimizing the squared and summed error of both Newtonian torque approximations Etot=Ea2+Eh2, over all time steps ts (LSQ):
(6)Ea=∑tsMahg−MaLSQ,Eh=∑tsMhhg−MhLSQ.In (4c) (IIIc), β11,12,21,22 (kg⋅m2) are multiplied with Θi (deg); thus, yi has the unit (kg⋅m2⋅deg). Since in (2) the left hand side has to be a torque (Nm=kg⋅m2/s2), the unit for ki is given by (1/deg⋅s2), for di is given by (1/deg⋅s), and for qi is given by (kg⋅m), since qi is multiplied with an acceleration.


## 2.3. Experiment

Experiments were conducted to (i) calculate the parameters of our assumed controllers, (ii) check the robustness of these controllers, and (iii) verify validity of the model. Five university students, four (1, 3–5) healthy and one (2) with occasionally light knee, hip, and back pain (i.e., rather “stiff” back and hip joint), of very similar physique participated in the experiment: all male, 26 ± 1.5 yrs, 69.4 ± 5.03 kg, 1.756 ± 0.03 m.The risk to fall for the healthy young adults is minimal, when perturbing stance with the BaMPer system; thus, no additional safety measures had to be taken. All participants gave informed consent. The experiments took place, all on the same date, in the “Schwartz Movement Analysis and Rehabilitation Laboratory,” Faculty of Health Sciences, Ben-Gurion University of the Negev, Beer-Sheva, Israel. The subjects took part on a purely voluntary basis. Studies on the BaMPer system were approved by the Helsinki Committee of Barzilai University Medical Center, Ashkelon, Israel (ClinicalTrials.gov registration numberNCT01439451).The perturbations were applied using the BalanceTutor© by MediTouch Ltd. which is based on the BaMPer system developed by Shapiro and Melzer [39]. Motion analysis was provided by the Ariel Performance Analysis Software (APAS version 10.100.0.1) running on a lab PC, utilizing two camcorders (Cannon NTC ZR100) with a frame rate of 59.94 fps. Six reflective LED markers (1 cm diameter) were placed on the participants’ front: one on each foot, placed directly above the ankle joint, two at the hip, on the most protruding points of each sides of pelvic bone, and another two at the shoulder, placed at the most protruding point of the clavicle. Additional markers were placed on the moving part of the treadmill in order to capture the perturbation, as well as a global reference marker on the lab wall.The participants were individually led into the lab and instructed to stand barefoot, upright, in hip-wide stance (on individually placed ground markers) on the BalanceTutor. They were asked to fold their hands behind their backs and keep them and the feet at their initial positions, throughout the whole experiment, if possible. Participants knew that they would receive lateral perturbations at their base of support, repeatedly during a period of a couple of minutes, executed through a shift of the laterally movable treadmill, to either the left or right.The set of at least 12 perturbations was of increasing magnitude (each to the left and right, for 10/14/18 cm distance with acceleration capped to either150cm/s2 or 200cm/s2) and was executed with no prior notice. Thus, perturbations were random in direction and time of onset, to mimic unexpected loss and recovery of balance. If arm movement or stepping was recorded, the trial was repeated once.


## 2.4. Average Controller and Prediction of Trials

From the set of valid experimental trials, we built a training (Str) and two test sets (Ste, Ste4). Training is done by optimizing the controller, individually for each trial and by then building the average over all obtained parameters. The final controller is then verified by predicting on the test sets. Due to strong qualitative similarities, we did not distinguish between the different perturbation magnitudes and directions, as well as the (in magnitude) different responses of the participants.


## 2.5. Overview: Data/Workflow

Figure2 depicts the data and workflow of the whole project.Figure 2
Data and workflow of the project. This sketch illustrates what was done, in which order and with the help of which tools. Large subscript A/H are used if the trajectory is meant, and small subscript a/h are used if only one value (e.g., sensor signals at a certain point in time).


## 3. Results

### 3.1. Experiment

#### 3.1.1. Perturbation

A set of 12 perturbations was executed in the experiment. They differ in magnitude and direction but not qualitatively. In Figure3, the motion signal used for perturbing the DIP model is shown for a test case. The test case is a trial by participant (5) (which is 14cm displacement and 200cm/sec2 acceleration trial; from now on referred to as (5)14200lt), as this is closest to the average in the experiment. Trial (5)14200lt is shown in Figure 4 as an example. The whole motion, from preperturbation (A), to perturbation (B, C), to initial participant response (D), to full balance recovery (E, F), is explained in the following.Figure 3
Filtered perturbation signal. For each trial, a qualitatively similar signal, containing change in position, velocity, and acceleration, was used. Quantitatively the perturbation signals differed in direction, thus to the left or right and in displacement and velocity.Figure 4
(a) Typical frames of the experiment, from a single trial, at the distinguished phases (A, B, C, D, E, F). (b) Corresponding trajectories of the hip and ankle angle, horizontal lines correspond to the phases. See also AppendixD for further reading.
(a)
(b)


#### 3.1.2. Results of the Experiment and Strategy-Dependent Variations

A simple segment length criterion (see AppendixA.2) identified two participants (1, 2), which deviated from the basic strategy, invalidating almost all their trials due to length changes of up to 15cm in the Leg segment and up to 6cm in the HAT segment. This was confirmed visually by the video, where in nearly all trials, participant (1) is extremely bending the knees and rotating the hip and spine horizontally, resulting in heavy A/P tilting. Participant (2) on the other hand acts extreme stiff around the lumbar spine and hip, resulting in a lift of one or both feet, see also Section 4.2.Loss of balance (i.e., falling) indeed never occurred, as participants used arm raising or stepping if needed, resulting in71/87 (81.6%) initially marked as “successful trials” (no arm movement and intentional stepping detected). Aside from that, within the participants (3, 4, 5) exhibiting a “proper” execution of the strategy (e.g., no A/P movement), only 5/36 trials (13.9%) had to be excluded.Participant (4), as can also be seen in Figure5 exhibited bending of the lumbar spine (just like participant (1) and in contrast to (2, 3, 5)). However in his case, this led to changes of the segment length within the error margin r, and thus, his trials were not excluded. Here, the hip angle excursions were up to 3.5 times higher than those of participants (3 and 5).Figure 5
Unfiltered angular trajectories for valid trials of the experiment. In the (a) plot are the perturbations to the left, and in the (b) are perturbations to the right. The solid line (blue) is the hip the dashed line (red) is the ankle angle trajectory. All the “extreme” hip angle peaks stem from participant (4).
(a)
(b)Within the valid trials, the ankle bends to a peak 5–7.5° while the “extended hip” (i.e., hip + lower spine) bends to a peak of up to 10° (34° for participant (4)). The higher variety within the hip stems from bending of the lumbar spine, which was accepted if segment lengths did not change significantly and thus lumped into the “hip” DOF. However, bending of the spine in the coronal plane, together with bending of the hip in the sagittal plane (often together with a rotation in the transversal plane), occurs in some trials, in all participants and thus, leading to greater segment length changes. Such trials were excluded from our data set as our planar model could not cope with them. In Figure5, valid trials of the experiment are plotted, while in Table 1, the mean of all max angles and angular velocities over all valid trials is given for each participant.Table 1
Mean and standard deviation of the maximal absolute angle and angular velocity for the ankle and hip for all valid trials.
Participant
Θamax
Θhmax
Θ̇amax
Θ̇hmax
(3)
5.26±1.94deg
10.38±5.42deg
19.88±7.09deg/s
60.29±26.98deg/s
(4)
3.73±1.55deg
22.40±8.76deg
24.00±8.73deg/s
106.7±39.45deg/s
(5)
5.62±0.9deg
7.45±2.57deg
25.94±4.08deg/s
50.14±10.42deg/s


### 3.2. Model and Control

#### 3.2.1. Reference Controller (I) and (IIa/b)

In Figure6, we are presenting a comparison of results obtained by using “no control” (I), “P(200)D(20)” (IIa), and “P(200)D(20) + 20 ms delay” (IIb), to the results obtained in the experiment. While, obviously, none of the approaches can reproduce the experimental data, it is still possible to stabilize the system using a common PD controller (which is making the model behave very similar to a SIP), however, only if there is no delay introduced in the sensor signals or the controller output. Note that delay on either the output or delay on the sensor signals leads to a similar result (Figure 6(b)) and that delays are inevitable in humans. Consequently, any simple PD control must fail to stabilize in the biologically restricted system.Figure 6
(a) Experiment trajectory in blue, “no control” (I) in red, and PD control (IIa) in dashed red. (b) Experiment trajectory in blue and PD control with 20 ms delay (IIb) in dashed red. First plot is the ankle trajectory, second is the hip, and third is the perturbation. While with no control applied, the system obviously collapses, also simple PD control, while being able to stabilize was not able to match the angular trajectories of the experiment. If a delay is introduced, the system becomes highly unstable. The PD controller used here is the MATLAB : Simulink build in PID controller withP=200, D=20, I=0, and N=100.
(a)
(b)


#### 3.2.2. Controller (IIIa/b)

Figure7 depicts the best fits of the first two approaches (IIIa, b) in comparison to the test case (5)14200lt. The parameters of the dynamic system (IIIa/b) are optimized such that Θa/hsimt are fitted to Θa/hexpt of a single experimental trial ((5)14200lt), which is done by unconstrained nonlinear optimization.Figure 7
(a) Controller IIIa, best fit, with an remaining error of 45.73deg/ts. (b) Controller IIIb, best fit, with an remaining error of 2.7 deg/ts. Both fitted to the same experimental trial (5)14200lt. First plot is the ankle trajectory, second is the hip, and third is the perturbation. While IIIa fails to stabilize (no torque coupling), even the introduction of simple one to one torque coupling (IIIb) is enough to (1) stabilize the system and (2) almost reproduce the experimental trajectories.
(a)
(b)The approach (IIIa) was invalidated, as it completely fails to reproduce the angular trajectories (not stabilizing the mechanical model at all).The approach (IIIb), which, by using a linear combination of both angles, is a lot closer to the sought solution, is nevertheless also invalid: the remaining error after optimization is big, and the controller is unable to stabilize the mechanical model.


#### 3.2.3. Controller (IIIc)

First, the weights fory1/2 (see (4c)) are calculated by an inverse dynamics approach via least square optimization. They are specific to the participant (not to the perturbation), that is, β5 is the same for all trials of participant (5), within small numerical differences. Thus, the weights can be calculated a priori, individually for each participant. The weights of participant (5) are β11=1.12,β12=β21=0.23,β22=0.06 (kg⋅m2). Note that it is a result of the simulation that β12=β21, within an accuracy 0.1%. The before mentioned unconstrained nonlinear optimization algorithm is used to find well-suited control parameters X=k1,d1,q1,delay1,k2,d2,q2,delay2.“Optimal” (i.e., local minimum) values for ((5)14200lt) are(7)k1=0.721deg⋅s2,d1=2.041deg⋅s,q1=0.83kg⋅m,delay1=0.06s,k2=0.661deg⋅s2,d2=1.71deg⋅s,q2=0.57kg⋅m,delay2=0.13s.The left plot in Figure8 depicts the “exact” reproduction of Mahg and Mhhg by the LSQ algorithm (used to find β, remember (5)) and the trajectories of the corresponding controller, utilizing the upper 2×2 part of β, in comparison to the experiments on the right.Figure 8
(a) The torques generated by the high gain PD:Mihg versus MiLSQ, i∈a,h with virtually MLSQ=Msim. y-axis is torque in Nm over time steps of the simulation on the x-axis. (b) Response of controller IIIc. The best fit, with an remaining error of 0.16 deg/ts. First plot is the ankle trajectory, second is the hip, and third is the perturbation. Finally, the controller IIIc with torque coupling according to the numerical inertia tensor β realized the almost exact reproduction of the experiments.
(a)
(b)


### 3.3. Average Controller “IIIcavrg”

Of the 31 valid trials, we took 10 for the training setStr and 11 for the test set Ste. Participant (4) is not part of training and test set due to his great (quantitative) deviation in strategy (extreme hip angles by bending the lower spine a lot more than other participants); he was treated separately by checking the performance, by the set Str parametrized average controller, on all his trials.In Table2, the average of all parameters after optimization on the training set is listed. The column “Error” is the average error made by the controller optimized for each specific trial in Str. This value is useful comparing the average controller error on the training and test set to the error made by the controller optimized to specific trials. The mean error thus is 0.4 deg/ts, while the best fit is trial (5)14200lt with 0.18 deg/ts and the worst fit is trial (3)14150lt with 0.58 deg/ts.Table 2
Average parameter setX for controller IIIc on the training set Str referred to by IIIcavrg.
Controller
Error
k1
d1
q1
delay1
k2
d2
q2
delay2
IIIc
0.35deg/s
1.181/deg⋅s2
1.971/deg⋅s2
0.45kg⋅m
0.06 s
0.851/deg⋅s2
1.731/deg⋅s
0.53kg⋅m
0.17 s
IIIc.1
0.54deg/ts
1.271/deg⋅s2
2.531/deg⋅s
—
0.06 s
1.191/deg⋅s2
2.061/deg⋅s
—
0.12 s
IIIc.2
0.55deg/ts
—
2.531/deg⋅s
—
0.06 s
—
2.051/deg⋅s
—
0.12 sNote that the average is taken over the controller parameter setX, not over the β values. The latter are participant specific and thus have to be changed according to the current participant in every simulation (see also Robustness Analysis and Discussion).When applying controller IIIc without the information about the disturbances (pi=0kg⋅m, IIIc.1) and additionally even without the explicit angular information (ki=01/deg⋅s2, IIIc.2), qualitatively nothing changes on the training set except for an increase in error by 100%. This means that, in principle, only angular velocity information is needed for stabilization and qualitative correct reproduction of the trials, all other information just enhance the fit and prediction quality.


#### 3.3.1. Training

We applied the average controller (parameters given in Table2), to the trials both in the training and the test set and evaluated the corresponding errors.The mean error ofIIIcavrg on Str was 1.35 deg/ts (assuming a mean trial length of 1.6 sec and thus 1600 timestepsts per angle making a total of 3200 ts for the trial). The greatest error made was on the trial (5)14200lt with an error of 2.17 deg/ts. Trial (3)10150lt was the best, with 0.87 deg/ts. The respective fits are shown in Figure 9.Figure 9
Best (a) and worst (b) prediction onStr. Plotted are the ankle and hip angle in deg. over time in seconds. The red/dashed lines are the experimental trials, and the blue/solid lines are the corresponding simulation of controller IIIcavrg. On the test set, the averaged controller ranged in quality of reproduction between those two trials. Obviously the prediction is good, since the test set trials were used to find the parameters X.
(a)
(b)


#### 3.3.2. Test

The average controller was tested onSte, the mean error on the test set was 2.24 deg/ts (again mean trial length = 1.6 sec), while trial (5)14150rt was the best fit with an error of 0.89 deg/ts, and the second best was (5)10150lt with 1.52 deg/ts and the greatest error was made on trial (3)18200lt with an error of 3.32 deg/ts, while the second greatest was 2.43 deg/ts in (3)18150lt. The average error without the worst two and the best two fits was 1.7 deg/ts. The best and worst fit are again shown in Figure 10.Figure 10
Best (a) and worst (b) prediction onSte. Plotted are the ankle and hip angle in deg. over time in seconds. The red/dashed lines are the experimental trials, and the blue/solid lines are the corresponding simulation of controller IIIcavrg. Both trials are qualitatively correct and the system is stabilized. The best prediction on the left shows a slight underestimation in the ankle and a slight overestimation in the hip at the end of trial, predicting the rest of the trajectory with high precision. The right (worst) prediction on the test set underestimates the experimental trajectory midtrial, while the beginning and end are predicted fine.
(a)
(b)


#### 3.3.3. Participant (4)

We furthermore compared the average controller’s performance to all nine trials of participant (4) (Ste4). While the time courses of the ankle angle fits are always reasonable, the great hip angle deflection cannot be reproduced by the controller fitted to the training set. The error in the least square step (calculating the ratios β) is already two magnitudes higher than with the other participants (usually between 20–200, for participant (4); however, errors are in the thousands). Also the error when applying the average controller is 1-2 magnitudes higher than for the other participants (mean error is 6.12 deg/ts, best fit is for trial (4)14150rt with 3.24 deg/ts, and worst fit is for (4)14150lt with 7.55 deg/ts). In comparison, the controller optimized to the specific trials of participant (4) has an mean error of 0.77 deg/ts, with the best fit for (4)14150rt with an error of 0.3 deg/ts and the worst fit for trial (4)18150lt with 1.14 deg/ts. Nevertheless, the average controller manages to reproduce the experiment trajectories qualitatively, which is shown in Figure 11, and the model restabilizes in response to the disturbance.Figure 11
Best (a) and worst (b) prediction for participant (4). Plotted are the ankle and hip angle in deg. over time in seconds. The red/dashed lines are the experimental trials, and the blue/solid lines are the corresponding simulation of controllerIIIcavrg. Without considering participant (4) in finding the control parameters X, the controller is able to predict his trials, however with a considerable amount of under and overestimation, but not to that point where qualitative change in behavior or instability is observed.
(a)
(b)


### 3.4. Robustness Analysis

Finally, we conducted a robustness analysis, by varying controller parameters of the controller IIIc optimized to a single trial (again (5)14200lt) manually, in order to assess the robustness of the simulated responses with respect to the parameters of the controller. Robustness in our context is defined as still showing qualitatively the same time course of the responses and stability (low remaining error between simulation and experiment at the end of balance recovery,t=tend).


#### 3.4.1. ParametersX

The controller IIIc has proven to be robust, not at last since it can cope with a range of perturbation strengths and varying kinematics in responses. Variations in the parametersk1,k2,q1,q2,delay1,delay2 of X more than 10% did not change qualitative behavior of the controller and still allowed for a good fit to the experimental curves to be made, even in the extreme case, where more than one parameter is varied by more than 10%. If only one parameter from k1,k2,q1,q2,delay1,delay2 is varied by 10%, the error goes up to about 0.55 deg/ts. By varying all parameters in X at once by ±10% (thus not changing the ratio between d1 and d2) the error grows from 0.18 deg/ts to 0.86 deg/ts, no loss of stability and no significant change in kinematics occur. The ratio between d1 and d2, that is, the ratio between the velocity contributions to both joints is the most critical one and should be varied by at most 5%; thus, if the ratio between the coupled torque generation of the controllers (which is mainly dependent on d1,d2) is changed, qualitative change and loss of stability is the consequence, for example, variation of d1 by 5% (for a fixed d2) leads to an error of 1.16 deg/ts, and a variation of 10% leads to qualitative change in behavior (especially in the hip) with an error of 2.4 deg/ts. For all trials, with a small enough remaining error (local minimum with an error ≤0.75deg/ts after optimization to a specific trial), the ratio was found to be the same, namely d1/d2=1.14±0.03. Thus, also if one varies all parameters of the ankle controller k1,d1,q1,delay1 by 10% without changing the hip controller parameters, the error goes to 2.62 deg/ts with qualitative change in response.Furthermore, in unconstrained nonlinear optimization, several (equivalent) local optima may be found(8)X∗=argminXEang=minX∑iΘA,iexp−ΘA,isim2+ΘH,iexp−ΘH,isim2,where i is the sampling points. So for example, optimization for (5)14200lt returned the set of parameters X1=0.721/deg⋅s2,2.041/deg⋅s,0.83kg⋅m,0.06s,0.661/deg⋅s2,1.71/deg⋅s,0.57kg⋅m,0.13s with a remaining error of 0.16 deg/ts and the second local minimum X2=1.401/deg⋅s2,1.791/deg⋅s,0.62kg⋅m,0.06s,1.231/deg⋅s2,1.571/deg⋅s,0.59kg⋅m,0.15s with a remaining error of 0.18 deg/ts, depending on initial conditions X0. However, note the ratio d1/d2=1.12 for X1 and 1.14 for X2 (please refer to Figure 8 to see the controller with X1), which is within the margin for both optimization results.As can be seen in Table2, the delays of the ankle were 60 ms in average (represented by delay1 of (2)) and the delays in the hip 171 ms in average (represented by delay2 of (2)). Even if the controller seems pretty robust in the delays (mentionable is that they were pretty robust if decreased but hit a limit very fast when increased at around +15%), all delays were found closely to these values in the optimization, with an at least twice as high delay in the hip as in the ankle. The balance compensatory responses in the hip are “reflex-like” responses. They are long latency responses that involved the central nervous system (CNS) and not only the spinal cord, as compared to the delays at the ankle joint.


#### 3.4.2. Ratiosβ

However, regarding the angular mixing ratiosβ (remember β is dependent on the participants’ individual physique, for example, inertia and geometry, see Table 3), the controller shows altered qualitative behavior with loss of stability if varied too much, and even for small variations, a remaining angular deviation from zero (or from the experimental curve) remains at the end of simulation (t=tend) which is too great. This was again a result firstly obtained from the construction of the average controller IIIcavrg, where in a first approach also the mean values of all β values in Str were used, while the DIP model still adapted segment data according to the current trial. This practice was leading to qualitatively wrong (and unstable) controller behavior (not fitting well to the experiment trajectory).Table 3
Participant specific ratiosβ.
Participant
β11
β12=β21
β22
β11/β12
β22/β21
(3)
1.34kg⋅m2
0.25kg⋅m2
0.06kg⋅m2
5.45
0.24
(4)
1.03kg⋅m2
0.20kg⋅m2
0.05kg⋅m2
5.06
0.26
(5)
1.12kg⋅m2
0.23kg⋅m2
0.06kg⋅m2
4.8
0.27Furthermore, we manually tried to determine robust ranges for the parametersβ or more precise for the ratios β11/β12 and β22/β21 and thus the ratios between mixing angles for the coupled torques, by varying the parameters β12=β21. For the test trial (5)14200lt, β12=β21 cannot exceed variations of ±2.5%, without undergoing qualitative changes (e.g., only positive/negative angles) and for variations greater than ±5% loss of stability. Analogous outcome for variation of β11 or β22. This is reflected in the errors for varying β11 or β22: for the original parameters, it is 0.18 deg/ts and for a variation of 2.5% it is 0.68 deg/ts, exhibiting slight under and over fitting of the angular trajectories. For a variation of 5%, it is 1.18 deg/ts with a remaining error at the end of balancing, and for a variation of 10%, it is 2.5 deg/ts, with loss of stability. If all values in β are changed by the same ±10%, the error grows to 0.87 deg/ts. Thus, the controller still works, if the ratios of angular mixing are left untouched (which are thus mainly depended on β and d1/2 in X).Remember, the error values alone are not sufficient to detect qualitative change. Only by looking at the trajectories one may determine if qualitative change has indeed occurred. The error is merely an indicator, which works good to determine if a fit or prediction was correct (low error values) but not if it was incorrect just because the error values are high(er).


## 3.1. Experiment

### 3.1.1. Perturbation

A set of 12 perturbations was executed in the experiment. They differ in magnitude and direction but not qualitatively. In Figure3, the motion signal used for perturbing the DIP model is shown for a test case. The test case is a trial by participant (5) (which is 14cm displacement and 200cm/sec2 acceleration trial; from now on referred to as (5)14200lt), as this is closest to the average in the experiment. Trial (5)14200lt is shown in Figure 4 as an example. The whole motion, from preperturbation (A), to perturbation (B, C), to initial participant response (D), to full balance recovery (E, F), is explained in the following.Figure 3
Filtered perturbation signal. For each trial, a qualitatively similar signal, containing change in position, velocity, and acceleration, was used. Quantitatively the perturbation signals differed in direction, thus to the left or right and in displacement and velocity.Figure 4
(a) Typical frames of the experiment, from a single trial, at the distinguished phases (A, B, C, D, E, F). (b) Corresponding trajectories of the hip and ankle angle, horizontal lines correspond to the phases. See also AppendixD for further reading.
(a)
(b)


### 3.1.2. Results of the Experiment and Strategy-Dependent Variations

A simple segment length criterion (see AppendixA.2) identified two participants (1, 2), which deviated from the basic strategy, invalidating almost all their trials due to length changes of up to 15cm in the Leg segment and up to 6cm in the HAT segment. This was confirmed visually by the video, where in nearly all trials, participant (1) is extremely bending the knees and rotating the hip and spine horizontally, resulting in heavy A/P tilting. Participant (2) on the other hand acts extreme stiff around the lumbar spine and hip, resulting in a lift of one or both feet, see also Section 4.2.Loss of balance (i.e., falling) indeed never occurred, as participants used arm raising or stepping if needed, resulting in71/87 (81.6%) initially marked as “successful trials” (no arm movement and intentional stepping detected). Aside from that, within the participants (3, 4, 5) exhibiting a “proper” execution of the strategy (e.g., no A/P movement), only 5/36 trials (13.9%) had to be excluded.Participant (4), as can also be seen in Figure5 exhibited bending of the lumbar spine (just like participant (1) and in contrast to (2, 3, 5)). However in his case, this led to changes of the segment length within the error margin r, and thus, his trials were not excluded. Here, the hip angle excursions were up to 3.5 times higher than those of participants (3 and 5).Figure 5
Unfiltered angular trajectories for valid trials of the experiment. In the (a) plot are the perturbations to the left, and in the (b) are perturbations to the right. The solid line (blue) is the hip the dashed line (red) is the ankle angle trajectory. All the “extreme” hip angle peaks stem from participant (4).
(a)
(b)Within the valid trials, the ankle bends to a peak 5–7.5° while the “extended hip” (i.e., hip + lower spine) bends to a peak of up to 10° (34° for participant (4)). The higher variety within the hip stems from bending of the lumbar spine, which was accepted if segment lengths did not change significantly and thus lumped into the “hip” DOF. However, bending of the spine in the coronal plane, together with bending of the hip in the sagittal plane (often together with a rotation in the transversal plane), occurs in some trials, in all participants and thus, leading to greater segment length changes. Such trials were excluded from our data set as our planar model could not cope with them. In Figure5, valid trials of the experiment are plotted, while in Table 1, the mean of all max angles and angular velocities over all valid trials is given for each participant.Table 1
Mean and standard deviation of the maximal absolute angle and angular velocity for the ankle and hip for all valid trials.
Participant
Θamax
Θhmax
Θ̇amax
Θ̇hmax
(3)
5.26±1.94deg
10.38±5.42deg
19.88±7.09deg/s
60.29±26.98deg/s
(4)
3.73±1.55deg
22.40±8.76deg
24.00±8.73deg/s
106.7±39.45deg/s
(5)
5.62±0.9deg
7.45±2.57deg
25.94±4.08deg/s
50.14±10.42deg/s


## 3.1.1. Perturbation

A set of 12 perturbations was executed in the experiment. They differ in magnitude and direction but not qualitatively. In Figure3, the motion signal used for perturbing the DIP model is shown for a test case. The test case is a trial by participant (5) (which is 14cm displacement and 200cm/sec2 acceleration trial; from now on referred to as (5)14200lt), as this is closest to the average in the experiment. Trial (5)14200lt is shown in Figure 4 as an example. The whole motion, from preperturbation (A), to perturbation (B, C), to initial participant response (D), to full balance recovery (E, F), is explained in the following.Figure 3
Filtered perturbation signal. For each trial, a qualitatively similar signal, containing change in position, velocity, and acceleration, was used. Quantitatively the perturbation signals differed in direction, thus to the left or right and in displacement and velocity.Figure 4
(a) Typical frames of the experiment, from a single trial, at the distinguished phases (A, B, C, D, E, F). (b) Corresponding trajectories of the hip and ankle angle, horizontal lines correspond to the phases. See also AppendixD for further reading.
(a)
(b)


## 3.1.2. Results of the Experiment and Strategy-Dependent Variations

A simple segment length criterion (see AppendixA.2) identified two participants (1, 2), which deviated from the basic strategy, invalidating almost all their trials due to length changes of up to 15cm in the Leg segment and up to 6cm in the HAT segment. This was confirmed visually by the video, where in nearly all trials, participant (1) is extremely bending the knees and rotating the hip and spine horizontally, resulting in heavy A/P tilting. Participant (2) on the other hand acts extreme stiff around the lumbar spine and hip, resulting in a lift of one or both feet, see also Section 4.2.Loss of balance (i.e., falling) indeed never occurred, as participants used arm raising or stepping if needed, resulting in71/87 (81.6%) initially marked as “successful trials” (no arm movement and intentional stepping detected). Aside from that, within the participants (3, 4, 5) exhibiting a “proper” execution of the strategy (e.g., no A/P movement), only 5/36 trials (13.9%) had to be excluded.Participant (4), as can also be seen in Figure5 exhibited bending of the lumbar spine (just like participant (1) and in contrast to (2, 3, 5)). However in his case, this led to changes of the segment length within the error margin r, and thus, his trials were not excluded. Here, the hip angle excursions were up to 3.5 times higher than those of participants (3 and 5).Figure 5
Unfiltered angular trajectories for valid trials of the experiment. In the (a) plot are the perturbations to the left, and in the (b) are perturbations to the right. The solid line (blue) is the hip the dashed line (red) is the ankle angle trajectory. All the “extreme” hip angle peaks stem from participant (4).
(a)
(b)Within the valid trials, the ankle bends to a peak 5–7.5° while the “extended hip” (i.e., hip + lower spine) bends to a peak of up to 10° (34° for participant (4)). The higher variety within the hip stems from bending of the lumbar spine, which was accepted if segment lengths did not change significantly and thus lumped into the “hip” DOF. However, bending of the spine in the coronal plane, together with bending of the hip in the sagittal plane (often together with a rotation in the transversal plane), occurs in some trials, in all participants and thus, leading to greater segment length changes. Such trials were excluded from our data set as our planar model could not cope with them. In Figure5, valid trials of the experiment are plotted, while in Table 1, the mean of all max angles and angular velocities over all valid trials is given for each participant.Table 1
Mean and standard deviation of the maximal absolute angle and angular velocity for the ankle and hip for all valid trials.
Participant
Θamax
Θhmax
Θ̇amax
Θ̇hmax
(3)
5.26±1.94deg
10.38±5.42deg
19.88±7.09deg/s
60.29±26.98deg/s
(4)
3.73±1.55deg
22.40±8.76deg
24.00±8.73deg/s
106.7±39.45deg/s
(5)
5.62±0.9deg
7.45±2.57deg
25.94±4.08deg/s
50.14±10.42deg/s


## 3.2. Model and Control

### 3.2.1. Reference Controller (I) and (IIa/b)

In Figure6, we are presenting a comparison of results obtained by using “no control” (I), “P(200)D(20)” (IIa), and “P(200)D(20) + 20 ms delay” (IIb), to the results obtained in the experiment. While, obviously, none of the approaches can reproduce the experimental data, it is still possible to stabilize the system using a common PD controller (which is making the model behave very similar to a SIP), however, only if there is no delay introduced in the sensor signals or the controller output. Note that delay on either the output or delay on the sensor signals leads to a similar result (Figure 6(b)) and that delays are inevitable in humans. Consequently, any simple PD control must fail to stabilize in the biologically restricted system.Figure 6
(a) Experiment trajectory in blue, “no control” (I) in red, and PD control (IIa) in dashed red. (b) Experiment trajectory in blue and PD control with 20 ms delay (IIb) in dashed red. First plot is the ankle trajectory, second is the hip, and third is the perturbation. While with no control applied, the system obviously collapses, also simple PD control, while being able to stabilize was not able to match the angular trajectories of the experiment. If a delay is introduced, the system becomes highly unstable. The PD controller used here is the MATLAB : Simulink build in PID controller withP=200, D=20, I=0, and N=100.
(a)
(b)


### 3.2.2. Controller (IIIa/b)

Figure7 depicts the best fits of the first two approaches (IIIa, b) in comparison to the test case (5)14200lt. The parameters of the dynamic system (IIIa/b) are optimized such that Θa/hsimt are fitted to Θa/hexpt of a single experimental trial ((5)14200lt), which is done by unconstrained nonlinear optimization.Figure 7
(a) Controller IIIa, best fit, with an remaining error of 45.73deg/ts. (b) Controller IIIb, best fit, with an remaining error of 2.7 deg/ts. Both fitted to the same experimental trial (5)14200lt. First plot is the ankle trajectory, second is the hip, and third is the perturbation. While IIIa fails to stabilize (no torque coupling), even the introduction of simple one to one torque coupling (IIIb) is enough to (1) stabilize the system and (2) almost reproduce the experimental trajectories.
(a)
(b)The approach (IIIa) was invalidated, as it completely fails to reproduce the angular trajectories (not stabilizing the mechanical model at all).The approach (IIIb), which, by using a linear combination of both angles, is a lot closer to the sought solution, is nevertheless also invalid: the remaining error after optimization is big, and the controller is unable to stabilize the mechanical model.


### 3.2.3. Controller (IIIc)

First, the weights fory1/2 (see (4c)) are calculated by an inverse dynamics approach via least square optimization. They are specific to the participant (not to the perturbation), that is, β5 is the same for all trials of participant (5), within small numerical differences. Thus, the weights can be calculated a priori, individually for each participant. The weights of participant (5) are β11=1.12,β12=β21=0.23,β22=0.06 (kg⋅m2). Note that it is a result of the simulation that β12=β21, within an accuracy 0.1%. The before mentioned unconstrained nonlinear optimization algorithm is used to find well-suited control parameters X=k1,d1,q1,delay1,k2,d2,q2,delay2.“Optimal” (i.e., local minimum) values for ((5)14200lt) are(7)k1=0.721deg⋅s2,d1=2.041deg⋅s,q1=0.83kg⋅m,delay1=0.06s,k2=0.661deg⋅s2,d2=1.71deg⋅s,q2=0.57kg⋅m,delay2=0.13s.The left plot in Figure8 depicts the “exact” reproduction of Mahg and Mhhg by the LSQ algorithm (used to find β, remember (5)) and the trajectories of the corresponding controller, utilizing the upper 2×2 part of β, in comparison to the experiments on the right.Figure 8
(a) The torques generated by the high gain PD:Mihg versus MiLSQ, i∈a,h with virtually MLSQ=Msim. y-axis is torque in Nm over time steps of the simulation on the x-axis. (b) Response of controller IIIc. The best fit, with an remaining error of 0.16 deg/ts. First plot is the ankle trajectory, second is the hip, and third is the perturbation. Finally, the controller IIIc with torque coupling according to the numerical inertia tensor β realized the almost exact reproduction of the experiments.
(a)
(b)


## 3.2.1. Reference Controller (I) and (IIa/b)

In Figure6, we are presenting a comparison of results obtained by using “no control” (I), “P(200)D(20)” (IIa), and “P(200)D(20) + 20 ms delay” (IIb), to the results obtained in the experiment. While, obviously, none of the approaches can reproduce the experimental data, it is still possible to stabilize the system using a common PD controller (which is making the model behave very similar to a SIP), however, only if there is no delay introduced in the sensor signals or the controller output. Note that delay on either the output or delay on the sensor signals leads to a similar result (Figure 6(b)) and that delays are inevitable in humans. Consequently, any simple PD control must fail to stabilize in the biologically restricted system.Figure 6
(a) Experiment trajectory in blue, “no control” (I) in red, and PD control (IIa) in dashed red. (b) Experiment trajectory in blue and PD control with 20 ms delay (IIb) in dashed red. First plot is the ankle trajectory, second is the hip, and third is the perturbation. While with no control applied, the system obviously collapses, also simple PD control, while being able to stabilize was not able to match the angular trajectories of the experiment. If a delay is introduced, the system becomes highly unstable. The PD controller used here is the MATLAB : Simulink build in PID controller withP=200, D=20, I=0, and N=100.
(a)
(b)


## 3.2.2. Controller (IIIa/b)

Figure7 depicts the best fits of the first two approaches (IIIa, b) in comparison to the test case (5)14200lt. The parameters of the dynamic system (IIIa/b) are optimized such that Θa/hsimt are fitted to Θa/hexpt of a single experimental trial ((5)14200lt), which is done by unconstrained nonlinear optimization.Figure 7
(a) Controller IIIa, best fit, with an remaining error of 45.73deg/ts. (b) Controller IIIb, best fit, with an remaining error of 2.7 deg/ts. Both fitted to the same experimental trial (5)14200lt. First plot is the ankle trajectory, second is the hip, and third is the perturbation. While IIIa fails to stabilize (no torque coupling), even the introduction of simple one to one torque coupling (IIIb) is enough to (1) stabilize the system and (2) almost reproduce the experimental trajectories.
(a)
(b)The approach (IIIa) was invalidated, as it completely fails to reproduce the angular trajectories (not stabilizing the mechanical model at all).The approach (IIIb), which, by using a linear combination of both angles, is a lot closer to the sought solution, is nevertheless also invalid: the remaining error after optimization is big, and the controller is unable to stabilize the mechanical model.


## 3.2.3. Controller (IIIc)

First, the weights fory1/2 (see (4c)) are calculated by an inverse dynamics approach via least square optimization. They are specific to the participant (not to the perturbation), that is, β5 is the same for all trials of participant (5), within small numerical differences. Thus, the weights can be calculated a priori, individually for each participant. The weights of participant (5) are β11=1.12,β12=β21=0.23,β22=0.06 (kg⋅m2). Note that it is a result of the simulation that β12=β21, within an accuracy 0.1%. The before mentioned unconstrained nonlinear optimization algorithm is used to find well-suited control parameters X=k1,d1,q1,delay1,k2,d2,q2,delay2.“Optimal” (i.e., local minimum) values for ((5)14200lt) are(7)k1=0.721deg⋅s2,d1=2.041deg⋅s,q1=0.83kg⋅m,delay1=0.06s,k2=0.661deg⋅s2,d2=1.71deg⋅s,q2=0.57kg⋅m,delay2=0.13s.The left plot in Figure8 depicts the “exact” reproduction of Mahg and Mhhg by the LSQ algorithm (used to find β, remember (5)) and the trajectories of the corresponding controller, utilizing the upper 2×2 part of β, in comparison to the experiments on the right.Figure 8
(a) The torques generated by the high gain PD:Mihg versus MiLSQ, i∈a,h with virtually MLSQ=Msim. y-axis is torque in Nm over time steps of the simulation on the x-axis. (b) Response of controller IIIc. The best fit, with an remaining error of 0.16 deg/ts. First plot is the ankle trajectory, second is the hip, and third is the perturbation. Finally, the controller IIIc with torque coupling according to the numerical inertia tensor β realized the almost exact reproduction of the experiments.
(a)
(b)


## 3.3. Average Controller “IIIcavrg”

Of the 31 valid trials, we took 10 for the training setStr and 11 for the test set Ste. Participant (4) is not part of training and test set due to his great (quantitative) deviation in strategy (extreme hip angles by bending the lower spine a lot more than other participants); he was treated separately by checking the performance, by the set Str parametrized average controller, on all his trials.In Table2, the average of all parameters after optimization on the training set is listed. The column “Error” is the average error made by the controller optimized for each specific trial in Str. This value is useful comparing the average controller error on the training and test set to the error made by the controller optimized to specific trials. The mean error thus is 0.4 deg/ts, while the best fit is trial (5)14200lt with 0.18 deg/ts and the worst fit is trial (3)14150lt with 0.58 deg/ts.Table 2
Average parameter setX for controller IIIc on the training set Str referred to by IIIcavrg.
Controller
Error
k1
d1
q1
delay1
k2
d2
q2
delay2
IIIc
0.35deg/s
1.181/deg⋅s2
1.971/deg⋅s2
0.45kg⋅m
0.06 s
0.851/deg⋅s2
1.731/deg⋅s
0.53kg⋅m
0.17 s
IIIc.1
0.54deg/ts
1.271/deg⋅s2
2.531/deg⋅s
—
0.06 s
1.191/deg⋅s2
2.061/deg⋅s
—
0.12 s
IIIc.2
0.55deg/ts
—
2.531/deg⋅s
—
0.06 s
—
2.051/deg⋅s
—
0.12 sNote that the average is taken over the controller parameter setX, not over the β values. The latter are participant specific and thus have to be changed according to the current participant in every simulation (see also Robustness Analysis and Discussion).When applying controller IIIc without the information about the disturbances (pi=0kg⋅m, IIIc.1) and additionally even without the explicit angular information (ki=01/deg⋅s2, IIIc.2), qualitatively nothing changes on the training set except for an increase in error by 100%. This means that, in principle, only angular velocity information is needed for stabilization and qualitative correct reproduction of the trials, all other information just enhance the fit and prediction quality.


### 3.3.1. Training

We applied the average controller (parameters given in Table2), to the trials both in the training and the test set and evaluated the corresponding errors.The mean error ofIIIcavrg on Str was 1.35 deg/ts (assuming a mean trial length of 1.6 sec and thus 1600 timestepsts per angle making a total of 3200 ts for the trial). The greatest error made was on the trial (5)14200lt with an error of 2.17 deg/ts. Trial (3)10150lt was the best, with 0.87 deg/ts. The respective fits are shown in Figure 9.Figure 9
Best (a) and worst (b) prediction onStr. Plotted are the ankle and hip angle in deg. over time in seconds. The red/dashed lines are the experimental trials, and the blue/solid lines are the corresponding simulation of controller IIIcavrg. On the test set, the averaged controller ranged in quality of reproduction between those two trials. Obviously the prediction is good, since the test set trials were used to find the parameters X.
(a)
(b)


### 3.3.2. Test

The average controller was tested onSte, the mean error on the test set was 2.24 deg/ts (again mean trial length = 1.6 sec), while trial (5)14150rt was the best fit with an error of 0.89 deg/ts, and the second best was (5)10150lt with 1.52 deg/ts and the greatest error was made on trial (3)18200lt with an error of 3.32 deg/ts, while the second greatest was 2.43 deg/ts in (3)18150lt. The average error without the worst two and the best two fits was 1.7 deg/ts. The best and worst fit are again shown in Figure 10.Figure 10
Best (a) and worst (b) prediction onSte. Plotted are the ankle and hip angle in deg. over time in seconds. The red/dashed lines are the experimental trials, and the blue/solid lines are the corresponding simulation of controller IIIcavrg. Both trials are qualitatively correct and the system is stabilized. The best prediction on the left shows a slight underestimation in the ankle and a slight overestimation in the hip at the end of trial, predicting the rest of the trajectory with high precision. The right (worst) prediction on the test set underestimates the experimental trajectory midtrial, while the beginning and end are predicted fine.
(a)
(b)


### 3.3.3. Participant (4)

We furthermore compared the average controller’s performance to all nine trials of participant (4) (Ste4). While the time courses of the ankle angle fits are always reasonable, the great hip angle deflection cannot be reproduced by the controller fitted to the training set. The error in the least square step (calculating the ratios β) is already two magnitudes higher than with the other participants (usually between 20–200, for participant (4); however, errors are in the thousands). Also the error when applying the average controller is 1-2 magnitudes higher than for the other participants (mean error is 6.12 deg/ts, best fit is for trial (4)14150rt with 3.24 deg/ts, and worst fit is for (4)14150lt with 7.55 deg/ts). In comparison, the controller optimized to the specific trials of participant (4) has an mean error of 0.77 deg/ts, with the best fit for (4)14150rt with an error of 0.3 deg/ts and the worst fit for trial (4)18150lt with 1.14 deg/ts. Nevertheless, the average controller manages to reproduce the experiment trajectories qualitatively, which is shown in Figure 11, and the model restabilizes in response to the disturbance.Figure 11
Best (a) and worst (b) prediction for participant (4). Plotted are the ankle and hip angle in deg. over time in seconds. The red/dashed lines are the experimental trials, and the blue/solid lines are the corresponding simulation of controllerIIIcavrg. Without considering participant (4) in finding the control parameters X, the controller is able to predict his trials, however with a considerable amount of under and overestimation, but not to that point where qualitative change in behavior or instability is observed.
(a)
(b)


## 3.3.1. Training

We applied the average controller (parameters given in Table2), to the trials both in the training and the test set and evaluated the corresponding errors.The mean error ofIIIcavrg on Str was 1.35 deg/ts (assuming a mean trial length of 1.6 sec and thus 1600 timestepsts per angle making a total of 3200 ts for the trial). The greatest error made was on the trial (5)14200lt with an error of 2.17 deg/ts. Trial (3)10150lt was the best, with 0.87 deg/ts. The respective fits are shown in Figure 9.Figure 9
Best (a) and worst (b) prediction onStr. Plotted are the ankle and hip angle in deg. over time in seconds. The red/dashed lines are the experimental trials, and the blue/solid lines are the corresponding simulation of controller IIIcavrg. On the test set, the averaged controller ranged in quality of reproduction between those two trials. Obviously the prediction is good, since the test set trials were used to find the parameters X.
(a)
(b)


## 3.3.2. Test

The average controller was tested onSte, the mean error on the test set was 2.24 deg/ts (again mean trial length = 1.6 sec), while trial (5)14150rt was the best fit with an error of 0.89 deg/ts, and the second best was (5)10150lt with 1.52 deg/ts and the greatest error was made on trial (3)18200lt with an error of 3.32 deg/ts, while the second greatest was 2.43 deg/ts in (3)18150lt. The average error without the worst two and the best two fits was 1.7 deg/ts. The best and worst fit are again shown in Figure 10.Figure 10
Best (a) and worst (b) prediction onSte. Plotted are the ankle and hip angle in deg. over time in seconds. The red/dashed lines are the experimental trials, and the blue/solid lines are the corresponding simulation of controller IIIcavrg. Both trials are qualitatively correct and the system is stabilized. The best prediction on the left shows a slight underestimation in the ankle and a slight overestimation in the hip at the end of trial, predicting the rest of the trajectory with high precision. The right (worst) prediction on the test set underestimates the experimental trajectory midtrial, while the beginning and end are predicted fine.
(a)
(b)


## 3.3.3. Participant (4)

We furthermore compared the average controller’s performance to all nine trials of participant (4) (Ste4). While the time courses of the ankle angle fits are always reasonable, the great hip angle deflection cannot be reproduced by the controller fitted to the training set. The error in the least square step (calculating the ratios β) is already two magnitudes higher than with the other participants (usually between 20–200, for participant (4); however, errors are in the thousands). Also the error when applying the average controller is 1-2 magnitudes higher than for the other participants (mean error is 6.12 deg/ts, best fit is for trial (4)14150rt with 3.24 deg/ts, and worst fit is for (4)14150lt with 7.55 deg/ts). In comparison, the controller optimized to the specific trials of participant (4) has an mean error of 0.77 deg/ts, with the best fit for (4)14150rt with an error of 0.3 deg/ts and the worst fit for trial (4)18150lt with 1.14 deg/ts. Nevertheless, the average controller manages to reproduce the experiment trajectories qualitatively, which is shown in Figure 11, and the model restabilizes in response to the disturbance.Figure 11
Best (a) and worst (b) prediction for participant (4). Plotted are the ankle and hip angle in deg. over time in seconds. The red/dashed lines are the experimental trials, and the blue/solid lines are the corresponding simulation of controllerIIIcavrg. Without considering participant (4) in finding the control parameters X, the controller is able to predict his trials, however with a considerable amount of under and overestimation, but not to that point where qualitative change in behavior or instability is observed.
(a)
(b)


## 3.4. Robustness Analysis

Finally, we conducted a robustness analysis, by varying controller parameters of the controller IIIc optimized to a single trial (again (5)14200lt) manually, in order to assess the robustness of the simulated responses with respect to the parameters of the controller. Robustness in our context is defined as still showing qualitatively the same time course of the responses and stability (low remaining error between simulation and experiment at the end of balance recovery,t=tend).


### 3.4.1. ParametersX

The controller IIIc has proven to be robust, not at last since it can cope with a range of perturbation strengths and varying kinematics in responses. Variations in the parametersk1,k2,q1,q2,delay1,delay2 of X more than 10% did not change qualitative behavior of the controller and still allowed for a good fit to the experimental curves to be made, even in the extreme case, where more than one parameter is varied by more than 10%. If only one parameter from k1,k2,q1,q2,delay1,delay2 is varied by 10%, the error goes up to about 0.55 deg/ts. By varying all parameters in X at once by ±10% (thus not changing the ratio between d1 and d2) the error grows from 0.18 deg/ts to 0.86 deg/ts, no loss of stability and no significant change in kinematics occur. The ratio between d1 and d2, that is, the ratio between the velocity contributions to both joints is the most critical one and should be varied by at most 5%; thus, if the ratio between the coupled torque generation of the controllers (which is mainly dependent on d1,d2) is changed, qualitative change and loss of stability is the consequence, for example, variation of d1 by 5% (for a fixed d2) leads to an error of 1.16 deg/ts, and a variation of 10% leads to qualitative change in behavior (especially in the hip) with an error of 2.4 deg/ts. For all trials, with a small enough remaining error (local minimum with an error ≤0.75deg/ts after optimization to a specific trial), the ratio was found to be the same, namely d1/d2=1.14±0.03. Thus, also if one varies all parameters of the ankle controller k1,d1,q1,delay1 by 10% without changing the hip controller parameters, the error goes to 2.62 deg/ts with qualitative change in response.Furthermore, in unconstrained nonlinear optimization, several (equivalent) local optima may be found(8)X∗=argminXEang=minX∑iΘA,iexp−ΘA,isim2+ΘH,iexp−ΘH,isim2,where i is the sampling points. So for example, optimization for (5)14200lt returned the set of parameters X1=0.721/deg⋅s2,2.041/deg⋅s,0.83kg⋅m,0.06s,0.661/deg⋅s2,1.71/deg⋅s,0.57kg⋅m,0.13s with a remaining error of 0.16 deg/ts and the second local minimum X2=1.401/deg⋅s2,1.791/deg⋅s,0.62kg⋅m,0.06s,1.231/deg⋅s2,1.571/deg⋅s,0.59kg⋅m,0.15s with a remaining error of 0.18 deg/ts, depending on initial conditions X0. However, note the ratio d1/d2=1.12 for X1 and 1.14 for X2 (please refer to Figure 8 to see the controller with X1), which is within the margin for both optimization results.As can be seen in Table2, the delays of the ankle were 60 ms in average (represented by delay1 of (2)) and the delays in the hip 171 ms in average (represented by delay2 of (2)). Even if the controller seems pretty robust in the delays (mentionable is that they were pretty robust if decreased but hit a limit very fast when increased at around +15%), all delays were found closely to these values in the optimization, with an at least twice as high delay in the hip as in the ankle. The balance compensatory responses in the hip are “reflex-like” responses. They are long latency responses that involved the central nervous system (CNS) and not only the spinal cord, as compared to the delays at the ankle joint.


### 3.4.2. Ratiosβ

However, regarding the angular mixing ratiosβ (remember β is dependent on the participants’ individual physique, for example, inertia and geometry, see Table 3), the controller shows altered qualitative behavior with loss of stability if varied too much, and even for small variations, a remaining angular deviation from zero (or from the experimental curve) remains at the end of simulation (t=tend) which is too great. This was again a result firstly obtained from the construction of the average controller IIIcavrg, where in a first approach also the mean values of all β values in Str were used, while the DIP model still adapted segment data according to the current trial. This practice was leading to qualitatively wrong (and unstable) controller behavior (not fitting well to the experiment trajectory).Table 3
Participant specific ratiosβ.
Participant
β11
β12=β21
β22
β11/β12
β22/β21
(3)
1.34kg⋅m2
0.25kg⋅m2
0.06kg⋅m2
5.45
0.24
(4)
1.03kg⋅m2
0.20kg⋅m2
0.05kg⋅m2
5.06
0.26
(5)
1.12kg⋅m2
0.23kg⋅m2
0.06kg⋅m2
4.8
0.27Furthermore, we manually tried to determine robust ranges for the parametersβ or more precise for the ratios β11/β12 and β22/β21 and thus the ratios between mixing angles for the coupled torques, by varying the parameters β12=β21. For the test trial (5)14200lt, β12=β21 cannot exceed variations of ±2.5%, without undergoing qualitative changes (e.g., only positive/negative angles) and for variations greater than ±5% loss of stability. Analogous outcome for variation of β11 or β22. This is reflected in the errors for varying β11 or β22: for the original parameters, it is 0.18 deg/ts and for a variation of 2.5% it is 0.68 deg/ts, exhibiting slight under and over fitting of the angular trajectories. For a variation of 5%, it is 1.18 deg/ts with a remaining error at the end of balancing, and for a variation of 10%, it is 2.5 deg/ts, with loss of stability. If all values in β are changed by the same ±10%, the error grows to 0.87 deg/ts. Thus, the controller still works, if the ratios of angular mixing are left untouched (which are thus mainly depended on β and d1/2 in X).Remember, the error values alone are not sufficient to detect qualitative change. Only by looking at the trajectories one may determine if qualitative change has indeed occurred. The error is merely an indicator, which works good to determine if a fit or prediction was correct (low error values) but not if it was incorrect just because the error values are high(er).


## 3.4.1. ParametersX

The controller IIIc has proven to be robust, not at last since it can cope with a range of perturbation strengths and varying kinematics in responses. Variations in the parametersk1,k2,q1,q2,delay1,delay2 of X more than 10% did not change qualitative behavior of the controller and still allowed for a good fit to the experimental curves to be made, even in the extreme case, where more than one parameter is varied by more than 10%. If only one parameter from k1,k2,q1,q2,delay1,delay2 is varied by 10%, the error goes up to about 0.55 deg/ts. By varying all parameters in X at once by ±10% (thus not changing the ratio between d1 and d2) the error grows from 0.18 deg/ts to 0.86 deg/ts, no loss of stability and no significant change in kinematics occur. The ratio between d1 and d2, that is, the ratio between the velocity contributions to both joints is the most critical one and should be varied by at most 5%; thus, if the ratio between the coupled torque generation of the controllers (which is mainly dependent on d1,d2) is changed, qualitative change and loss of stability is the consequence, for example, variation of d1 by 5% (for a fixed d2) leads to an error of 1.16 deg/ts, and a variation of 10% leads to qualitative change in behavior (especially in the hip) with an error of 2.4 deg/ts. For all trials, with a small enough remaining error (local minimum with an error ≤0.75deg/ts after optimization to a specific trial), the ratio was found to be the same, namely d1/d2=1.14±0.03. Thus, also if one varies all parameters of the ankle controller k1,d1,q1,delay1 by 10% without changing the hip controller parameters, the error goes to 2.62 deg/ts with qualitative change in response.Furthermore, in unconstrained nonlinear optimization, several (equivalent) local optima may be found(8)X∗=argminXEang=minX∑iΘA,iexp−ΘA,isim2+ΘH,iexp−ΘH,isim2,where i is the sampling points. So for example, optimization for (5)14200lt returned the set of parameters X1=0.721/deg⋅s2,2.041/deg⋅s,0.83kg⋅m,0.06s,0.661/deg⋅s2,1.71/deg⋅s,0.57kg⋅m,0.13s with a remaining error of 0.16 deg/ts and the second local minimum X2=1.401/deg⋅s2,1.791/deg⋅s,0.62kg⋅m,0.06s,1.231/deg⋅s2,1.571/deg⋅s,0.59kg⋅m,0.15s with a remaining error of 0.18 deg/ts, depending on initial conditions X0. However, note the ratio d1/d2=1.12 for X1 and 1.14 for X2 (please refer to Figure 8 to see the controller with X1), which is within the margin for both optimization results.As can be seen in Table2, the delays of the ankle were 60 ms in average (represented by delay1 of (2)) and the delays in the hip 171 ms in average (represented by delay2 of (2)). Even if the controller seems pretty robust in the delays (mentionable is that they were pretty robust if decreased but hit a limit very fast when increased at around +15%), all delays were found closely to these values in the optimization, with an at least twice as high delay in the hip as in the ankle. The balance compensatory responses in the hip are “reflex-like” responses. They are long latency responses that involved the central nervous system (CNS) and not only the spinal cord, as compared to the delays at the ankle joint.


## 3.4.2. Ratiosβ

However, regarding the angular mixing ratiosβ (remember β is dependent on the participants’ individual physique, for example, inertia and geometry, see Table 3), the controller shows altered qualitative behavior with loss of stability if varied too much, and even for small variations, a remaining angular deviation from zero (or from the experimental curve) remains at the end of simulation (t=tend) which is too great. This was again a result firstly obtained from the construction of the average controller IIIcavrg, where in a first approach also the mean values of all β values in Str were used, while the DIP model still adapted segment data according to the current trial. This practice was leading to qualitatively wrong (and unstable) controller behavior (not fitting well to the experiment trajectory).Table 3
Participant specific ratiosβ.
Participant
β11
β12=β21
β22
β11/β12
β22/β21
(3)
1.34kg⋅m2
0.25kg⋅m2
0.06kg⋅m2
5.45
0.24
(4)
1.03kg⋅m2
0.20kg⋅m2
0.05kg⋅m2
5.06
0.26
(5)
1.12kg⋅m2
0.23kg⋅m2
0.06kg⋅m2
4.8
0.27Furthermore, we manually tried to determine robust ranges for the parametersβ or more precise for the ratios β11/β12 and β22/β21 and thus the ratios between mixing angles for the coupled torques, by varying the parameters β12=β21. For the test trial (5)14200lt, β12=β21 cannot exceed variations of ±2.5%, without undergoing qualitative changes (e.g., only positive/negative angles) and for variations greater than ±5% loss of stability. Analogous outcome for variation of β11 or β22. This is reflected in the errors for varying β11 or β22: for the original parameters, it is 0.18 deg/ts and for a variation of 2.5% it is 0.68 deg/ts, exhibiting slight under and over fitting of the angular trajectories. For a variation of 5%, it is 1.18 deg/ts with a remaining error at the end of balancing, and for a variation of 10%, it is 2.5 deg/ts, with loss of stability. If all values in β are changed by the same ±10%, the error grows to 0.87 deg/ts. Thus, the controller still works, if the ratios of angular mixing are left untouched (which are thus mainly depended on β and d1/2 in X).Remember, the error values alone are not sufficient to detect qualitative change. Only by looking at the trajectories one may determine if qualitative change has indeed occurred. The error is merely an indicator, which works good to determine if a fit or prediction was correct (low error values) but not if it was incorrect just because the error values are high(er).


## 4. Discussion

### 4.1. Generality of the Controller

The controller IIIc is not only used to reproduce trials (to which it was parametrized to), but also to predict trials not used in training (remember training error of 1.35deg/ts). The quality of prediction for this test set is absolutely within the range of the fitting to the training set (even if the numerical error value on the test set with 2.24 deg/ts or 1.7 deg/ts is higher than in training, the fits are very good as can be seen, comparing Figures 9 and 10).Notable is that the training consisted of trials from two participants, which were used regardless of their magnitude, over the full range of perturbations tested in the experiment, also ignoring the direction. As this was not only possible, but also leading to good results, the strong correlation between trials of different magnitudes and participants could be verified: the same control scheme, with one parameter set could be used. The full range of reactions to the perturbation, tolerable by the controller, is also illustrated by Figure5. We are certain that for a greater number of trials, participants, and magnitudes, this mean result would converge to a smaller error than possible now, under the prior mentioned premises, and to a certainly much more accurate prediction if only one magnitude of perturbations is used with many participants. It would also be possible to extend the found control law such that it considers the perturbation strength to a higher degree and apply different magnitudes of reaction accordingly.However, we have chosen this advance to prove some generality of the controller (e.g., robustness over a broad range of trials). Even for participant (4), where a variation of strategy could be observed (much greater angles in the hip due to extended spinal bending in l1-l5), the average controller (IIIcavrg) was able to qualitatively reproduce the behavior. If participant (4), however, would have been part of training, the controller would fit much better to these trials, but also worse for the other two participants which did not bend the lower spine to that extend.The extreme numerical values of the errors on the set of participant (4),Ste4, show, if again comparing Figures 9, 10, and 11, that the error value is/can be used only as a proof of correct response and only as an indicator for possibly incorrect response.


### 4.2. Limits of This Study

First of all, we have to state that an uncertainty quantification (UQ) was not performed in a mathematical rigorous manner.The present study contains several potential sources of errors which could lead to wrong conclusions. Most importantly, the choice of an appropriate model is crucial to analyse the measured data. Moreover, model parameters have to be found and are subject to an uncertain quantification. Further, the experimental data itself is also known to be a typical source of errors.Some sources of errors are easy to mention, for example, the data recording had an error in the marker tracking of approximately 1 cm, which lead to maximum errors in the calculated segment lengths of up to 2 cm. Please refer to AppendixA.2 (invalid trials and exceeder of the “basic” strategy) for a further discussion.Other potential biases are not so trivial to judge. In particular, the level of detail of the model is a crucial first-hand choice. For this study, we chose a double-inverted pendulum (DIP) model. However, simpler single-inverted pendulum (SIP) or more complex models are also used in literature for representing human upright stance balance Günther et al. [30]; Günther and Wagner [31]. Among the more complex models, triple-inverted pendulum (TIP) models are used for lateral stance control, but it exist also models incorporating even a higher complexity, like multiple multibody systems, taking the spinal movement into account, for example, Rupp et al. [40].For the case of simpler models, like the SIP model, they can be eliminated due to two reasons: First, the experimental results of our study show that all participants exhibited roughly the same magnitude of movement in the hip and ankle joint following the perturbation. A SIP would depict only one of these degrees of freedom, totally ignoring the dynamics in the other. Second, the simulated controller IIa (simple PD control) is not reproducing the experimental results and also fails to stabilize, if small delays are introduced (IIb, Figure6). This controller behaves roughly like a SIP with almost all movement occurring in the ankle joint.More interestingly from our perspective would be the case of using more complex models for the analysis. First, currently, some trials had to be excluded from the analysis due to the length criterion (see AppendixA.2) based on anterior/posterior movement. Using a 3D model including 3D joints, the variation of the balance strategy (extension of the basic planar strategy) could have been pictured. By that, the biological system benefits by (i) increasing robustness by further lowering the COM. Furthermore, by using a 3D TIP, the hip and ankle would also move in the A/P direction and the knee joint could be included for even greater increase of robustness too. (ii) The biological system further benefits by a higher flexibility, which is equivalent to a situational dependent change of the balance strategy. Second, all participants exhibited movement in (at least) the lower spine. Some subjects were very stiff in the lower spine region. In some, however, one could see with the naked eye that the lower spine bending contributed more to lateral angular excursions of the trunk than rotations around the hip joints. A model considering a flexible upper segment, in total or in part, would improve the reproduction and depiction of this motion.The last source of uncertainty we want to mention is the commonly known fact that unconstrained nonlinear optimization may end up in several local minima that may or may not be roughly equivalent (see also Results). So better solutions cannot be fully excluded.


### 4.3. Restabilization Strategy

In the present work, we were able to establish and review key factors of balance. Therefore, we want to especially discuss the need for coupling of torques and the role of delays imposed on balance control. We were also able to draw a line of reasoning between our findings and the eigenmovement theory [37], which is discussed together with the predicted delays at the second part of this section.The necessity for torque coupling, where one joint torque needs to take into account the other joint torques in order to account for their changes and influence on the other joints, was proposed and reviewed by earlier studies [9, 37] and was again confirmed by our work. We realize this by sensor signal coupling, utilizing β in (4c) and thus in the generation of motor actions (2). We found that in the latter (2), the most critical parameters are the damping parameters d1 and d2 as they influence torque generation the most (ẏ has the highest influence on the torques, see Figure 12). Especially the ratio between d1 and d2 was found to be roughly d1/d2=1.14±0.03 within our trials and to ensure a certain ratio between the two torques in the ankle and the hip, which was preestablished by the ratios β, variations of one of the d parameters while fixing the other should not exceed 5%. Aside from that critical relation, the possibility of variation of the parameters X is intuitive since a slight change of reaction strength and time (i.e., onset) should not be a problem. If they would, stability in an uncertain system like the upright standing human body would not be possible. Thus, the observed robustness in X is no surprise. Also, the smaller robustness in β, which contains information about the inertia and geometry (which are rather fixed properties of the mechanobiological system), is clear. Already slight change in these ratios results in qualitative changes and instability. Thus, with β and the relation d1/d2, which are together defining the coupling of torques, being the least robust, the necessity of (correct) coupling and relations between torques of different joints becomes obvious. The simulation yielded β12=β21, which states that into each joint torque, the same amount of information about the other joint (angle and angular velocity) is taken into account.Figure 12
ControllerIIIcavrg in the test trial ((5)14200lt). (a) M1 (undelayed MA) full and split into components. (b) M2 (undelayed MH) full and split into components. As can be easily seen, the primary influence on the applied torques always stems from the velocity contributions, while the other two contributions (displacement and perturbation acceleration) just slightly alter the final torques.
(a)
(b)For our motor task, the found parameters ofX and β further impose that y is mainly generated by utilizing ankle angle information (y1: β11 is roughly 5 times larger than β12 and y2: β21 is roughly 4 times larger than β22) and that the torque in the ankle is 5 times larger than in the hip as well (again the ratio between β11 and β21). Numerical values of the fixed relation between torques in different joints were given for A/P balance in [9], with ratios for hip to knee at 1.5 and ankle to knee at 2. Further, it was stated that biarticulated muscles “have the architectural capacity to elicit fixed proportions of two torques about the joints they cross.” According to [9], the delay between ankle and hip could be the hip “waiting” for the ankle reaction to start, in order to react correctly.The predicted delays between kinematic signals and joint torques for controller IIIc (about60ms for ankle and 170ms for hip) are intriguing, because 60ms in particular seems to be a good reflexion of a neuromuscular delay between changes in muscle proprioceptive signals (muscle spindles, Golgi tendon organs) and corresponding changes in motor output, that is, muscle force or joint torque, respectively.Propagation velocities in afferent and efferent, mammalian motor nerves are about100m/s (typ I) or 50m/s (typ II) ([41]). The distance between muscle and motoneuron is roughly 0.5m for human leg muscles and between muscle and brain about 1.5m. The fastest conceivable way of delayed sensor information feedback from mechanical actuator or joint sensors back to the actuator in a biological system is via the monosynaptic pathway muscle-motoneuron-muscle. According to the above-mentioned distances, the total monosynaptic delay would thus be about 10ms (and 30ms in the brain case). Another synaptic delay is added at the neuromuscular junction (between the axon’s end terminal and muscle surface: about 1ms Katz and Miledi [42]). A third contribution to neuromuscular delay seems to be the prevailing one in humans: the electromechanical delay (EMD) as the time passing between the onset of a finite change in electrical muscle surface stimulus beyond the junction, measurable, for example, by a surface electromyography (EMG), to the onset of a measurable, finite change in muscle force ([43]). EMD was found in human’s quadriceps muscle to be 50ms ([44]), when starting from inactivity. Signal processing in neural nets like the motoneuron pool in the spinal cord or in the cerebellum introduces additional contributions to feedback loop delays.Monosynaptic delay in leg muscles via the motoneuron pool is thus about60ms and matches well the delay between joint angle proprioception and ankle torque as found here to be effective for restabilization. This fact is an indication of the neuromuscular design of human beings to be well adapted to the mechanical requirements of regulating human upright posture which is a salient feature in the animal world due to the unique mechanical design of the human foot. A hip joint torque contribution effective for the same requirements is, however, even more delayed by further 110ms. Are there indications that may help to understand this further delay?A delay between distinct subfeatures of motor commands was found to be effective in a different, yet related, motor task: voluntary trunk bending of humans during upright standing ([37, 45]). More precisely, Alexandrov et al. [37] had found for voluntary trunk bending that (i) a superposition of A- and H-eigenmovement is optimal for fulfilling the bending movement and stability at once and that (ii) a delay between the onset of both eigenmovements is essential (A-eigenmovement must start earlier than H-eigenmovement). A- and H-eigenmovements are two motor commands based on differently weighted combinations in joint angular shifts (the eigenvectors of the mechanical, three-link system). In our case, the task is different, somehow an inverse problem: the body has to be restabilized after a disturbance, whereas in the experiments of Alexandrov et al. [45], humans bent their trunks voluntarily. In our case, obviously a delay between the joint torque contributions within a near-A-eigenmovement (i.e., ankle and hip contributions) must become effective (difference between 170 ms and 60 ms).In the first instance, this finding of a delay between local joint torque contributions when compensating a lateral disturbance is clearly different from time sequencing the onset of A- and H-eigenmovements. However, the approach of Alexandrov et al. [37] to motor synthesis based on utilizing the eigenvectors of the mechanical system may be appropriate as a general background with explanatory potential, if particularly combined with ideas ([46]) about the mutual relation of sensory signals and motor commands within the three-dimensional, redundant neuromusculoskeletal system and their potential processing in the central nervous system. Our findings seem to provide some further support for use of the eigenvector approach. Since, the restabilization strategy (2) with (4c) identified in our study is very similar to the A-eigenmovement identified by ((10) in [37]) in trunk bending (“A” refers to the ankle joint). To see this, we have to compare our joint torques Mads and Mhds, respectively, to the first and third components of their vector uA, that is, we first have to compare the ratios of our weightings (Table 3) of ankle and hip angles in (4c) β11/β12≈5 as well as β21/β24≈4—determining Mads and Mhds in (2) by k1·y1 and k2·y2 (k1≈k2≈0.7 s2/deg, see (7)), respectively—to the ratio wA,1/wA,3≈3 of the first and third components in the A-eigenvector ωA in ((10) in [37]). Second, due to k1≈k2, a ratio Mads/Mhds≈β11/β21≈5 of the angle-dependent joint torque contributions in Mads and Mhds for the same angle deflections in ankle and hip compares to the ratio of first and third components uA,1/uA,3≈4.7 in the left equation (11) in [37]. Altogether, we may thus call our identified strategy a “near”-A-eigenmovement because the ratios β11/β12 and β21/β22 of the ankle and hip angle-dependent contributions to the joint torques Mads and Mhds, respectively, in this strategy are slightly distorted as compared to wA,1/wA,3 in the A-eigenmovement (see left equation (10) in [37]), whereas eventually regarding the ratio of the angle-dependent ankle and hip torque contributions themselves (compare β11/β21 and β12/β22 to uA,1/uA,3 in the left equation (11) in [37], there is almost no difference to the A-eigenmovement in terms of the joint torque ratio Mads/Mhds.Velocity-dependent contributions are an implicit part of the eigenvector-based approach in Alexandrov et al. [37]. They are explicitly reflected in the restabilization strategy (2) with (4c) for compensating lateral disturbances as treated here. The ratio of linearly velocity-dependent contributions to Mads and Mhds is approximately d1/d2⋅β11/β21≈5.5 times higher for Mads than for Mhds, which thus makes little difference to the ratio of angle-dependent contributions because d1/d2≈1. With angular velocities that are typically about the same in the hip and the ankle (see Table 1 with maximum angular velocities), the ratio of velocity-dependent contributions to Mads and Mhds is also similar to the ratio of the angle-dependent parts. Remarkably, the velocity-dependent torque contributions (i) dominate the whole restabilization strategy (2) with (4c) by far (compare magnitudes of angle- and velocity-dependent contributions in Figure 12) and (ii) the strategy works even without any “knowledge” of the control system about the disturbance (qi=0) and the angular position (ki=0). Using the word “work” implies that this strategy restabilized the model in the simulations of all trials available, including even those by participant (4) exhibiting exceptional angular amplitudes, albeit the intermittent and final states deviating from corresponding measured ones more than when using angular information (ki≠0). All this indicates that sensor information about rates of mechanical state variables is of fundamental importance for human upright stance. Positional information takes a back seat as does explicit disturbance sensation.As the velocity-dependent torque contributions are time-delayed, they do not represent mechanical friction or “damping” but are rather forecasts (see McMahon [47], p.154–156) of an angle-dependent torque contribution that corresponds to a first-order prediction of a change in joint angle ΔΘj=di/ki⋅Θ̇j (j=a,h) predicted to occur for a given angular velocity Θ̇j that will, however, not take mechanical effect before an instant delayis later. Such velocity-dependent torque contributions compensate for potential deflections in the future that would be much larger than will actually occur as long as their compensatory effect is in fact stabilizing. Here, these torque contributions take mechanical effect delayed by times comparable to the typical time to reach the maximum angular deflections (100…200ms). Introducing motor actions like (2) which are not due to locally instantaneous mechanical effects but based on delayed sensor signals of mechanical variables fed into an information transmission channel to be decoded and transformed into a proportional actuator force later on is information-costly ([48]) but enables access to a greatly enlarged diversity of movement solutions.A rather simple control system based on the system’s mechanical eigenvectors and one or another tailored time delay can be evidently set up for initiating voluntary movements or compensating significant disturbances of anN-link-inverted pendulum without escalating the movement system into resonance disaster. It should be checked whether a recently suggested stabilization strategy for a TIP model simulating quiet human stance in the sagittal plane ([31]) is in fact similar, related, or even more or less identical to the eigenvector approach. It seems an obvious thing to take further steps and broaden the check for how general this idea of motor control based on mechanical eigenvectors is. Since, to the best of our knowledge and somehow surprisingly, this approach has solely been explicitly applied to human stance [9, 10, 13, 49–51] in the field of model-based, biomechanical movement synthesis of multilink systems so far. We cannot see obvious reasons, however, why cyclic or goal-directed movement should not be a promising field alike for eigenvector-based investigations to their motor control.


## 4.1. Generality of the Controller

The controller IIIc is not only used to reproduce trials (to which it was parametrized to), but also to predict trials not used in training (remember training error of 1.35deg/ts). The quality of prediction for this test set is absolutely within the range of the fitting to the training set (even if the numerical error value on the test set with 2.24 deg/ts or 1.7 deg/ts is higher than in training, the fits are very good as can be seen, comparing Figures 9 and 10).Notable is that the training consisted of trials from two participants, which were used regardless of their magnitude, over the full range of perturbations tested in the experiment, also ignoring the direction. As this was not only possible, but also leading to good results, the strong correlation between trials of different magnitudes and participants could be verified: the same control scheme, with one parameter set could be used. The full range of reactions to the perturbation, tolerable by the controller, is also illustrated by Figure5. We are certain that for a greater number of trials, participants, and magnitudes, this mean result would converge to a smaller error than possible now, under the prior mentioned premises, and to a certainly much more accurate prediction if only one magnitude of perturbations is used with many participants. It would also be possible to extend the found control law such that it considers the perturbation strength to a higher degree and apply different magnitudes of reaction accordingly.However, we have chosen this advance to prove some generality of the controller (e.g., robustness over a broad range of trials). Even for participant (4), where a variation of strategy could be observed (much greater angles in the hip due to extended spinal bending in l1-l5), the average controller (IIIcavrg) was able to qualitatively reproduce the behavior. If participant (4), however, would have been part of training, the controller would fit much better to these trials, but also worse for the other two participants which did not bend the lower spine to that extend.The extreme numerical values of the errors on the set of participant (4),Ste4, show, if again comparing Figures 9, 10, and 11, that the error value is/can be used only as a proof of correct response and only as an indicator for possibly incorrect response.


## 4.2. Limits of This Study

First of all, we have to state that an uncertainty quantification (UQ) was not performed in a mathematical rigorous manner.The present study contains several potential sources of errors which could lead to wrong conclusions. Most importantly, the choice of an appropriate model is crucial to analyse the measured data. Moreover, model parameters have to be found and are subject to an uncertain quantification. Further, the experimental data itself is also known to be a typical source of errors.Some sources of errors are easy to mention, for example, the data recording had an error in the marker tracking of approximately 1 cm, which lead to maximum errors in the calculated segment lengths of up to 2 cm. Please refer to AppendixA.2 (invalid trials and exceeder of the “basic” strategy) for a further discussion.Other potential biases are not so trivial to judge. In particular, the level of detail of the model is a crucial first-hand choice. For this study, we chose a double-inverted pendulum (DIP) model. However, simpler single-inverted pendulum (SIP) or more complex models are also used in literature for representing human upright stance balance Günther et al. [30]; Günther and Wagner [31]. Among the more complex models, triple-inverted pendulum (TIP) models are used for lateral stance control, but it exist also models incorporating even a higher complexity, like multiple multibody systems, taking the spinal movement into account, for example, Rupp et al. [40].For the case of simpler models, like the SIP model, they can be eliminated due to two reasons: First, the experimental results of our study show that all participants exhibited roughly the same magnitude of movement in the hip and ankle joint following the perturbation. A SIP would depict only one of these degrees of freedom, totally ignoring the dynamics in the other. Second, the simulated controller IIa (simple PD control) is not reproducing the experimental results and also fails to stabilize, if small delays are introduced (IIb, Figure6). This controller behaves roughly like a SIP with almost all movement occurring in the ankle joint.More interestingly from our perspective would be the case of using more complex models for the analysis. First, currently, some trials had to be excluded from the analysis due to the length criterion (see AppendixA.2) based on anterior/posterior movement. Using a 3D model including 3D joints, the variation of the balance strategy (extension of the basic planar strategy) could have been pictured. By that, the biological system benefits by (i) increasing robustness by further lowering the COM. Furthermore, by using a 3D TIP, the hip and ankle would also move in the A/P direction and the knee joint could be included for even greater increase of robustness too. (ii) The biological system further benefits by a higher flexibility, which is equivalent to a situational dependent change of the balance strategy. Second, all participants exhibited movement in (at least) the lower spine. Some subjects were very stiff in the lower spine region. In some, however, one could see with the naked eye that the lower spine bending contributed more to lateral angular excursions of the trunk than rotations around the hip joints. A model considering a flexible upper segment, in total or in part, would improve the reproduction and depiction of this motion.The last source of uncertainty we want to mention is the commonly known fact that unconstrained nonlinear optimization may end up in several local minima that may or may not be roughly equivalent (see also Results). So better solutions cannot be fully excluded.


## 4.3. Restabilization Strategy

In the present work, we were able to establish and review key factors of balance. Therefore, we want to especially discuss the need for coupling of torques and the role of delays imposed on balance control. We were also able to draw a line of reasoning between our findings and the eigenmovement theory [37], which is discussed together with the predicted delays at the second part of this section.The necessity for torque coupling, where one joint torque needs to take into account the other joint torques in order to account for their changes and influence on the other joints, was proposed and reviewed by earlier studies [9, 37] and was again confirmed by our work. We realize this by sensor signal coupling, utilizing β in (4c) and thus in the generation of motor actions (2). We found that in the latter (2), the most critical parameters are the damping parameters d1 and d2 as they influence torque generation the most (ẏ has the highest influence on the torques, see Figure 12). Especially the ratio between d1 and d2 was found to be roughly d1/d2=1.14±0.03 within our trials and to ensure a certain ratio between the two torques in the ankle and the hip, which was preestablished by the ratios β, variations of one of the d parameters while fixing the other should not exceed 5%. Aside from that critical relation, the possibility of variation of the parameters X is intuitive since a slight change of reaction strength and time (i.e., onset) should not be a problem. If they would, stability in an uncertain system like the upright standing human body would not be possible. Thus, the observed robustness in X is no surprise. Also, the smaller robustness in β, which contains information about the inertia and geometry (which are rather fixed properties of the mechanobiological system), is clear. Already slight change in these ratios results in qualitative changes and instability. Thus, with β and the relation d1/d2, which are together defining the coupling of torques, being the least robust, the necessity of (correct) coupling and relations between torques of different joints becomes obvious. The simulation yielded β12=β21, which states that into each joint torque, the same amount of information about the other joint (angle and angular velocity) is taken into account.Figure 12
ControllerIIIcavrg in the test trial ((5)14200lt). (a) M1 (undelayed MA) full and split into components. (b) M2 (undelayed MH) full and split into components. As can be easily seen, the primary influence on the applied torques always stems from the velocity contributions, while the other two contributions (displacement and perturbation acceleration) just slightly alter the final torques.
(a)
(b)For our motor task, the found parameters ofX and β further impose that y is mainly generated by utilizing ankle angle information (y1: β11 is roughly 5 times larger than β12 and y2: β21 is roughly 4 times larger than β22) and that the torque in the ankle is 5 times larger than in the hip as well (again the ratio between β11 and β21). Numerical values of the fixed relation between torques in different joints were given for A/P balance in [9], with ratios for hip to knee at 1.5 and ankle to knee at 2. Further, it was stated that biarticulated muscles “have the architectural capacity to elicit fixed proportions of two torques about the joints they cross.” According to [9], the delay between ankle and hip could be the hip “waiting” for the ankle reaction to start, in order to react correctly.The predicted delays between kinematic signals and joint torques for controller IIIc (about60ms for ankle and 170ms for hip) are intriguing, because 60ms in particular seems to be a good reflexion of a neuromuscular delay between changes in muscle proprioceptive signals (muscle spindles, Golgi tendon organs) and corresponding changes in motor output, that is, muscle force or joint torque, respectively.Propagation velocities in afferent and efferent, mammalian motor nerves are about100m/s (typ I) or 50m/s (typ II) ([41]). The distance between muscle and motoneuron is roughly 0.5m for human leg muscles and between muscle and brain about 1.5m. The fastest conceivable way of delayed sensor information feedback from mechanical actuator or joint sensors back to the actuator in a biological system is via the monosynaptic pathway muscle-motoneuron-muscle. According to the above-mentioned distances, the total monosynaptic delay would thus be about 10ms (and 30ms in the brain case). Another synaptic delay is added at the neuromuscular junction (between the axon’s end terminal and muscle surface: about 1ms Katz and Miledi [42]). A third contribution to neuromuscular delay seems to be the prevailing one in humans: the electromechanical delay (EMD) as the time passing between the onset of a finite change in electrical muscle surface stimulus beyond the junction, measurable, for example, by a surface electromyography (EMG), to the onset of a measurable, finite change in muscle force ([43]). EMD was found in human’s quadriceps muscle to be 50ms ([44]), when starting from inactivity. Signal processing in neural nets like the motoneuron pool in the spinal cord or in the cerebellum introduces additional contributions to feedback loop delays.Monosynaptic delay in leg muscles via the motoneuron pool is thus about60ms and matches well the delay between joint angle proprioception and ankle torque as found here to be effective for restabilization. This fact is an indication of the neuromuscular design of human beings to be well adapted to the mechanical requirements of regulating human upright posture which is a salient feature in the animal world due to the unique mechanical design of the human foot. A hip joint torque contribution effective for the same requirements is, however, even more delayed by further 110ms. Are there indications that may help to understand this further delay?A delay between distinct subfeatures of motor commands was found to be effective in a different, yet related, motor task: voluntary trunk bending of humans during upright standing ([37, 45]). More precisely, Alexandrov et al. [37] had found for voluntary trunk bending that (i) a superposition of A- and H-eigenmovement is optimal for fulfilling the bending movement and stability at once and that (ii) a delay between the onset of both eigenmovements is essential (A-eigenmovement must start earlier than H-eigenmovement). A- and H-eigenmovements are two motor commands based on differently weighted combinations in joint angular shifts (the eigenvectors of the mechanical, three-link system). In our case, the task is different, somehow an inverse problem: the body has to be restabilized after a disturbance, whereas in the experiments of Alexandrov et al. [45], humans bent their trunks voluntarily. In our case, obviously a delay between the joint torque contributions within a near-A-eigenmovement (i.e., ankle and hip contributions) must become effective (difference between 170 ms and 60 ms).In the first instance, this finding of a delay between local joint torque contributions when compensating a lateral disturbance is clearly different from time sequencing the onset of A- and H-eigenmovements. However, the approach of Alexandrov et al. [37] to motor synthesis based on utilizing the eigenvectors of the mechanical system may be appropriate as a general background with explanatory potential, if particularly combined with ideas ([46]) about the mutual relation of sensory signals and motor commands within the three-dimensional, redundant neuromusculoskeletal system and their potential processing in the central nervous system. Our findings seem to provide some further support for use of the eigenvector approach. Since, the restabilization strategy (2) with (4c) identified in our study is very similar to the A-eigenmovement identified by ((10) in [37]) in trunk bending (“A” refers to the ankle joint). To see this, we have to compare our joint torques Mads and Mhds, respectively, to the first and third components of their vector uA, that is, we first have to compare the ratios of our weightings (Table 3) of ankle and hip angles in (4c) β11/β12≈5 as well as β21/β24≈4—determining Mads and Mhds in (2) by k1·y1 and k2·y2 (k1≈k2≈0.7 s2/deg, see (7)), respectively—to the ratio wA,1/wA,3≈3 of the first and third components in the A-eigenvector ωA in ((10) in [37]). Second, due to k1≈k2, a ratio Mads/Mhds≈β11/β21≈5 of the angle-dependent joint torque contributions in Mads and Mhds for the same angle deflections in ankle and hip compares to the ratio of first and third components uA,1/uA,3≈4.7 in the left equation (11) in [37]. Altogether, we may thus call our identified strategy a “near”-A-eigenmovement because the ratios β11/β12 and β21/β22 of the ankle and hip angle-dependent contributions to the joint torques Mads and Mhds, respectively, in this strategy are slightly distorted as compared to wA,1/wA,3 in the A-eigenmovement (see left equation (10) in [37]), whereas eventually regarding the ratio of the angle-dependent ankle and hip torque contributions themselves (compare β11/β21 and β12/β22 to uA,1/uA,3 in the left equation (11) in [37], there is almost no difference to the A-eigenmovement in terms of the joint torque ratio Mads/Mhds.Velocity-dependent contributions are an implicit part of the eigenvector-based approach in Alexandrov et al. [37]. They are explicitly reflected in the restabilization strategy (2) with (4c) for compensating lateral disturbances as treated here. The ratio of linearly velocity-dependent contributions to Mads and Mhds is approximately d1/d2⋅β11/β21≈5.5 times higher for Mads than for Mhds, which thus makes little difference to the ratio of angle-dependent contributions because d1/d2≈1. With angular velocities that are typically about the same in the hip and the ankle (see Table 1 with maximum angular velocities), the ratio of velocity-dependent contributions to Mads and Mhds is also similar to the ratio of the angle-dependent parts. Remarkably, the velocity-dependent torque contributions (i) dominate the whole restabilization strategy (2) with (4c) by far (compare magnitudes of angle- and velocity-dependent contributions in Figure 12) and (ii) the strategy works even without any “knowledge” of the control system about the disturbance (qi=0) and the angular position (ki=0). Using the word “work” implies that this strategy restabilized the model in the simulations of all trials available, including even those by participant (4) exhibiting exceptional angular amplitudes, albeit the intermittent and final states deviating from corresponding measured ones more than when using angular information (ki≠0). All this indicates that sensor information about rates of mechanical state variables is of fundamental importance for human upright stance. Positional information takes a back seat as does explicit disturbance sensation.As the velocity-dependent torque contributions are time-delayed, they do not represent mechanical friction or “damping” but are rather forecasts (see McMahon [47], p.154–156) of an angle-dependent torque contribution that corresponds to a first-order prediction of a change in joint angle ΔΘj=di/ki⋅Θ̇j (j=a,h) predicted to occur for a given angular velocity Θ̇j that will, however, not take mechanical effect before an instant delayis later. Such velocity-dependent torque contributions compensate for potential deflections in the future that would be much larger than will actually occur as long as their compensatory effect is in fact stabilizing. Here, these torque contributions take mechanical effect delayed by times comparable to the typical time to reach the maximum angular deflections (100…200ms). Introducing motor actions like (2) which are not due to locally instantaneous mechanical effects but based on delayed sensor signals of mechanical variables fed into an information transmission channel to be decoded and transformed into a proportional actuator force later on is information-costly ([48]) but enables access to a greatly enlarged diversity of movement solutions.A rather simple control system based on the system’s mechanical eigenvectors and one or another tailored time delay can be evidently set up for initiating voluntary movements or compensating significant disturbances of anN-link-inverted pendulum without escalating the movement system into resonance disaster. It should be checked whether a recently suggested stabilization strategy for a TIP model simulating quiet human stance in the sagittal plane ([31]) is in fact similar, related, or even more or less identical to the eigenvector approach. It seems an obvious thing to take further steps and broaden the check for how general this idea of motor control based on mechanical eigenvectors is. Since, to the best of our knowledge and somehow surprisingly, this approach has solely been explicitly applied to human stance [9, 10, 13, 49–51] in the field of model-based, biomechanical movement synthesis of multilink systems so far. We cannot see obvious reasons, however, why cyclic or goal-directed movement should not be a promising field alike for eigenvector-based investigations to their motor control.



---

*Source: 4767624-2018-05-02.xml*